Hi,
On 6 May 2014 11:05, Andrew Baldwin <andrew.baldwin@xxxxxxxxxxxxx> wrote:
That CarverRichards number now looks what I would personally expect
from the theory!
I'd agree with that. This is noteworthy - we're agreeing on something! Had
to happen eventually :)
:)
I don't like the Nikolai results, but I also don't
like the fact that I have two completely different implementations of
this code (from the sim_all software and the compareTroells.py
script).
The different versions of code is a good point. I copied and pasted what I
think is your function here:
http://www.nmr-relax.com/api/3.1/lib.dispersion.ns_cpmg_2site_expanded-pysrc.html
And have added this as the fun mode (see function nikolaiFun). My Nikolai
version came from 'reiko_ishima_v1.tar' and 'Exact_Rex_mod_v1.m' as i'm
sure
he'll explain shortly. I ran my python version against the matlab original
to confirm that the transcription was good.
Note in the grid search in this version of testdisp, the grid is no longer
going over R2E, as the funmode is R2g=R2e. Note also that I've jigged up
the
Pb in the attached code to cover higher pbs (we now go to 40%) to increase
the apparent errors in the CR equation. This is well outside anywhere
experimentally accessible, but fun to observe.
Troels, it's worth noting this R2g=R2e. As most people will use this
in the field, it would be of great value to have a 'B14' model variant
whereby R20A == R20B ;)
So your function is already faster by a factor of 2 which is good. Please
can you check the function to make sure i've wired your definitions up to
my
function correctly? When I use a narrow grid search, the function seems to
give the exact result so I think it's probably right.
Here's the o/p:
####
Relative times of the various approaches:
Baldwin 0.194165 1.0
Numerical 21.562522 111.052568692
Meiboom 0.070515 0.363170499318
Carver 0.105395 0.542811526279
compareTroells.py:396: RuntimeWarning: invalid value encountered in log
R2eff=(1/Tc)*numpy.log(intensity0/intensity); # we did not factor
out (Ka+Kb)/2 here
Nikolai 0.615217 3.16852676847
compareTroells.py:274: RuntimeWarning: invalid value encountered in log
R2eff=-inv_relax_time * numpy.log(Mx)
NikolaiFun 0.376738 1.94029819998
NikolaiLong 81.447843 419.477470193
NikolaiLong 82.098139 422.82666289
Maximum error over grid (s-1):
Meiboom: 237214256.582
Baldwin: 6.57168541807e-10
CarverRichards: 29.9287668504
Nikolai fun: 2964.41300775
Nikolai dougle (9): 142.46676955
Nikolai long double (18): 0.063449549192819563320223
Nikolai long double (23): 5.7775487944482434059452e-7
####
The tcp in this function in relax is defined as 1/(4*nu_CPMG), but in
your script you have "tcp=1./(2*nu_cpmg);". Apart from that, it looks
ok. Changing to the correct tcp gives:
Relative times of the various approaches:
Baldwin 0.159632 1.0
Numerical 19.974315 125.127261451
Meiboom 0.056497 0.353920266613
Carver 0.080027 0.501321790117
Nikolai 0.453068 2.83820286659
NikolaiFun 0.284052 1.77941766062
NikolaiLong 50.020116 313.34642177
NikolaiLong 50.430594 315.917823494
Maximum error over grid (s-1):
Meiboom: 237214256.582
Baldwin: 6.57168541807e-10
CarverRichards: 29.9287668504
Nikolai fun: 104.735381561
Nikolai dougle (9): 142.46676955
Nikolai long double (18): 0.063449549192819563320223
Nikolai long double (23): 5.7775487944482434059452e-7
That's far more reasonable for Nikolai's fun!
I was hoping that the modification would make it more stable. If anything
the reverse is true. I would have included a note to this effect in the
paper (and still will if you can show it's better). Note the error is very
temperamental. You can miss it even by setting up the grid search
differently. Give it a narrow grid, and you won't see the error. Only
relatively specific, but seemingly hard to predict sets of parameters (all
with large deltaOs) lead it to explode.
It is worth looking at all deviations rather than just the single
maximum R2eff difference. It could be that there is an edge case in
the grid whereby a maths domain error occurs. Or an trig function
overflow. The median deviation to avoid such outliers might also be
of worth. By switching in your code all numpy.max() calls with
numpy.median(), and making sure the "Nikolai long double (*)" results
are pre-converted to a numpy float64 array, I then see an incredibly
different picture:
Relative times of the various approaches:
Baldwin 0.152725 1.0
Numerical 19.596446 128.3119725
Meiboom 0.05604 0.366934031756
Carver 0.077236 0.505719430349
Nikolai 0.43944 2.87732853167
NikolaiFun 0.273977 1.7939237191
NikolaiLong 50.243433 328.979754461
NikolaiLong 50.784841 332.524740547
Median error over grid (s-1):
Meiboom: 4.0064828013
Baldwin: 3.9630521087e-12
CarverRichards: 0.00661151728704
Nikolai fun: 6.62225829728e-12
Nikolai dougle (9): 6.70041799822e-12
Nikolai long double (18): 2.62099543458e-12
Nikolai long double (23): 2.27758174227e-12
!
So this is much clearer - outliers are the major difference you are
seeing. A study of such outliers and their cause would be quite
interesting. But in most of these cases, simple collapsed equations
for each model should be derivable at the edges which then produce the
expected R2eff over the entire parameter space rather than
computational artefacts in some corners. Anyway, this result makes it
crystal clear that the different models should not be judged by the
single maximum outlier over the entire parameter space.
Regards,
Edward
Re: [sr #3154] Implementation of Baldwin (2014) B14 model - 2-site exact solution model for all time scales