About special functions reogranization

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

About special functions reogranization

Michele Ginesi
Dear all,
I'm Michele Ginesi, the student working on special functions for GSoC2017.
During the last period I saw on this mailing list some discussions about
the reorganization of special functions in Octave. I saw there is the
idea to use the standard math library of C++ or GSL library.
Unfortunately, as I explained in the last post on my blog [1], I think
it is not the case to fully substitute the old (but still the more
accurate) Fortran code present in the Amos library.
For example, for what concern Bessel functions, I would  suggest to
"unlock" the Amos, since the bug related to them [2] is caused by the
fact that if the magnitude of the input is too large, the output is
simply not computed (you can find more details on [1], as well as other
thoughts relative to incomplete gamma and beta functions).
I hope you can give me some opinions and/or suggestions.
Thank you for your kind attention.

[1] http://gsocspecfun.blogspot.it/2017/07/gnu-scientific-library.html
[2] https://savannah.gnu.org/bugs/?func=detailitem&item_id=48316

--
Michele Ginesi

Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: About special functions reogranization

Rik-4
On 07/13/2017 09:00 AM, [hidden email] wrote:
Subject:
About special functions reogranization
From:
Michele [hidden email]
Date:
07/13/2017 01:49 AM
To:
[hidden email]
List-Post:
[hidden email]
Content-Transfer-Encoding:
7bit
Precedence:
list
MIME-Version:
1.0
Message-ID:
[hidden email]
Content-Type:
text/plain; charset=utf-8; format=flowed
Message:
3

Dear all,
I'm Michele Ginesi, the student working on special functions for GSoC2017.
During the last period I saw on this mailing list some discussions about the reorganization of special functions in Octave. I saw there is the idea to use the standard math library of C++ or GSL library.
Unfortunately, as I explained in the last post on my blog [1], I think it is not the case to fully substitute the old (but still the more accurate) Fortran code present in the Amos library.
For example, for what concern Bessel functions, I would  suggest to "unlock" the Amos, since the bug related to them [2] is caused by the fact that if the magnitude of the input is too large, the output is simply not computed (you can find more details on [1], as well as other thoughts relative to incomplete gamma and beta functions).
I hope you can give me some opinions and/or suggestions.
Thank you for your kind attention.

[1] http://gsocspecfun.blogspot.it/2017/07/gnu-scientific-library.html
[2] https://savannah.gnu.org/bugs/?func=detailitem&item_id=48316


Michele,

I read your summary of the state of the bessel, gammainc, and betainc functions with interest.

For the Bessel functions, you state that Matlab allows a complex alpha parameter, but this is not the case (http://www.mathworks.com/help/matlab/ref/besselj.html).  This doesn't detract much from your argument.  The real issue which prevents using the GSL library is that the input itself must be real, while AMOS and Matlab both support a complex input Z.

For gammainc, you assert that there is no "scaled" option.  But, there is an unnormalized version of the function gsl_sf_gamma_inc.  It wouldn't be too difficult to transform the unnormalized output to the scaled output defined by Matlab.  It does involve multiplication by a scaling factor, but underflow potential, which is what the scaled option in Matlab is seeking to prevent, looks much more unlikely in this case.  However, you're still probably right to create your own implementation given that the GSL library requires X to be positive for all versions of the function.  The current implementation of gammainc is in Fortran.  Will you be writing in Fortran, C++, or Octave?  I doubt you can get the performance out of an Octave m-file that you can out of the other two compiled languages, unless there are very few iterative loops in the algorithm.

I agree with your conclusion that betainc is insufficient.  You mention writing an m-file to replace what is now Fortran code.  Again, my concern here is that the performance will not be high enough unless the algorithm does not use iterative loops.

Happy coding,
Rik
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: About special functions reogranization

Michele Ginesi
On 07/14/2017 09:04 PM, Rik wrote:
On 07/13/2017 09:00 AM, [hidden email] wrote:
Subject:
About special functions reogranization
From:
Michele [hidden email]
Date:
07/13/2017 01:49 AM
To:
[hidden email]
List-Post:
[hidden email]
Content-Transfer-Encoding:
7bit
Precedence:
list
MIME-Version:
1.0
Message-ID:
[hidden email]
Content-Type:
text/plain; charset=utf-8; format=flowed
Message:
3

Dear all,
I'm Michele Ginesi, the student working on special functions for GSoC2017.
During the last period I saw on this mailing list some discussions about the reorganization of special functions in Octave. I saw there is the idea to use the standard math library of C++ or GSL library.
Unfortunately, as I explained in the last post on my blog [1], I think it is not the case to fully substitute the old (but still the more accurate) Fortran code present in the Amos library.
For example, for what concern Bessel functions, I would  suggest to "unlock" the Amos, since the bug related to them [2] is caused by the fact that if the magnitude of the input is too large, the output is simply not computed (you can find more details on [1], as well as other thoughts relative to incomplete gamma and beta functions).
I hope you can give me some opinions and/or suggestions.
Thank you for your kind attention.

[1] http://gsocspecfun.blogspot.it/2017/07/gnu-scientific-library.html
[2] https://savannah.gnu.org/bugs/?func=detailitem&item_id=48316


Michele,

I read your summary of the state of the bessel, gammainc, and betainc functions with interest.

For the Bessel functions, you state that Matlab allows a complex alpha parameter, but this is not the case (http://www.mathworks.com/help/matlab/ref/besselj.html).  This doesn't detract much from your argument.  The real issue which prevents using the GSL library is that the input itself must be real, while AMOS and Matlab both support a complex input Z.

For gammainc, you assert that there is no "scaled" option.  But, there is an unnormalized version of the function gsl_sf_gamma_inc.  It wouldn't be too difficult to transform the unnormalized output to the scaled output defined by Matlab.  It does involve multiplication by a scaling factor, but underflow potential, which is what the scaled option in Matlab is seeking to prevent, looks much more unlikely in this case.  However, you're still probably right to create your own implementation given that the GSL library requires X to be positive for all versions of the function.  The current implementation of gammainc is in Fortran.  Will you be writing in Fortran, C++, or Octave?  I doubt you can get the performance out of an Octave m-file that you can out of the other two compiled languages, unless there are very few iterative loops in the algorithm.

I agree with your conclusion that betainc is insufficient.  You mention writing an m-file to replace what is now Fortran code.  Again, my concern here is that the performance will not be high enough unless the algorithm does not use iterative loops.

Happy coding,
Rik

Dear Rik,

thanks for your comments. I already worked on gammainc, you can see the details in [1]. The important thing is that Marco, Nir and I implemented it as an .m file, with a subroutine in C++ (you can find this implementation in the bookmark "gammainc" of my repository [2]). Every special function has more or less the same structure: the domain is split into subdomains and different subfunctions are called. Each subfunction implements an algorithm, the most used ones being series expansions and continued fractions. I think that I could continue my work (e.g. betainc) in this way, finding algorithms that are accurate in their domains. Then, with maybe the help of the community, start to translate the subfunctions from .m files to C++ (I'm not autonomous with C++). In this way, at the end of GSoC we will have accurate algorithms, that in the next time can become also fast.

[1] http://gsocspecfun.blogspot.it/2017/06/gammainc.html
[2] https://bitbucket.org/M_Ginesi/octave

--
Michele Ginesi
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: About special functions reogranization

Robert T. Short
On 07/19/2017 01:20 AM, Michele wrote:
On 07/14/2017 09:04 PM, Rik wrote:
On 07/13/2017 09:00 AM, [hidden email] wrote:
Subject:
About special functions reogranization
From:
Michele [hidden email]
Date:
07/13/2017 01:49 AM
To:
[hidden email]
List-Post:
[hidden email]
Content-Transfer-Encoding:
7bit
Precedence:
list
MIME-Version:
1.0
Message-ID:
[hidden email]
Content-Type:
text/plain; charset=utf-8; format=flowed
Message:
3

Dear all,
I'm Michele Ginesi, the student working on special functions for GSoC2017.
During the last period I saw on this mailing list some discussions about the reorganization of special functions in Octave. I saw there is the idea to use the standard math library of C++ or GSL library.
Unfortunately, as I explained in the last post on my blog [1], I think it is not the case to fully substitute the old (but still the more accurate) Fortran code present in the Amos library.
For example, for what concern Bessel functions, I would  suggest to "unlock" the Amos, since the bug related to them [2] is caused by the fact that if the magnitude of the input is too large, the output is simply not computed (you can find more details on [1], as well as other thoughts relative to incomplete gamma and beta functions).
I hope you can give me some opinions and/or suggestions.
Thank you for your kind attention.

[1] http://gsocspecfun.blogspot.it/2017/07/gnu-scientific-library.html
[2] https://savannah.gnu.org/bugs/?func=detailitem&item_id=48316


Michele,

I read your summary of the state of the bessel, gammainc, and betainc functions with interest.

For the Bessel functions, you state that Matlab allows a complex alpha parameter, but this is not the case (http://www.mathworks.com/help/matlab/ref/besselj.html).  This doesn't detract much from your argument.  The real issue which prevents using the GSL library is that the input itself must be real, while AMOS and Matlab both support a complex input Z.

For gammainc, you assert that there is no "scaled" option.  But, there is an unnormalized version of the function gsl_sf_gamma_inc.  It wouldn't be too difficult to transform the unnormalized output to the scaled output defined by Matlab.  It does involve multiplication by a scaling factor, but underflow potential, which is what the scaled option in Matlab is seeking to prevent, looks much more unlikely in this case.  However, you're still probably right to create your own implementation given that the GSL library requires X to be positive for all versions of the function.  The current implementation of gammainc is in Fortran.  Will you be writing in Fortran, C++, or Octave?  I doubt you can get the performance out of an Octave m-file that you can out of the other two compiled languages, unless there are very few iterative loops in the algorithm.

I agree with your conclusion that betainc is insufficient.  You mention writing an m-file to replace what is now Fortran code.  Again, my concern here is that the performance will not be high enough unless the algorithm does not use iterative loops.

Happy coding,
Rik

Dear Rik,

thanks for your comments. I already worked on gammainc, you can see the details in [1]. The important thing is that Marco, Nir and I implemented it as an .m file, with a subroutine in C++ (you can find this implementation in the bookmark "gammainc" of my repository [2]). Every special function has more or less the same structure: the domain is split into subdomains and different subfunctions are called. Each subfunction implements an algorithm, the most used ones being series expansions and continued fractions. I think that I could continue my work (e.g. betainc) in this way, finding algorithms that are accurate in their domains. Then, with maybe the help of the community, start to translate the subfunctions from .m files to C++ (I'm not autonomous with C++). In this way, at the end of GSoC we will have accurate algorithms, that in the next time can become also fast.

[1] http://gsocspecfun.blogspot.it/2017/06/gammainc.html
[2] https://bitbucket.org/M_Ginesi/octave

--
Michele Ginesi

Just FYI, attached is a list of MATLAB special functions and the arguments and parameters that they take.  It is possible that there is an error in the table, but I took some time to put this together so there shouldn't be many errors.

Bob


ListOfSpecialFunctions.pdf (38K) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: About special functions reogranization

Rik-4
In reply to this post by Michele Ginesi
On 07/19/2017 01:20 AM, Michele wrote:
On 07/14/2017 09:04 PM, Rik wrote:
On 07/13/2017 09:00 AM, [hidden email] wrote:
Subject:
About special functions reogranization
From:
Michele [hidden email]
Date:
07/13/2017 01:49 AM
To:
[hidden email]
List-Post:
[hidden email]
Content-Transfer-Encoding:
7bit
Precedence:
list
MIME-Version:
1.0
Message-ID:
[hidden email]
Content-Type:
text/plain; charset=utf-8; format=flowed
Message:
3

Dear all,
I'm Michele Ginesi, the student working on special functions for GSoC2017.
During the last period I saw on this mailing list some discussions about the reorganization of special functions in Octave. I saw there is the idea to use the standard math library of C++ or GSL library.
Unfortunately, as I explained in the last post on my blog [1], I think it is not the case to fully substitute the old (but still the more accurate) Fortran code present in the Amos library.
For example, for what concern Bessel functions, I would  suggest to "unlock" the Amos, since the bug related to them [2] is caused by the fact that if the magnitude of the input is too large, the output is simply not computed (you can find more details on [1], as well as other thoughts relative to incomplete gamma and beta functions).
I hope you can give me some opinions and/or suggestions.
Thank you for your kind attention.

[1] http://gsocspecfun.blogspot.it/2017/07/gnu-scientific-library.html
[2] https://savannah.gnu.org/bugs/?func=detailitem&item_id=48316


Michele,

I read your summary of the state of the bessel, gammainc, and betainc functions with interest.

For the Bessel functions, you state that Matlab allows a complex alpha parameter, but this is not the case (http://www.mathworks.com/help/matlab/ref/besselj.html).  This doesn't detract much from your argument.  The real issue which prevents using the GSL library is that the input itself must be real, while AMOS and Matlab both support a complex input Z.

For gammainc, you assert that there is no "scaled" option.  But, there is an unnormalized version of the function gsl_sf_gamma_inc.  It wouldn't be too difficult to transform the unnormalized output to the scaled output defined by Matlab.  It does involve multiplication by a scaling factor, but underflow potential, which is what the scaled option in Matlab is seeking to prevent, looks much more unlikely in this case.  However, you're still probably right to create your own implementation given that the GSL library requires X to be positive for all versions of the function.  The current implementation of gammainc is in Fortran.  Will you be writing in Fortran, C++, or Octave?  I doubt you can get the performance out of an Octave m-file that you can out of the other two compiled languages, unless there are very few iterative loops in the algorithm.

I agree with your conclusion that betainc is insufficient.  You mention writing an m-file to replace what is now Fortran code.  Again, my concern here is that the performance will not be high enough unless the algorithm does not use iterative loops.

Happy coding,
Rik

Dear Rik,

thanks for your comments. I already worked on gammainc, you can see the details in [1]. The important thing is that Marco, Nir and I implemented it as an .m file, with a subroutine in C++ (you can find this implementation in the bookmark "gammainc" of my repository [2]). Every special function has more or less the same structure: the domain is split into subdomains and different subfunctions are called. Each subfunction implements an algorithm, the most used ones being series expansions and continued fractions. I think that I could continue my work (e.g. betainc) in this way, finding algorithms that are accurate in their domains. Then, with maybe the help of the community, start to translate the subfunctions from .m files to C++ (I'm not autonomous with C++). In this way, at the end of GSoC we will have accurate algorithms, that in the next time can become also fast.

This is a good strategy.  Concentrate on correctness first, and performance second.  It is usually not too difficult to translate an algorithm, once it has been written, from one programming language to another.  And the C++ used by Octave has a lot of features of the m-file language itself which will make the transition even easier.

--Rik
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: About special functions reogranization

NJank
In reply to this post by Robert T. Short
On Jul 19, 2017 10:28 AM, "Robert T. Short" 

Just FYI, attached is a list of MATLAB special functions and the arguments and parameters that they take.  It is possible that there is an error in the table, but I took some time to put this together so there shouldn't be many errors.

Bob

Nice list. Matlab recently changed how some of their functions handle NaN and empty inputs (scalar and n-dimensional for both). It would be worth checking on the input/output of those for compatibility. 
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: About special functions reogranization

Robert T. Short
On 07/19/2017 08:14 AM, Nicholas Jankowski wrote:
On Jul 19, 2017 10:28 AM, "Robert T. Short" 

Just FYI, attached is a list of MATLAB special functions and the arguments and parameters that they take.  It is possible that there is an error in the table, but I took some time to put this together so there shouldn't be many errors.

Bob

Nice list. Matlab recently changed how some of their functions handle NaN and empty inputs (scalar and n-dimensional for both). It would be worth checking on the input/output of those for compatibility. 

I will look at it.  Thanks.  I have 2017.1 available and we should have .2 soon.  This list was done with 2014.

BTW, I agree that correctness is more important than performance.  I have been looking at algorithms for Marcum Q functions are related.  There are methods that purport to be extremely accurate but sloooooow, so you can take this to extreme.  Tests are also critical.

Bob

Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: About special functions reogranization

NJank
And I CC'd the wrong list... It's been that kind of a day.

On Jul 19, 2017 2:53 PM, "Nicholas Jankowski" <[hidden email]> wrote:
Forgot to cc the list:

On Jul 19, 2017 2:52 PM, Nicholas Jankowski <[hidden email]> wrote:
Ok, thanks. I'm sans MATLAB for a little while or I'd help throw a few compat. tests together. basically I've been throwing Nan and [] at each input, then combinations of them for multi input functions, then a couple n-dimensional input forms ( nan(2,2), ones(2,0,1), etc) I think I had stumbled into their general handling scheme for the empties when I redid the statistics functions [1]. Don't know how much of that carries through for other functions.


Loading...