mailRe: r25411 - in /trunk: specific_analyses/relax_disp/ test_suite/system_tests/ user_functions/


Others Months | Index by Date | Thread Index
>>   [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Header


Content

Posted by Edward d'Auvergne on August 29, 2014 - 12:03:
Could you create a simple relax script and log file in the test suite
directories with the results of a super-massive Monte Carlo simulation
run?  I.e. a really insane number of Monte Carlo simulations and then
a clear print out of both the R2eff and I0 errors.  I really cannot
follow the numbers in these tests.  There are too many points where
bugs could be hidden.  It would be good to also have a simple relax
script and log file for the covariance matrix estimate as well.

Cheers,

Edward

On 29 August 2014 11:59, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:
You may want to look here:

relax -s Relax_disp.test_estimate_r2eff_err_methods -d

2014-08-29 11:57 GMT+02:00 Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>:
Hi Edward.

There is something totally wrong with the C, Jacobian.
Errors are estimated to:

37.619 17.290 25.616 16.036 16.164 32.826 22.920 21.462 7.777 145.309
36.884 9.116 6.199 7.018 sum= 402.235

Which is much different to:
0.041 0.040 0.040 0.054 0.041 0.044 0.042 0.037 0.034 0.043 0.013
0.018 0.007 0.010 sum= 0.462

You can see how the error estimation develops in:
verify_estimate_r2eff_err_compare_mc

You will see, that just 50 monte carlo simulations is better than 
estimating.

Best
Troels


2014-08-29 11:51 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
Hi,

I saw the results from that 'hidden' system test and was wondering
what was happening?  The Jacobian of the chi-squared function should
remove the factor of 2, as it has a factor of minus two.  But it also
includes the difference between the measured and back-calculated peak
intensities divided by the variance as well.  So why does this
Jacobian, which is much closer to the 2000 MC simulations, not work?
I cannot understand this as it is totally illogical.  If your error
estimate is closer to the real thing, then you should get closer to
the real optimisation results.

Do you have a log file somewhere which contains the results from the
2000 MC simulations?  It might be worth creating a file which compares
this, or even more simulations, 100,000 for example, to the covariance
technique.  Once the error estimate technique is functional and
debugged, then we can work out why the models are optimisating
differently.  These two problems need to be separated and solved
independently, otherwise you can encounter the common yet fatal coding
problem of two opposing bugs partially cancelling out their effects.

Regards,

Edward

On 29 August 2014 11:01, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> 
wrote:
Hi Edward.

Would it be possible to have both?

The exponential Jacobian, and the chi2 Jacobian.

My tests last night showed something weird.

Using the chi2 Jacobian, the errors come closer to the ones reported
my MC calculations.
The direct jacobian would have double error on R2eff.

But when fitting for R1rho models, using the errors from the direct
jacobian, was much better in agreement with
MC error fitting.

The parameters from chi2 Jacobian, was worse.

See verify_r1rho_kjaergaard_missing_r1() in systemtest for comparison.

Look at the 'kex' parameter!

# Compare values.
if spin_id == ':52@N':
    if param == 'r1':
        if model == MODEL_NOREX:
            if r2eff_estimate == 'direct':
                self.assertAlmostEqual(value, 1.46138805)
            elif r2eff_estimate == 'MC2000':
                self.assertAlmostEqual(value, 1.46328102)
            elif r2eff_estimate == 'chi2':
                self.assertAlmostEqual(value, 1.43820629)
        elif model == MODEL_DPL94:
            if r2eff_estimate == 'direct':
                self.assertAlmostEqual(value, 1.44845742)
            elif r2eff_estimate == 'MC2000':
                self.assertAlmostEqual(value, 1.45019848)
            elif r2eff_estimate == 'chi2':
                self.assertAlmostEqual(value, 1.44666512)
        elif model == MODEL_TP02:
            if r2eff_estimate == 'direct':
                self.assertAlmostEqual(value, 1.54354392)
            elif r2eff_estimate == 'MC2000':
                self.assertAlmostEqual(value, 1.54352369)
            elif r2eff_estimate == 'chi2':
                self.assertAlmostEqual(value, 1.55964020)
        elif model == MODEL_TAP03:
            if r2eff_estimate == 'direct':
                self.assertAlmostEqual(value, 1.54356410)
            elif r2eff_estimate == 'MC2000':
                self.assertAlmostEqual(value, 1.54354367)
            elif r2eff_estimate == 'chi2':
                self.assertAlmostEqual(value, 1.55967157)
        elif model == MODEL_MP05:
            if r2eff_estimate == 'direct':
                self.assertAlmostEqual(value, 1.54356416)
            elif r2eff_estimate == 'MC2000':
                self.assertAlmostEqual(value, 1.54354372)
            elif r2eff_estimate == 'chi2':
                self.assertAlmostEqual(value, 1.55967163)
        elif model == MODEL_NS_R1RHO_2SITE:
            if r2eff_estimate == 'direct':
                self.assertAlmostEqual(value, 1.41359221, 5)
            elif r2eff_estimate == 'MC2000':
                self.assertAlmostEqual(value, 1.41321968, 5)
            elif r2eff_estimate == 'chi2':
                self.assertAlmostEqual(value, 1.36303129, 5)

    elif param == 'r2':
        if model == MODEL_NOREX:
            if r2eff_estimate == 'direct':
                self.assertAlmostEqual(value, 11.48392439)
            elif r2eff_estimate == 'MC2000':
                self.assertAlmostEqual(value, 11.48040934)
            elif r2eff_estimate == 'chi2':
                self.assertAlmostEqual(value, 11.47224488)
        elif model == MODEL_DPL94:
            if r2eff_estimate == 'direct':
                self.assertAlmostEqual(value, 10.15688372, 6)
            elif r2eff_estimate == 'MC2000':
                self.assertAlmostEqual(value, 10.16304887, 6)
            elif r2eff_estimate == 'chi2':
                self.assertAlmostEqual(value, 9.20037797, 6)
        elif model == MODEL_TP02:
            if r2eff_estimate == 'direct':
                self.assertAlmostEqual(value, 9.72654896, 6)
            elif r2eff_estimate == 'MC2000':
                self.assertAlmostEqual(value, 9.72772726, 6)
            elif r2eff_estimate == 'chi2':
                self.assertAlmostEqual(value, 9.53948340, 6)
        elif model == MODEL_TAP03:
            if r2eff_estimate == 'direct':
                self.assertAlmostEqual(value, 9.72641887, 6)
            elif r2eff_estimate == 'MC2000':
                self.assertAlmostEqual(value, 9.72759374, 6)
            elif r2eff_estimate == 'chi2':
                self.assertAlmostEqual(value, 9.53926913, 6)
        elif model == MODEL_MP05:
            if r2eff_estimate == 'direct':
                self.assertAlmostEqual(value, 9.72641723, 6)
            elif r2eff_estimate == 'MC2000':
                self.assertAlmostEqual(value, 9.72759220, 6)
            elif r2eff_estimate == 'chi2':
                self.assertAlmostEqual(value, 9.53926778, 6)
        elif model == MODEL_NS_R1RHO_2SITE:
            if r2eff_estimate == 'direct':
                self.assertAlmostEqual(value, 9.34531535, 5)
            elif r2eff_estimate == 'MC2000':
                self.assertAlmostEqual(value, 9.34602793, 5)
            elif r2eff_estimate == 'chi2':
                self.assertAlmostEqual(value, 9.17631409, 5)

# For all other parameters.
else:
# Get the value.
value = getattr(cur_spin, param)

# Print value.
print("%-10s %-6s %-6s %3.8f" % ("Parameter:", param, "Value:", value))

# Compare values.
if spin_id == ':52@N':
if param == 'phi_ex':
    if model == MODEL_DPL94:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 0.07599563)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 0.07561937)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 0.12946061)

elif param == 'pA':
    if model == MODEL_TP02:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 0.88827040)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 0.88807487)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 0.87746233)
    elif model == MODEL_TAP03:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 0.88828922)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 0.88809318)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 0.87747558)
    elif model == MODEL_MP05:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 0.88828924)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 0.88809321)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 0.87747562)
    elif model == MODEL_NS_R1RHO_2SITE:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 0.94504369, 6)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 0.94496541, 6)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 0.92084707, 6)

elif param == 'dw':
    if model == MODEL_TP02:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 1.08875840, 6)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 1.08765638, 6)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 1.09753230, 6)
    elif model == MODEL_TAP03:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 1.08837238, 6)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 1.08726698, 6)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 1.09708821, 6)
    elif model == MODEL_MP05:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 1.08837241, 6)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 1.08726706, 6)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 1.09708832, 6)
    elif model == MODEL_NS_R1RHO_2SITE:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 1.56001812, 5)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 1.55833321, 5)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 1.36406712, 5)

elif param == 'kex':
    if model == MODEL_DPL94:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 4460.43711569, 2)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 4419.03917195, 2)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 6790.22736344, 2)
    elif model == MODEL_TP02:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 4921.28602757, 3)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 4904.70144883, 3)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 5146.20306591, 3)
    elif model == MODEL_TAP03:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 4926.42963491, 3)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 4909.86877150, 3)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 5152.51105814, 3)
    elif model == MODEL_MP05:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 4926.44236315, 3)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 4909.88110195, 3)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 5152.52097111, 3)
    elif model == MODEL_NS_R1RHO_2SITE:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 5628.66061488, 2)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 5610.20221435, 2)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 5643.34067090, 2)

elif param == 'chi2':
    if model == MODEL_NOREX:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 848.42016907, 5)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 3363.95829122, 5)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 5976.49946726, 5)
    elif model == MODEL_DPL94:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 179.47041241)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 710.24767560)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 612.72616697, 5)
    elif model == MODEL_TP02:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 29.33882530, 6)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 114.47142772, 6)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 250.50838162, 5)
    elif model == MODEL_TAP03:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 29.29050673, 6)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 114.27987534)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 250.04050719, 5)
    elif model == MODEL_MP05:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 29.29054301, 6)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 114.28002272)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 250.04077478, 5)
    elif model == MODEL_NS_R1RHO_2SITE:
        if r2eff_estimate == 'direct':
            self.assertAlmostEqual(value, 34.44010543, 6)
        elif r2eff_estimate == 'MC2000':
            self.assertAlmostEqual(value, 134.14368365)
        elif r2eff_estimate == 'chi2':
            self.assertAlmostEqual(value, 278.55121388, 5)

2014-08-29 9:49 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
Hi Troels,

I've now converted the target_functions.relax_fit.jacobian() function
to be the Jacobian of the chi-squared function rather than the
Jacobian of the exponential function.  This should match your
specific_analyses.relax_disp.estimate_r2eff.func_exp_chi2_grad()
function.  I mixed up the two because the Levenberg-Marquardt
algorithm in minfx requires the Jacobian of the exponential, and it's
been 8 years since I last derived and implemented a Jacobian.

Regards,

Edward



On 28 August 2014 21:43,  <tlinnet@xxxxxxxxxxxxx> wrote:
Author: tlinnet
Date: Thu Aug 28 21:43:13 2014
New Revision: 25411

URL: http://svn.gna.org/viewcvs/relax?rev=25411&view=rev
Log:
Reverted the logic, that the chi2 Jacobian should be used.

Instead, the direct Jacobian exponential is used instead.

When fitting with the estimated errors from the Direct Jacobian, the 
results are MUCH better, and comparable
to 2000 Monte-Carlo simulations.

task #7822(https://gna.org/task/index.php?7822): Implement user 
function to estimate R2eff and associated errors for exponential curve 
fitting.

Modified:
    trunk/specific_analyses/relax_disp/estimate_r2eff.py
    trunk/test_suite/system_tests/relax_disp.py
    trunk/user_functions/relax_disp.py

Modified: trunk/specific_analyses/relax_disp/estimate_r2eff.py
URL: 
http://svn.gna.org/viewcvs/relax/trunk/specific_analyses/relax_disp/estimate_r2eff.py?rev=25411&r1=25410&r2=25411&view=diff
==============================================================================
--- trunk/specific_analyses/relax_disp/estimate_r2eff.py        
(original)
+++ trunk/specific_analyses/relax_disp/estimate_r2eff.py        Thu 
Aug 28 21:43:13 2014
@@ -90,7 +90,7 @@
     return jacobian_matrix_exp_chi2


-def estimate_r2eff_err(chi2_jacobian=True, spin_id=None, epsrel=0.0, 
verbosity=1):
+def estimate_r2eff_err(chi2_jacobian=False, spin_id=None, epsrel=0.0, 
verbosity=1):
     """This will estimate the R2eff and i0 errors from the covariance 
matrix Qxx.  Qxx is calculated from the Jacobian matrix and the 
optimised parameters.

     @keyword chi2_jacobian: If the Jacobian derived from the chi2 
function, should be used instead of the Jacobian from the exponential 
function.

Modified: trunk/test_suite/system_tests/relax_disp.py
URL: 
http://svn.gna.org/viewcvs/relax/trunk/test_suite/system_tests/relax_disp.py?rev=25411&r1=25410&r2=25411&view=diff
==============================================================================
--- trunk/test_suite/system_tests/relax_disp.py (original)
+++ trunk/test_suite/system_tests/relax_disp.py Thu Aug 28 21:43:13 
2014
@@ -2744,13 +2744,13 @@
         self.interpreter.minimise.execute(min_algor='Newton', 
constraints=False, verbosity=1)

         # Estimate R2eff errors.
-        
self.interpreter.relax_disp.r2eff_err_estimate(chi2_jacobian=False)
+        
self.interpreter.relax_disp.r2eff_err_estimate(chi2_jacobian=True)

         # Run the analysis.
         relax_disp.Relax_disp(pipe_name=ds.pipe_name, 
pipe_bundle=ds.pipe_bundle, results_dir=result_dir_name, 
models=MODELS, grid_inc=GRID_INC, mc_sim_num=MC_NUM, modsel=MODSEL)

         # Verify the data.
-        self.verify_r1rho_kjaergaard_missing_r1(models=MODELS, 
result_dir_name=result_dir_name, r2eff_estimate='direct')
+        self.verify_r1rho_kjaergaard_missing_r1(models=MODELS, 
result_dir_name=result_dir_name, r2eff_estimate='chi2')


     def test_estimate_r2eff_err_auto(self):
@@ -2849,7 +2849,7 @@
         relax_disp.Relax_disp(pipe_name=pipe_name, 
pipe_bundle=pipe_bundle, results_dir=result_dir_name, models=MODELS, 
grid_inc=GRID_INC, mc_sim_num=MC_NUM, exp_mc_sim_num=EXP_MC_NUM, 
modsel=MODSEL, r1_fit=r1_fit)

         # Verify the data.
-        self.verify_r1rho_kjaergaard_missing_r1(models=MODELS, 
result_dir_name=result_dir_name, r2eff_estimate='chi2')
+        self.verify_r1rho_kjaergaard_missing_r1(models=MODELS, 
result_dir_name=result_dir_name, r2eff_estimate='direct')


     def test_estimate_r2eff_err_methods(self):

Modified: trunk/user_functions/relax_disp.py
URL: 
http://svn.gna.org/viewcvs/relax/trunk/user_functions/relax_disp.py?rev=25411&r1=25410&r2=25411&view=diff
==============================================================================
--- trunk/user_functions/relax_disp.py  (original)
+++ trunk/user_functions/relax_disp.py  Thu Aug 28 21:43:13 2014
@@ -636,7 +636,7 @@
 uf.title_short = "Estimate R2eff errors."
 uf.add_keyarg(
     name = "chi2_jacobian",
-    default = True,
+    default = False,
     py_type = "bool",
     desc_short = "use of chi2 Jacobian",
     desc = "If the Jacobian derived from the chi2 function, should be 
used instead of the Jacobian from the exponential function."


_______________________________________________
relax (http://www.nmr-relax.com)

This is the relax-commits mailing list
relax-commits@xxxxxxx

To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-commits

_______________________________________________
relax (http://www.nmr-relax.com)

This is the relax-devel mailing list
relax-devel@xxxxxxx

To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-devel



Related Messages


Powered by MHonArc, Updated Fri Aug 29 12:20:17 2014