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