New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DOC: update RESTORE tutorial to use new noise estimation technique #463

Merged
merged 1 commit into from Nov 17, 2014

Conversation

Projects
None yet
5 participants
@arokem
Member

arokem commented Nov 13, 2014

No description provided.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Nov 17, 2014

Have you checked the output of the tutorial? Is it of the same quality as before?

@arokem

This comment has been minimized.

Member

arokem commented Nov 17, 2014

Yes. The noise estimates with this method are slightly smaller than my
previous estimate, but it seems to have a similar ultimate effect on the
ellipsoid maps and on the FA histograms.

On Mon, Nov 17, 2014 at 8:03 AM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Have you checked the output of the tutorial? Is it of the same quality as
before?


Reply to this email directly or view it on GitHub
#463 (comment).

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Nov 17, 2014

Okay good to hear.

Garyfallidis added a commit that referenced this pull request Nov 17, 2014

Merge pull request #463 from arokem/restore-noise
DOC: update RESTORE tutorial to use new noise estimation technique

@Garyfallidis Garyfallidis merged commit ed6838d into nipy:master Nov 17, 2014

1 check passed

continuous-integration/travis-ci The Travis CI build passed
Details
@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 17, 2014

A reason why it might be smaller is because no correction is made for the noise type. A small example to illustrate my point

import numpy as np
from dipy.denoise.noise_estimate import estimate_sigma
n1 = np.random.randn(100,100,100) * 50 + 25 # std is 50, mean is 25
n2 = np.random.randn(100,100,100) * 50 + 25
noise = np.sqrt(n1**2 + n2**2) # rician noise

sigma = estimate_sigma(noise)

# std was really 50 in fact 
print(sigma[0])
print(noise.std())

# apply correction factor based on mean, std and type of noise 
np.sqrt(noise.var() / _xi(25,50,1)) # mean=25, std=50, N=1
# answer is 52.323814196283301, much closer to 50

and this _xi function uses scipy.special.hyp1f1, which is numerically unstable in scipy, with various fixes in various PR for like forever.

from scipy.special import hyp1f1
eta=800.
sigma=20.
N=1
hyp1f1(-0.5, N, -eta**2/(2*sigma**2))

will return nan, and it's easy to have in the b0 a value of 800 with low noise.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 25, 2014

Maybe we should collect a list of the hyp1f1 fixes to see if we could improve it with a small amount of effort?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 25, 2014

I already have a custom fork with those fixes merged https://github.com/samuelstjean/scipy
I was using that before, it doesn't give nans where the current scipy does (yay!), but I found scipy 0.14 to be unstable by something like 1e14 for some large negative values (no idea about the fork, I stopped using it since apparently nobody wants to compile and install a custom scipy).

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 25, 2014

How much work do you think it would be to get those fixes merged upstream?

Is it practical to carry our own fixed copy temporarily?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 25, 2014

No idea, those fixes are all in pending PRs already, I just merged them myself for convenience, there is always a lot of PR in scipy, so maybe they just got forgotten. The issue lies in fortran code, so one would need to compile some functions to use these fixes.

While mac and linux users are probably covered by the setup.py without any effort, probably nobody on windows has the required libraries/compilers installed (by default).

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 25, 2014

Would it be practical to implement our own using Cython do you think?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 25, 2014

That would be impossible, see #468 as to why it's a numerical nightmare. For the mathematicians amongst us, see instead https://en.wikipedia.org/wiki/Confluent_hypergeometric_function.

That being said, I am currently using for my own needs the implementation in the gsl wrapped in cython. While it is super easy to use, it's also licensed under the gpl and probably can't be used, but then again cvxopt is also gpl (http://cvxopt.org/copyright.html) and also already used (although I just realized it was gpl while looking for suitesparse, and that's probably not permitted...).

Anyway, suitesparse is a wrapper under the BSD license for gpl code (https://github.com/njsmith/scikits-sparse), and they just say respect whatever license, so maybe we could do a legal trick like that.

@matthew-brett

This comment has been minimized.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 25, 2014

Well, it is a direct implementation of the function, but I highly doubt it holds numerically. I ran the code in octave, and just looking at the comments and using an unstable provided example, I have

octave gsl version
octave-cli:23> hyperg_1F1(-4.5, 1, -75)
ans = 6767147.04168941

wolfram alpha
confluent hypergeometric function of the first kind, a=-4.5, b=1, x=-75
6.76715×10^6

kummer.m function version
octave:16> kummer(-4.5,1,-75)
ans = 14380208.6642492

So... it's probably gonna blow up in a bunch of places.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 25, 2014

Ah yes, I didn't see how simple that code was.

How about porting the Fortran from your patched version of scipy to Cython? Is that feasible?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 25, 2014

I don't know anything/enough about fortran to do that, it's mostly numerical stability fixes involving gamma functions and their variants apparently according to a quick diff
samuelstjean/scipy@gertingold:master...master

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 25, 2014

I'm guessing you've merged scipy/scipy#3022 ?

I noticed the comment on that PR that there's a lot of work to do to make the test_mpmath tests pass.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 25, 2014

Seems like the only one actually, there are a bunch of bugs reports not PR after all. Anyway, I will try to check if this fixes those errors https://github.com/scipy/scipy/issues?q=hyp1f1

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 25, 2014

Just tested my branch, and there is some instability remaining

octave-cli:4> hyperg_1F1(50,100,0.01)
ans = 1.0050
wolfram alpha also agrees on this one

In [7]: from scipy.special import hyp1f1

In [8]: hyp1f1(50,100,0.01)
out[8]: 1774259915037.1118


octave-cli:3> hyperg_1F1(0.01,150,-4)
ans = 0.99974

In [5]: hyp1f1(0.01,150,-4)
Out[5]: -80933155.990513608

So, even with that fix, some issues sadly remains. On a positive note, this one works though

octave-cli:4> hyperg_1F1(0.5,1.5,-1000)
ans = 0.028025

In [2]: hyp1f1(0.5,1.5,-1000)
Out[2]: 0.028024956081989651

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 25, 2014

Yes, looks like that code is hard to work with. How about using the mpmath implementation, but using floats / complex?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 25, 2014

I tried mpmath before, but it's unusably slow. Something running in 2 minutes took approximately 6 hours, since it's symbolic mathematics it's much more slower.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 25, 2014

Sure, I mean porting the mpmath code to use floats / complex.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 25, 2014

I don't see what you mean by that actually, mpmath doesn't accept numpy arrays and runs symbolic calculations. It also already works with floats (and probably complex), but they use their own representation for numbers and whatnot, which makes it unlikely we can pull code from the lib.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 26, 2014

I mean to take the code for _hyp1f1 in mpmath, and change it to use actual floats, integers, complex instead of the mpath versions of these. When this is working we can then port this to Cython, optimize, and use it to iterate over arrays in the same way that scipy uses the current Fortran function.

It would be a fair amount of work, but tractable it seems to me - maybe a week or so. We could then contribute that back to scipy.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 26, 2014

I think it's gonna be even worse, since the specialized version like hyp1f1, hyp2f1 and so on are simple calls to the general hypergeometric function, which they are subcases of.

All of the hyper* functions call ctx.hypersum, which seems to be the general formula for the infinite sum fo hypergeometric function. Since they evaluate it at the very end I guess by keeping the symbolic formulation, it amounts to one (slow) evaluation at the end, while evaluating each iteration numerically is (probably) gonna blow up.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 26, 2014

I'm looking here: https://github.com/fredrik-johansson/mpmath/blob/master/mpmath/functions/hypergeometric.py#L311

This uses the hypercomb function at https://github.com/fredrik-johansson/mpmath/blob/master/mpmath/functions/hypergeometric.py#L59

It also uses hypsum, as you say, which, for the float / complex context is:

In [9]: mpmath.fp.hypsum??
Type:        instancemethod
String form: <bound method FPContext.hypsum of <mpmath.ctx_fp.FPContext object at 0x10326b950>>
File:        /Users/mb312/dev_trees/mpmath/mpmath/ctx_fp.py
Definition:  mpmath.fp.hypsum(ctx, p, q, types, coeffs, z, maxterms=6000, **kwargs)
Source:
    def hypsum(ctx, p, q, types, coeffs, z, maxterms=6000, **kwargs):
        coeffs = list(coeffs)
        num = range(p)
        den = range(p,p+q)
        tol = ctx.eps
        s = t = 1.0
        k = 0
        while 1:
            for i in num: t *= (coeffs[i]+k)
            for i in den: t /= (coeffs[i]+k)
            k += 1; t /= k; t *= z; s += t
            if abs(t) < tol:
                return s
            if k > maxterms:
                raise ctx.NoConvergence
@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 26, 2014

This seems to be evaluating the function in a direct way, much like matlab code. The difference as to why this one works lies in using these special types, since they do not store limited floating point precision number.

They are able to accumulate them and give an evaluation at the end completely, while the numerical version will have a division resulting in a nan/inf or other unstable term midway in the summation, giving a false result. Chaging that to numerical code will results in instability, which mpmath tackles by not doing numerical evaluation, but at the cost of speed.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 26, 2014

I can imagine that might be the case, that the gain in precision is from using some arbitrary precision intermediate step, but I can't see where that would happen - can you point to the part of the code that is not using the input float precision, for the floating point (mpmath.fp) context - as in :

In [11]: mpmath.fp.hyp1f1(0.01,150,-4)
Out[11]: 0.9997368389767753
@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 26, 2014

Not sure I understand all that, but this inner function hypsum you posted only runs on floating points number and not their fancy symbolic representation? If that's the case, I'm surprised it works and that every other library seems to struggle with the function.

Joking aside, they seem to be precomputing a huge list of coefficients for the infinite sum and sum over them until they add up to almost nothing. While that is most certainly not efficient due to the precomputing, if I got that analysis right and you think it would be doable in a reasonable computational time, it might be worth a shot.

They go up to precomputing 6000 coefficients, maybe we could go by increments of 100 and count less stuff by declaring convergence before reaching 6000 terms, saving a huge ass time in the process.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 28, 2014

So I went ahead and redid a small version of the sum approach limited to the 1f1 case, and surpringsingly enough it works sometimes, but suffer from slow convergence and doesn't always work (increasing iterations or reducing tolerance doesn't work either).

scipy
In [31]: sci1f1(-4.5,1,-75)
Out[31]: 6767147.0416893596

mine
In [20]: hyp1f1(-4.5,1,-75)
Out[20]: 3591155.3292986257

mpmath
In [21]: mpmath.hyp1f1(-4.5,1,-75)
Out[21]: mpf('6767147.0416894006')

But my implementation works on cases where scipy does not
scipy
In [8]: hyp1f1(50,100,0.01)
out[8]: 1774259915037.1118

mine
In [26]: hyp1f1(50,100,0.01)
Out[26]: 1.0050126452421175

gsl
In [27]: h1f1(50,100,0.01)
Out[27]: array(1.0050126452420993)

Go figure... Seems like we can't have everything at the same time. Here is the code if anyone is interested to have a look.


def hyp1f1(a, b, N, maxterms=6000, eps=1e-8):

    s = 1.
    t = 1.
    k = 0

    while True:

        t *= (a + k)
        t /= (b + k)
        k += 1
        t /= k
        t *= N
        s += t

        if abs(t) < eps:
            return s

        if k > maxterms:
            raise ValueError()
@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 28, 2014

So, if you have a good port of the mpmath algorithm, you should either be getting the same numbers, or be able to explain why you aren't getting the same numbers (mpmath using some unacceptably slow implementation etc).

I think your hyp1f1 is just the last part (hypsum) of the mpmath implementation, and doesn't have the first part _hyp1f1, so maybe that is why you aren't getting the same accuracy. I suspect you will have to port the whole algorithm for that.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 28, 2014

I knew that, hence why I set the general formulation to explicitely work as 1F1 instead of pFq like it was intended originally. I actually have no idea what they are doing to make it work. The whole code checks if the magnitude of Z (number of coils N in our case) is greater than 7, and if not just evaluates the function directly. for N=1 I get mag(N) = 1, N=48 mag=6 and only at N=64 mag=7, which does not explain the not working precision.

Still vote on including the gsl version.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 28, 2014

If you are interested, you really should either port the whole algorithm, or work out exactly why you get different answers to mpmath. Getting an accurate hyp1f1 would be a benefit to the whole scipy ecosystem, not just us, so I think it's worth putting in the effort to do that.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 28, 2014

Well, we are currently also sitting on a re-implementation of NODDI since it also uses hyp1f1, and a bunch of other stuff also, but I don't see the dark magic used by mpmath. There has ot be some more to it that numerical versions can't do, or else I can't see why scipy wouldn't already have used it (it's 10 lines for the general version after all) instead of using a fortran version using another algorithm (there are mathematical papers specifically adressing numerical ways to compute this kind of series)

The fact that the other C implementation of ~2000 lines has 1900 lines dedicated to corner cases is mostly what makes me think it's not possible to do it easily numerically.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 29, 2014

So you are pretty sure that if I do the port myself I won't get the same precision as mpmath?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 29, 2014

Feel free to try, but unless I missed something in the call from _hyp1f1 to hypsum, the direct evaluation is bound to fail in cases where the function converges slowly. Scipy uses an algorithm based on a bunch of log-gamma evaluation, which yields nan in some cases, but there are many ways to evaluate the function near problematic cases, which is the main problem as each cases relies on a different method/algorithm depending on the parameters supplied.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 29, 2014

OK - I played with this for a while. It would take some work to port completely, but I am now pretty sure that mpmath isn't doing any arbitrary precision stuff when in the floating point context:

https://gist.github.com/matthew-brett/0596950032608d32f247

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 29, 2014

Now not using mpmath at all, precision seems to be as for mpmath function:

In [1]: run hyp1f1_port.py

In [2]: hyp1f1(-4.5,1,-75), hyp1f1(50,100,0.01), hyp1f1(0.01,150,-4)
Out[2]: (array(6767147.041689399), 1.0050126452421462, 0.9997368389767753)
@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 29, 2014

Well, my hat to you sir, as this mean they could remove all that not that stable fortran code for all their hypergeometric functions. We could still only keep an internal version of hyp1f1 fixed, as it would be faster and we can discard those check for complex data, I'll try to write something from that effect using that gist.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Jan 5, 2015

Found something else : turns out hte cephes library also implements the hyp1f1 function, and is included in scipy. For some reason, they simply chose to not use it, but rely on something else for this particular one, but use cephes for a bunch of other hypergeometric series function.

https://github.com/scipy/scipy/blob/master/scipy/special/generate_ufuncs.py at line 101

Any idea why/if we could easily build the support for a new ufunc? With the hope that it actually works, since they don't use it for some reason.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Jan 7, 2015

Sam - why not investigate the cephes hyp1f1 function to see if it would be a better choice?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Jan 7, 2015

I posted on the scipy mailing list, if they don't use it, it must have a major flaw or something... Let's hope not and we'll see, it's the last hurdle actually to having adaptive noise estimation for all noise cases, and I don't think any other method has that so far.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Jan 7, 2015

I would not assume anything, always check. If I was doing this, I would work out how to use the cephes implementation with the current scipy mpmat test suite (I think that is skipped at the moment because of the failures).

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Jan 7, 2015

I'm trying to wrap it in cython actually meanwhile, with weird errors at
runtime alas.
Le 2015-01-07 08:33, Matthew Brett a écrit :

Sam - why not investigate the cephes hyp1f1 function to see if it
would be a better choice?


Reply to this email directly or view it on GitHub
#463 (comment).

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Jan 28, 2015

Well, as previously mentionned I tried wrapping the files with cython, only to get a bunch of errors on import. Any suggestion as to where I could find someone to help debug that? Maybe the lib os too old to compile without patches on modern system, the last update dates back from 2000.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Jan 31, 2015

So, after all the cephes version is also somewhat not working in all cases.

### our current gsl version
In [17]: scilpy1f1(-0.5,12.,-800**2/24**2)
Out[17]: array(9.773432897666439)

### scipy
In [19]: h1f1(-0.5,12.,-800**2/24**2)
Out[19]: nan

### cephes 
In [20]: hyp1f1(-0.5,12.,-800**2/24**2)
Out[20]: -2.125043669054867e+212

### wolfram 
hypergeometric1f1(-0.5,12.,-800**2/24**2)
9.773432897666456957821593082159404575552490725623272008191994

Maybe we should move that somewhere else, it's messy to find it back everythiing I try something. Is it possible to fork a topic in two?

And also, so we stay up to date, this actually serves the purpose of implementing both NODDI and adaptive noise estimation, stabilisation algorithms + some other stuff relevant to mic.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Jan 31, 2015

Oh, there is #468 which does that and I had forgotten about, might be easier to take the discussion there instead of this closed PR.

@argriffing

This comment has been minimized.

argriffing commented Aug 6, 2015

I've merged that scipy PR into the master branch, so now it looks like this:

>>> import numpy as np
>>> import scipy
>>> scipy.__version__
'0.17.0.dev0+89066e8'
>>> np.__version__
'1.10.0.dev0+9232200'
>>> from scipy.special import hyp1f1
>>> hyp1f1(-0.5,12.,-800**2/24**2)
9.7773014400072249

It only agrees with gsl and wolfram to a few digits, but it's better than nan or -inf.

@argriffing

This comment has been minimized.

argriffing commented Aug 6, 2015

That uses integer division on Python 2; it looks like the accuracy is good with true division.

>>> hyp1f1(-0.5,12.,-800.0**2/24.0**2)
9.7734328976664777
@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Aug 6, 2015

Yes I saw it on the original PR, thanks! It's surprising it was left unattended for so long, but maybe that will fix our problem back in #675 with any luck. Not having a gsl dependency is easier after all if the scipy version works reliably with this fix.

@argriffing

This comment has been minimized.

argriffing commented Aug 6, 2015

It's surprising it was left unattended for so long

It's not so surprising to me, for anyone in academia or industry or government there are pretty much only disincentives to spend time on it :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment