Re: r25254 - /trunk/specific_analyses/relax_disp/estimate_r2eff.py -- August 26, 2014 - 09:50

Hi Troels,

This is not quite right, the Jacobian, gradient and Hessian are
different constructs.  I'll use the model-free case as an example as I
have fully documented this in the relax manual.  The Jacobian is:

http://www.nmr-relax.com/manual/Ri_theta_gradients.html
http://www.nmr-relax.com/manual/Ri_prime_theta_gradients.html

For model-free, instead of {intensity, time} points the data is {Ri,
None}, i.e. the relaxation data which has no X-axis correspondence.
So it is slightly different.  Anyway, these partial derivatives, when
in an array, form the Jacobian matrix.  The gradient is different -
you need the chi-squared gradient:

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

This equation is the same for the exponential curve-fitting, just that
Ri is replaced by the intensities I and the sum is over the time
points.  The last part of this equation are the Jacobian matrix
elements.  I am implementing this slowly in the C modules, but it
still needs some debugging.  If you wish to use this code, feel free
to help look and see where the current bug might be.

I will not code the Hessian into the C module, but you are free to
copy the steps I am performing to implement these yourself.  The
calculation of d2_chi2_d_r2eff_d_r2eff (which in the current relax
notation would be d2chi2_dr2eff2), would go into the
target_functions/exponential.c file as the exponential_dr2eff2()
function.  Note that this is not the chi2 Hessian element
d2chi2_dr2eff2, but instead the exponential curve Hessian (which is
probably a second order Jacobian).  You need to translate this to the
chi2 Hessian via:

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

See the target_functions/chi2.py file for the implementation.  The
final implementation would be a translation of this into the c_chi2.c
file.

Regards,

Edward



On 25 August 2014 20:03,  <tlinnet@xxxxxxxxxxxxx> wrote:
> Author: tlinnet
> Date: Mon Aug 25 20:03:40 2014
> New Revision: 25254
>
> URL: http://svn.gna.org/viewcvs/relax?rev=25254&view=rev
> Log:
> Initial try to form the Jacobian and Hessian matrix for exponential decay.
>
> This can be tried with systemtest:
> relax -s Relax_disp.test_estimate_r2eff_error
>
> 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=25254&r1=25253&r2=25254&view=diff
> ==============================================================================
> --- trunk/specific_analyses/relax_disp/estimate_r2eff.py        (original)
> +++ trunk/specific_analyses/relax_disp/estimate_r2eff.py        Mon Aug 25 
> 20:03:40 2014
> @@ -24,7 +24,7 @@
>
>  # Python module imports.
>  from copy import deepcopy
> -from numpy import asarray, diag, dot, exp, inf, log, sqrt, sum, zeros
> +from numpy import asarray, array, diag, dot, exp, inf, log, sqrt, sum, 
> transpose, zeros
>  from minfx.generic import generic_minimise
>  import sys
>  from warnings import warn
> @@ -112,7 +112,8 @@
>              self.b = array([   0., -200.,    0.])
>
>          else:
> -            self.min_algor = 'simplex'
> +            #self.min_algor = 'simplex'
> +            self.min_algor = 'newton'
>              self.min_options = ()
>              self.A = None
>              self.b = None
> @@ -122,6 +123,8 @@
>
>          # Define function to minimise for minfx.
>          self.func = self.func_exp
> +        self.dfunc = self.func_exp_grad
> +        self.d2func = self.func_exp_hess
>
>
>      def calc_exp(self, times=None, r2eff=None, i0=None):
> @@ -208,6 +211,81 @@
>
>          # Calculate and return the chi-squared value.
>          return self.calc_exp_chi2(r2eff=r2eff, i0=i0)
> +
> +
> +    def func_exp_grad(self, params):
> +        """Target function for the gradient (Jacobian matrix) for 
> exponential fit in minfx.
> +
> +        @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.
> +        @rtype:         numpy array
> +        """
> +
> +        # Scaling.
> +        if self.scaling_flag:
> +            params = dot(params, self.scaling_matrix)
> +
> +        # 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*times*(-i0*exp(-r2eff*times) + 
> values)*exp(-r2eff*times)/errors**2
> +        d_chi2_d_r2eff = 2.0 * i0 * self.relax_times * ( -i0 * exp( -r2eff 
> * self.relax_times) + self.values) * exp( -r2eff * self.relax_times ) / 
> self.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
> +        d_chi2_d_i0 = - 2.0 * ( -i0 * exp( -r2eff * self.relax_times) + 
> self.values) * exp( -r2eff * self.relax_times) / self.errors**2
> +
> +        # Define Jacobian as m rows with function derivatives and n 
> columns of parameters.
> +        jacobian_matrix = transpose(array( [d_chi2_d_r2eff , d_chi2_d_i0] 
> ) )
> +
> +        # Return Jacobian matrix.
> +        return jacobian_matrix
> +
> +
> +    def func_exp_hess(self, params):
> +        """Target function for the gradient (Hessian matrix) for 
> exponential fit in minfx.
> +
> +        @param params:  The vector of parameter values.
> +        @type params:   numpy rank-1 float array
> +        @return:        The Hessian matrix with 'm' rows of function 
> derivatives per '4' columns of second order derivatives.
> +        @rtype:         numpy array
> +        """
> +
> +        # Scaling.
> +        if self.scaling_flag:
> +            params = dot(params, self.scaling_matrix)
> +
> +        # 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
> +
> +        # Calculate the Hessian. The second-order partial derivatives.
> +        #d2_chi2_d_r2eff_d_r2eff = 2.0*i0*times**2*(2*i0*exp(-r2eff*times) 
> - values)*exp(-r2eff*times)/errors**2
> +        d2_chi2_d_r2eff_d_r2eff = 2.0 * i0 * self.relax_times**2 * ( 2 * 
> i0 * exp( -r2eff * self.relax_times) - self.values) * exp( -r2eff * 
> self.relax_times) / self.errors**2
> +
> +        #d2_chi2_d_r2eff_d_i0 = -2.0*times*(2*i0*exp(-r2eff*times) - 
> values)*exp(-r2eff*times)/errors**2
> +        d2_chi2_d_r2eff_d_i0 = -2.0 * self.relax_times * (2. * i0 * exp( 
> -r2eff * self.relax_times) - self.values) * exp( -r2eff * self.relax_times) 
> / self.errors**2
> +
> +        #d2_chi2_d_i0_d_r2eff = -2.0*times*(2*i0*exp(-r2eff*times) - 
> values)*exp(-r2eff*times)/errors**2
> +        d2_chi2_d_i0_d_r2eff = -2.0 * self.relax_times * (2. * i0 * exp( 
> -r2eff * self.relax_times) - self.values) * exp( -r2eff * self.relax_times) 
> / self.errors**2
> +
> +        #d2_chi2_d_i0_d_i0 = 2.0*exp(-2*r2eff*times)/errors**2
> +        d2_chi2_d_i0_d_i0 = 2.0 * exp( -2. * r2eff * self.relax_times) / 
> self.errors**2
> +
> +        # Form hessian.
> +        hessian_matrix = transpose(array( [d2_chi2_d_r2eff_d_r2eff, 
> d2_chi2_d_r2eff_d_i0, d2_chi2_d_i0_d_r2eff, d2_chi2_d_i0_d_i0] ) )
> +
> +        # Return Jacobian matrix.
> +        return hessian_matrix
>
>
>      def func_exp_general(self, params, times, intensities):
> @@ -544,7 +622,7 @@
>      x0 = asarray(exp_class.estimate_x0_exp(intensities=values, 
> times=times))
>
>      # Minimise with minfx.
> -    results_minfx = generic_minimise(func=exp_class.func, args=(), x0=x0, 
> min_algor=exp_class.min_algor, min_options=exp_class.min_options, 
> func_tol=exp_class.func_tol, grad_tol=exp_class.grad_tol, 
> maxiter=exp_class.max_iterations, A=exp_class.A, b=exp_class.b, 
> full_output=True, print_flag=exp_class.verbosity)
> +    results_minfx = generic_minimise(func=exp_class.func, 
> dfunc=exp_class.dfunc, d2func=exp_class.d2func, args=(), x0=x0, 
> min_algor=exp_class.min_algor, min_options=exp_class.min_options, 
> func_tol=exp_class.func_tol, grad_tol=exp_class.grad_tol, 
> maxiter=exp_class.max_iterations, A=exp_class.A, b=exp_class.b, 
> full_output=True, print_flag=exp_class.verbosity)
>
>      # Unpack
>      param_vector, chi2, iter_count, f_count, g_count, h_count, warning = 
> results_minfx
>
>
> _______________________________________________
> 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