Re: r25411 - in /trunk: specific_analyses/relax_disp/ test_suite/system_tests/ user_functions/ -- August 29, 2014 - 11:57

Hi,

I just looked at the system test and to me it looks like the '-2'
setting is in the section with:

self.interpreter.relax_disp.r2eff_err_estimate(chi2_jacobian=False)

And the '-1' setting is in the section with:

self.interpreter.relax_disp.r2eff_err_estimate(chi2_jacobian=True)

Is this correct?

Cheers,

Edward




On 29 August 2014 11:30, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:
> Hi Edward.
>
> The new Jacobian, does not do the Job:
>
> -2, is the implemented C code chi2 Jacobian.
> -1, is:
>
> -------
>     # Make partial derivative, with respect to r2eff.
>     d_chi2_d_r2eff = 2.0 * i0 * times * ( -i0 * exp( -r2eff * times) +
> values) * exp( -r2eff * times ) / errors**2
>
>     # Make partial derivative, with respect to i0.
>     d_chi2_d_i0 = - 2.0 * ( -i0 * exp( -r2eff * times) + values) *
> exp( -r2eff * times) / errors**2
> ---------
>
> Results from:
> Relax_disp.verify_estimate_r2eff_err_compare_mc
>
>       -2 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
>       -1 0.052 0.023 0.034 0.021 0.020 0.041 0.030 0.028 0.011 0.163
> 0.048 0.012 0.009 0.010 sum= 0.502
>        0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
> 0.000 0.000 0.000 0.000 sum= 0.000
>       10 0.041 0.052 0.040 0.042 0.040 0.037 0.027 0.017 0.027 0.043
> 0.008 0.015 0.010 0.006 sum= 0.405
>       20 0.029 0.036 0.037 0.040 0.040 0.028 0.050 0.035 0.035 0.042
> 0.017 0.018 0.010 0.007 sum= 0.424
>       30 0.036 0.038 0.036 0.054 0.042 0.036 0.049 0.034 0.034 0.038
> 0.014 0.018 0.010 0.008 sum= 0.447
>       40 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
>
> 2014-08-29 11:01 GMT+02:00 Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx>:
> > 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