Author: tlinnet
Date: Thu Aug 28 19:58:11 2014
New Revision: 25404
URL: http://svn.gna.org/viewcvs/relax?rev=25404&view=rev
Log:
Cleaned up code in R2eff estimate module, by making each function independent
of class.
This is to give a better overview, how the different functions connect
together.
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
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=25404&r1=25403&r2=25404&view=diff
==============================================================================
--- trunk/specific_analyses/relax_disp/estimate_r2eff.py (original)
+++ trunk/specific_analyses/relax_disp/estimate_r2eff.py Thu Aug 28
19:58:11 2014
@@ -60,7 +60,7 @@
def func_exp_chi2_grad(params=None, times=None, values=None, errors=None):
- """Target function for the gradient (Jacobian matrix) of func_exp_chi2()
to minfx, for exponential fit .
+ """Return the gradient (Jacobian matrix) of func_exp_chi2() for
exponential fit .
@param params: The vector of parameter values.
@type params: numpy rank-1 float array
@@ -90,7 +90,7 @@
return jacobian_matrix_exp_chi2
-def estimate_r2eff_err(chi2_jacobian=False, spin_id=None, epsrel=0.0,
verbosity=1):
+def estimate_r2eff_err(chi2_jacobian=True, 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.
@@ -361,9 +361,6 @@
self.times = times
self.values = values
self.errors = errors
-
- # Create the structure for holding the back-calculated R2eff values
(matching the dimensions of the values structure).
- self.back_calc = deepcopy(self.values)
def set_settings_leastsq(self, ftol=1e-15, xtol=1e-15, maxfev=10000000,
factor=100.0):
@@ -532,11 +529,17 @@
return Kw
- def func_exp_grad(self, params):
- """Target function for the gradient (Jacobian matrix) of func_exp to
minfx, for exponential fit .
+ def func_exp_grad(self, params=None, times=None, values=None,
errors=None):
+ """The gradient (Jacobian matrix) of func_exp for Co-variance
calculation.
@param params: The vector of parameter values.
@type params: numpy rank-1 float array
+ @keyword times: The time points.
+ @type times: numpy array
+ @param values: The measured values.
+ @type values: numpy array
+ @param errors: The standard deviation of the measured intensity
values per time point.
+ @type errors: numpy array
@return: The Jacobian matrix with 'm' rows of function
derivatives per 'n' columns of parameters.
@rtype: numpy array
"""
@@ -546,68 +549,67 @@
i0 = params[1]
# Make partial derivative, with respect to r2eff.
- d_exp_d_r2eff = -i0 * self.times * exp(-r2eff * self.times)
+ d_exp_d_r2eff = -i0 * times * exp(-r2eff * times)
# Make partial derivative, with respect to i0.
- d_exp_d_i0 = exp(-r2eff * self.times)
+ d_exp_d_i0 = exp(-r2eff * times)
# Define Jacobian as m rows with function derivatives and n columns
of parameters.
- self.jacobian_matrix_exp = transpose(array( [d_exp_d_r2eff ,
d_exp_d_i0] ) )
-
- # Take the sum, to send to minfx.
- sum_d_exp_d_r2eff = sum( d_exp_d_r2eff )
- sum_d_exp_d_i0 = sum( d_exp_d_i0 )
-
- # Define Jacobian as m rows with function derivatives and n columns
of parameters.
- sum_jacobian_matrix_exp = transpose(array( [sum_d_exp_d_r2eff ,
sum_d_exp_d_i0] ) )
+ jacobian_matrix_exp = transpose(array( [d_exp_d_r2eff , d_exp_d_i0]
) )
# Return Jacobian matrix.
- return sum_jacobian_matrix_exp
-
-
- def func_exp_chi2(self, params):
+ return jacobian_matrix_exp
+
+
+ def func_exp_chi2(self, params=None, times=None, values=None,
errors=None):
"""Target function for minimising chi2 in minfx, for exponential fit.
@param params: The vector of parameter values.
@type params: numpy rank-1 float array
- @return: The chi-squared value.
+ @keyword times: The time points.
+ @type times: numpy array
+ @param values: The measured values.
+ @type values: numpy array
+ @param errors: The standard deviation of the measured intensity
values per time point.
+ @type errors: numpy array
+ @return: The chi2 value.
@rtype: float
"""
# Calculate.
- self.back_calc[:] = self.func_exp(params=params, times=self.times, )
+ back_calc = self.func_exp(params=params, times=times)
# Return the total chi-squared value.
- chi2 = chi2_rankN(data=self.values, back_calc_vals=self.back_calc,
errors=self.errors)
+ chi2 = chi2_rankN(data=values, back_calc_vals=back_calc,
errors=errors)
# Calculate and return the chi-squared value.
return chi2
- def func_exp_chi2_grad(self, params):
+ def func_exp_chi2_grad(self, params=None, times=None, values=None,
errors=None):
"""Target function for the gradient (Jacobian matrix) of
func_exp_chi2() to minfx, for exponential fit .
@param params: The vector of parameter values.
@type params: numpy rank-1 float array
- @return: The Jacobian matrix with 'm' rows of function
derivatives per 'n' columns of parameters.
+ @keyword times: The time points.
+ @type times: numpy array
+ @param values: The measured values.
+ @type values: numpy array
+ @param errors: The standard deviation of the measured intensity
values per time point.
+ @type errors: numpy array
+ @return: The Jacobian matrix with 'm' rows of function
derivatives per 'n' columns of parameters, which have been summed together.
@rtype: numpy array
"""
- # Unpack the parameter values.
- r2eff = params[0]
- i0 = params[1]
-
- # Calculate gradient according to chi2.
- # See:
http://wiki.nmr-relax.com/Calculate_jacobian_hessian_matrix_in_sympy_exponential_decay
-
- # Make partial derivative, with respect to r2eff.
- d_chi2_d_r2eff = 2.0 * i0 * self.times * ( -i0 * exp( -r2eff *
self.times) + self.values) * exp( -r2eff * self.times ) / self.errors**2
-
- # Make partial derivative, with respect to i0.
- d_chi2_d_i0 = - 2.0 * ( -i0 * exp( -r2eff * self.times) +
self.values) * exp( -r2eff * self.times) / self.errors**2
-
- # Define Jacobian as m rows with function derivatives and n columns
of parameters.
- self.jacobian_matrix_exp_chi2 = transpose(array( [d_chi2_d_r2eff ,
d_chi2_d_i0] ) )
+ # Get the Jacobian.
+ exp_chi2_grad = func_exp_chi2_grad(params=params, times=times,
values=values, errors=errors)
+
+ # Transpose back, to get rows.
+ exp_chi2_grad_t = transpose(exp_chi2_grad)
+
+ # Extract vectors:
+ d_chi2_d_r2eff = exp_chi2_grad_t[0]
+ d_chi2_d_i0 = exp_chi2_grad_t[1]
# Take the sum, to send to minfx.
sum_d_chi2_d_r2eff = sum( d_chi2_d_r2eff )
@@ -835,7 +837,7 @@
# If 'sigma'/erros describes one standard deviation errors of the
input data points. The estimated covariance in 'pcov' is based on these
values.
absolute_sigma = True
else:
- func = E.func_exp_general
+ func = E.func_exp_residual
absolute_sigma = False
# All args to function. Params are packed out through function, then
other parameters.
@@ -942,13 +944,17 @@
x0 = asarray( E.estimate_x0_exp(times=E.times, values=E.values) )
if E.c_code:
+ # Minimise with C code.
+
# Initialise the function to minimise.
scaling_list = [1.0, 1.0]
setup(num_params=len(x0), num_times=len(E.times), values=E.values,
sd=E.errors, relax_times=E.times, scaling_matrix=scaling_list)
+ # Define function to minimise for minfx.
t_func = func_wrapper
t_dfunc = dfunc_wrapper
t_d2func = d2func_wrapper
+ args=()
else:
# Minimise with minfx.
@@ -956,9 +962,11 @@
t_func = E.func_exp_chi2
t_dfunc = E.func_exp_chi2_grad
t_d2func = None
+ # All args to function. Params are packed out through function, then
other parameters.
+ args=(E.times, E.values, E.errors)
# Minimise.
- results_minfx = generic_minimise(func=t_func, dfunc=t_dfunc,
d2func=t_d2func, args=(), x0=x0, min_algor=E.min_algor,
min_options=E.min_options, func_tol=E.func_tol, grad_tol=E.grad_tol,
maxiter=E.max_iterations, A=E.A, b=E.b, full_output=True, print_flag=0)
+ results_minfx = generic_minimise(func=t_func, dfunc=t_dfunc,
d2func=t_d2func, args=args, x0=x0, min_algor=E.min_algor,
min_options=E.min_options, func_tol=E.func_tol, grad_tol=E.grad_tol,
maxiter=E.max_iterations, A=E.A, b=E.b, full_output=True, print_flag=0)
# Unpack
param_vector, chi2, iter_count, f_count, g_count, h_count, warning =
results_minfx
@@ -971,21 +979,13 @@
# Calculate the direct exponential Jacobian matrix from C code.
jacobian_matrix_exp = transpose(asarray( jacobian(param_vector) ) )
- # Compare with python code.
- #E.func_exp_grad(param_vector)
- #jacobian_matrix_exp2 = E.jacobian_matrix_exp
- #print jacobian_matrix_exp - jacobian_matrix_exp2
else:
if E.chi2_jacobian:
- # Call class, to store value.
- E.func_exp_chi2_grad(param_vector)
- jacobian_matrix_exp = E.jacobian_matrix_exp_chi2
+ # Use the chi2 Jacobian.
+ jacobian_matrix_exp = func_exp_chi2_grad(params=param_vector,
times=E.times, values=E.values, errors=E.errors)
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
+ # Use the direct Jacobian.
+ jacobian_matrix_exp = E.func_exp_grad(params=param_vector,
times=E.times, values=E.values, errors=E.errors)
# Get the co-variance
if E.chi2_jacobian:
@@ -1001,9 +1001,11 @@
eye_mat = eye(E.errors.shape[0])
# Get the residual vector K.
- K = E.func_exp_residual(params=param_vector, times=E.times,
values=E.values)
-
- weights = 1. / E.errors**2
+ #K = E.func_exp_residual(params=param_vector, times=E.times,
values=E.values)
+ #weights = 1. / E.errors**2
+
+ # Equal to:
+ K = E.func_exp_weighted_residual(params=param_vector, times=E.times,
values=E.values, errors=E.errors)
# Now form the weights matrix, with errors down the diagonal.
W = multiply(weights, eye_mat)
@@ -1021,7 +1023,6 @@
# Scale co-variance.
pcov = pcov * S02
-
# To compute one standard deviation errors on the parameters, take the
square root of the diagonal covariance.
param_vector_error = sqrt(diag(pcov))
r25404 - /trunk/specific_analyses/relax_disp/estimate_r2eff.py