mailRe: r25379 - in /trunk: specific_analyses/relax_disp/estimate_r2eff.py test_suite/system_tests/relax_disp.py


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

Header


Content

Posted by Edward d'Auvergne on August 28, 2014 - 16:26:
Hi,

Just a hint, the answer should look a lot like the chi-squared
gradient (http://www.nmr-relax.com/manual/chi_squared_gradient.html),
but probably without the sum and a few other differences.

Regards,

Edward



On 28 August 2014 15:41, Edward d'Auvergne <edward@xxxxxxxxxxxxx> wrote:

Hi,

No, that's not correct.  Try performing the maths yourself and try to
derive the chi-squared partial derivative.  You will see that it's a
little different.

Regards,

Edward



On 28 August 2014 15:38, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> 
wrote:

Hi Edward.

I is just target_functions/c_chi2.c where you dont sum the elements, but
return as array.

Best
Troels

2014-08-28 15:30 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:

Hi Troels,

Could you derive the chi-squared Jacobian?  Maybe the Jacobian I have
been using is not correct - this is the one required for the
Levenberg-Marquardt optimisation algorithm.  Because the chi-squared
is squared, its derivative will have a factor of 2 out the front, just
like the gradient:

http://www.nmr-relax.com/manual/chi_squared_gradient.html

It might be useful to add a Jacobian section to this part of the
manual with the equations.

Cheers,

Edward



On 28 August 2014 15:14,  <tlinnet@xxxxxxxxxxxxx> wrote:

Author: tlinnet
Date: Thu Aug 28 15:14:16 2014
New Revision: 25379

URL: http://svn.gna.org/viewcvs/relax?rev=25379&view=rev
Log:
Modified systemtest test Relax_disp.test_estimate_r2eff_err_methods() to 
show the difference between using the direct function Jacobian, or the 
chi2 function Jacobian.

Added also the functionality to the estimate R2eff module, to switch 
between using the different Jacobians.

The results show, that R2eff can be estimated better.

----------------------
The results are:

R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 431.0.
r2eff=8.646/8.646 r2eff_err=0.0348/0.0692 i0=202664.191/202664.191 
i0_err=699.6443/712.4201

R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 651.2.
r2eff=10.377/10.377 r2eff_err=0.0403/0.0810 i0=206049.558/206049.558 
i0_err=776.4215/782.1833

R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 800.5.
r2eff=10.506/10.506 r2eff_err=0.0440/0.0853 i0=202586.332/202586.332 
i0_err=763.9678/758.7052

R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 984.0.
r2eff=10.903/10.903 r2eff_err=0.0476/0.0922 i0=203455.021/203455.021 
i0_err=837.8779/828.7280

R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1341.1.
r2eff=10.684/10.684 r2eff_err=0.0446/0.0853 i0=218670.412/218670.412 
i0_err=850.0210/830.9558

R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1648.5.
r2eff=10.501/10.501 r2eff_err=0.0371/0.0742 i0=206502.512/206502.512 
i0_err=794.0523/772.9843

R1rho at 799.8 MHz, for offset=124.247 ppm and dispersion point 1341.1.
r2eff=11.118/11.118 r2eff_err=0.0413/0.0827 i0=216447.241/216447.241 
i0_err=784.6562/788.0384

R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 800.5.
r2eff=7.866/7.866 r2eff_err=0.0347/0.0695 i0=211869.715/211869.715 
i0_err=749.2776/763.6930

R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1341.1.
r2eff=9.259/9.259 r2eff_err=0.0331/0.0661 i0=217703.151/217703.151 
i0_err=682.2137/685.5838

R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1648.5.
r2eff=9.565/9.565 r2eff_err=0.0373/0.0745 i0=211988.939/211988.939 
i0_err=839.0313/827.0373

R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 800.5.
r2eff=3.240/3.240 r2eff_err=0.0127/0.0253 i0=214417.382/214417.382 
i0_err=595.8865/613.4378

R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 1341.1.
r2eff=5.084/5.084 r2eff_err=0.0177/0.0352 i0=226358.691/226358.691 
i0_err=660.5314/655.7670

R1rho at 799.8 MHz, for offset=179.768 ppm and dispersion point 1341.1.
r2eff=2.208/2.208 r2eff_err=0.0091/0.0178 i0=228620.553/228620.553 
i0_err=564.8353/560.0873

R1rho at 799.8 MHz, for offset=241.459 ppm and dispersion point 1341.1.
r2eff=1.711/1.711 r2eff_err=0.0077/0.0155 i0=224087.486/224087.486 
i0_err=539.4300/546.4217

Fitting with minfx to: 52V @N
-----------------------------

min_algor='Newton', c_code=True, constraints=False, chi2_jacobian?=False
------------------------------------------------------------------------

R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 431.0, 
with 4 time points. r2eff=8.646 r2eff_err=0.0692, i0=202664.2, 
i0_err=712.4201, chi2=3.758.
R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 651.2, 
with 5 time points. r2eff=10.377 r2eff_err=0.0810, i0=206049.6, 
i0_err=782.1833, chi2=27.291.
R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 800.5, 
with 5 time points. r2eff=10.506 r2eff_err=0.0853, i0=202586.3, 
i0_err=758.7052, chi2=13.357.
R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 984.0, 
with 5 time points. r2eff=10.903 r2eff_err=0.0922, i0=203455.0, 
i0_err=828.7280, chi2=33.632.
R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1341.1, 
with 5 time points. r2eff=10.684 r2eff_err=0.0853, i0=218670.4, 
i0_err=830.9558, chi2=35.818.
R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1648.5, 
with 5 time points. r2eff=10.501 r2eff_err=0.0742, i0=206502.5, 
i0_err=772.9843, chi2=7.356.
R1rho at 799.8 MHz, for offset=124.247 ppm and dispersion point 1341.1, 
with 5 time points. r2eff=11.118 r2eff_err=0.0827, i0=216447.2, 
i0_err=788.0384, chi2=15.587.
R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 800.5, 
with 5 time points. r2eff=7.866 r2eff_err=0.0695, i0=211869.7, 
i0_err=763.6930, chi2=14.585.
R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1341.1, 
with 5 time points. r2eff=9.259 r2eff_err=0.0661, i0=217703.2, 
i0_err=685.5838, chi2=79.498.
R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1648.5, 
with 5 time points. r2eff=9.565 r2eff_err=0.0745, i0=211988.9, 
i0_err=827.0373, chi2=0.447.
R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 800.5, 
with 5 time points. r2eff=3.240 r2eff_err=0.0253, i0=214417.4, 
i0_err=613.4378, chi2=1.681.
R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 1341.1, 
with 5 time points. r2eff=5.084 r2eff_err=0.0352, i0=226358.7, 
i0_err=655.7670, chi2=23.170.
R1rho at 799.8 MHz, for offset=179.768 ppm and dispersion point 1341.1, 
with 5 time points. r2eff=2.208 r2eff_err=0.0178, i0=228620.6, 
i0_err=560.0873, chi2=7.794.
R1rho at 799.8 MHz, for offset=241.459 ppm and dispersion point 1341.1, 
with 5 time points. r2eff=1.711 r2eff_err=0.0155, i0=224087.5, 
i0_err=546.4217, chi2=21.230.

Fitting with minfx to: 52V @N
-----------------------------

min_algor='BFGS', c_code=False, constraints=False, chi2_jacobian?=True
----------------------------------------------------------------------

R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 431.0, 
with 4 time points. r2eff=8.646 r2eff_err=0.0524, i0=202664.2, 
i0_err=1239.0827, chi2=3.758.
R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 651.2, 
with 5 time points. r2eff=10.377 r2eff_err=0.0228, i0=206049.6, 
i0_err=178.1907, chi2=27.291.
R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 800.5, 
with 5 time points. r2eff=10.506 r2eff_err=0.0345, i0=202586.3, 
i0_err=705.7630, chi2=13.357.
R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 984.0, 
with 5 time points. r2eff=10.903 r2eff_err=0.0206, i0=203455.0, 
i0_err=186.0857, chi2=33.632.
R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1341.1, 
with 5 time points. r2eff=10.684 r2eff_err=0.0198, i0=218670.4, 
i0_err=165.0420, chi2=35.818.
R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1648.5, 
with 5 time points. r2eff=10.501 r2eff_err=0.0407, i0=206502.5, 
i0_err=321.3685, chi2=7.356.
R1rho at 799.8 MHz, for offset=124.247 ppm and dispersion point 1341.1, 
with 5 time points. r2eff=11.118 r2eff_err=0.0301, i0=216447.2, 
i0_err=248.9394, chi2=15.587.
R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 800.5, 
with 5 time points. r2eff=7.866 r2eff_err=0.0280, i0=211869.7, 
i0_err=259.8845, chi2=14.585.
R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1341.1, 
with 5 time points. r2eff=9.259 r2eff_err=0.0108, i0=217703.2, 
i0_err=88.1514, chi2=79.498.
R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1648.5, 
with 5 time points. r2eff=9.565 r2eff_err=0.1630, i0=211988.9, 
i0_err=2054.6615, chi2=0.447.
R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 800.5, 
with 5 time points. r2eff=3.240 r2eff_err=0.0485, i0=214417.4, 
i0_err=611.7573, chi2=1.681.
R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 1341.1, 
with 5 time points. r2eff=5.084 r2eff_err=0.0124, i0=226358.7, 
i0_err=122.7341, chi2=23.170.
R1rho at 799.8 MHz, for offset=179.768 ppm and dispersion point 1341.1, 
with 5 time points. r2eff=2.208 r2eff_err=0.0086, i0=228620.6, 
i0_err=219.4208, chi2=7.794.
R1rho at 799.8 MHz, for offset=241.459 ppm and dispersion point 1341.1, 
with 5 time points. r2eff=1.711 r2eff_err=0.0101, i0=224087.5, 
i0_err=166.9081, chi2=21.230.

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

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=25379&r1=25378&r2=25379&view=diff
==============================================================================
--- trunk/specific_analyses/relax_disp/estimate_r2eff.py        
(original)
+++ trunk/specific_analyses/relax_disp/estimate_r2eff.py        Thu Aug 
28 15:14:16 2014
@@ -175,7 +175,7 @@
                     print(print_string),


-def multifit_covar(J=None, epsrel=0.0, errors=None):
+def multifit_covar(J=None, epsrel=0.0, errors=None, use_weights=True):
     """This is the implementation of the multifit covariance.

     This is inspired from GNU Scientific Library (GSL).
@@ -184,9 +184,15 @@

     The parameter 'epsrel' is used to remove linear-dependent columns 
when J is rank deficient.

+    The weighting matrix 'W', is a square symmetric matrix. For 
independent measurements, this is a diagonal matrix. Larger values 
indicate greater significance.  It is formed by multiplying the supplied 
errors as 1./errors^2 with an Identity matrix::
+
+        W = I.(1/errors^2)
+
+    If 'use_weights' is set to 'False', the errors are set to 1.0.
+
     The covariance matrix is given by::

-        covar = (J^T J)^{-1} ,
+        covar = (J^T.W.J)^{-1} ,

     and is computed by QR decomposition of J with column-pivoting. Any 
columns of R which satisfy::

@@ -224,6 +230,8 @@
     @type epsrel:           float
     @keyword errors:        The standard deviation of the measured 
intensity values per time point.
     @type errors:           numpy array
+    @keyword use_weights:   If the supplied weights should be used.
+    @type use_weights:      bool
     @return:                The co-variance matrix
     @rtype:                 square numpy array
     """
@@ -237,6 +245,10 @@
     # Now form the error matrix, with errors down the diagonal.
     weights = 1. / errors**2

+    if use_weights == False:
+        weights[:] = 1.0
+
+    # Form weight matrix.
     W = multiply(weights, eye_mat)

     # The covariance matrix (sometimes referred to as the 
variance-covariance matrix), Qxx, is defined as:
@@ -344,7 +356,7 @@
         self.factor = factor


-    def set_settings_minfx(self, scaling_matrix=None, 
min_algor='simplex', c_code=True, constraints=False, func_tol=1e-25, 
grad_tol=None, max_iterations=10000000):
+    def set_settings_minfx(self, scaling_matrix=None, 
min_algor='simplex', c_code=True, constraints=False, 
chi2_jacobian=False, func_tol=1e-25, grad_tol=None, 
max_iterations=10000000):
         """Setup options to minfx.

         @keyword scaling_matrix:    The square and diagonal scaling 
matrix.
@@ -355,6 +367,8 @@
         @type c_code:               bool
         @keyword constraints:       If constraints should be used.
         @type constraints:          bool
+        @keyword chi2_jacobian:     If the chi2 Jacobian should be used.
+        @type chi2_jacobian:        bool
         @keyword func_tol:          The function tolerance which, when 
reached, terminates optimisation.  Setting this to None turns of the 
check.
         @type func_tol:             None or float
         @keyword grad_tol:          The gradient tolerance which, when 
reached, terminates optimisation.  Setting this to None turns of the 
check.
@@ -366,6 +380,7 @@
         # Store variables.
         self.scaling_matrix = scaling_matrix
         self.c_code = c_code
+        self.chi2_jacobian = chi2_jacobian

         # Scaling initialisation.
         self.scaling_flag = False
@@ -561,7 +576,7 @@
         return 1. / self.errors * (self.func_exp(self.times, *params) - 
self.values)


-def estimate_r2eff(method='minfx', min_algor='simplex', c_code=True, 
constraints=False, spin_id=None, ftol=1e-15, xtol=1e-15, 
maxfev=10000000, factor=100.0, verbosity=1):
+def estimate_r2eff(method='minfx', min_algor='simplex', c_code=True, 
constraints=False, chi2_jacobian=False, spin_id=None, ftol=1e-15, 
xtol=1e-15, maxfev=10000000, factor=100.0, verbosity=1):
     """Estimate r2eff and errors by exponential curve fitting with 
scipy.optimize.leastsq or minfx.

     THIS IS ONLY FOR TESTING.
@@ -583,10 +598,12 @@
     @type method:               string
     @keyword min_algor:         The minimisation algorithm
     @type min_algor:            string
+    @keyword c_code:            If optimise with C code.
+    @type c_code:               bool
     @keyword constraints:       If constraints should be used.
     @type constraints:          bool
-    @keyword c_code:            If optimise with C code.
-    @type c_code:               bool
+    @keyword chi2_jacobian:     If the chi2 Jacobian should be used.
+    @type chi2_jacobian:        bool
     @keyword spin_id:           The spin identification string.
     @type spin_id:              str
     @keyword ftol:              The function tolerance for the relative 
error desired in the sum of squares, parsed to leastsq.
@@ -661,7 +678,7 @@
                 top += 2
             subsection(file=sys.stdout, text="Fitting with %s to: 
%s"%(method, spin_string), prespace=top)
             if method == 'minfx':
-                subsection(file=sys.stdout, text="min_algor='%s', 
c_code=%s, constraints=%s"%(min_algor, c_code, constraints), prespace=0)
+                subsection(file=sys.stdout, text="min_algor='%s', 
c_code=%s, constraints=%s, chi2_jacobian?=%s"%(min_algor, c_code, 
constraints, chi2_jacobian), prespace=0)

         # Loop over each spectrometer frequency and dispersion point.
         for exp_type, frq, offset, point, ei, mi, oi, di in 
loop_exp_frq_offset_point(return_indices=True):
@@ -692,7 +709,7 @@

             elif method == 'minfx':
                 # Set settings.
-                E.set_settings_minfx(min_algor=min_algor, 
c_code=c_code, constraints=constraints)
+                E.set_settings_minfx(min_algor=min_algor, 
c_code=c_code, chi2_jacobian=chi2_jacobian, constraints=constraints)

                 # Acquire results.
                 results = minimise_minfx(E=E)
@@ -737,7 +754,7 @@
                 point_info = "%s at %3.1f MHz, for offset=%3.3f ppm and 
dispersion point %-5.1f, with %i time points." % (exp_type, frq/1E6, 
offset, point, len(times))
                 print_strings.append(point_info)

-                par_info = "r2eff=%3.3f r2eff_err=%3.3f, i0=%6.1f, 
i0_err=%3.3f, chi2=%3.3f.\n" % ( r2eff, r2eff_err, i0, i0_err, chi2)
+                par_info = "r2eff=%3.3f r2eff_err=%3.4f, i0=%6.1f, 
i0_err=%3.4f, chi2=%3.3f.\n" % ( r2eff, r2eff_err, i0, i0_err, chi2)
                 print_strings.append(par_info)

                 if E.verbosity >= 2:
@@ -912,14 +929,24 @@
         #jacobian_matrix_exp2 = E.jacobian_matrix_exp
         #print jacobian_matrix_exp - jacobian_matrix_exp2
     else:
-        # Call class, to store value.
-        E.func_exp_grad(param_vector)
-        jacobian_matrix_exp = E.jacobian_matrix_exp
-        #E.func_exp_chi2_grad(param_vector)
-        #jacobian_matrix_exp = E.jacobian_matrix_exp_chi2
+        if E.chi2_jacobian:
+            # Call class, to store value.
+            E.func_exp_chi2_grad(param_vector)
+            jacobian_matrix_exp = E.jacobian_matrix_exp_chi2
+        else:
+            # Call class, to store value.
+            E.func_exp_grad(param_vector)
+            jacobian_matrix_exp = E.jacobian_matrix_exp
+            #E.func_exp_chi2_grad(param_vector)
+            #jacobian_matrix_exp = E.jacobian_matrix_exp_chi2

     # Get the co-variance
-    pcov = multifit_covar(J=jacobian_matrix_exp, errors=E.errors)
+    if E.chi2_jacobian:
+        use_weights = False
+    else:
+        use_weights = True
+
+    pcov = multifit_covar(J=jacobian_matrix_exp, errors=E.errors, 
use_weights=use_weights)

     # To compute one standard deviation errors on the parameters, take 
the square root of the diagonal covariance.
     param_vector_error = sqrt(diag(pcov))

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=25379&r1=25378&r2=25379&view=diff
==============================================================================
--- trunk/test_suite/system_tests/relax_disp.py (original)
+++ trunk/test_suite/system_tests/relax_disp.py Thu Aug 28 15:14:16 2014
@@ -2946,12 +2946,13 @@


         # Now do it manually.
-        estimate_r2eff(method='scipy.optimize.leastsq')
-        estimate_r2eff(method='minfx', min_algor='simplex', 
c_code=True, constraints=False)
-        estimate_r2eff(method='minfx', min_algor='simplex', 
c_code=False, constraints=False)
-        estimate_r2eff(method='minfx', min_algor='BFGS', c_code=True, 
constraints=False)
-        estimate_r2eff(method='minfx', min_algor='BFGS', c_code=False, 
constraints=False)
+        #estimate_r2eff(method='scipy.optimize.leastsq')
+        #estimate_r2eff(method='minfx', min_algor='simplex', 
c_code=True, constraints=False)
+        #estimate_r2eff(method='minfx', min_algor='simplex', 
c_code=False, constraints=False)
+        #estimate_r2eff(method='minfx', min_algor='BFGS', c_code=True, 
constraints=False)
+        #estimate_r2eff(method='minfx', min_algor='BFGS', c_code=False, 
constraints=False)
         estimate_r2eff(method='minfx', min_algor='Newton', c_code=True, 
constraints=False)
+        estimate_r2eff(method='minfx', min_algor='BFGS', c_code=False, 
constraints=False, chi2_jacobian=True)


     def test_exp_fit(self):


_______________________________________________
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 10:00:16 2014