Re: r25230 - in /trunk: target_functions/ test_suite/shared_data/curve_fitting/profiling/ -- August 25, 2014 - 16:42

Hi,

I've just created the target_functions.relax_disp.jacobian() function
at r25249 (http://article.gmane.org/gmane.science.nmr.relax.scm/22999).
This can be called with the optimised parameter vector as an argument
and then the covariance matrix calculated quite simply as described
in:

https://www.gnu.org/software/gsl/manual/html_node/Computing-the-covariance-matrix-of-best-fit-parameters.html

The equation is:

covar = (J^T J)^{-1}

A function to do this could be added to the relax library.  The input
argument is the Jacobian matrix, and the output is the covariance
matrix.  The QR decomposition in the above link is simply a faster
technique to find the inverse.

Regards,

Edward


P. S.  It shouldn't be too hard to extend the c_chi.c file to include
the chi-squared gradient to then have access to faster optimisation
algorithms.


On 25 August 2014 12:05, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:
> Yes, I am just after the error.
>
> Monte-Carlo simulations kill the "screening of data" of R1rho
> analysis, when R2eff error estimation is so slow.
>
> Best
> Troels
>
> 2014-08-25 11:44 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
> > If the issue of "fast" is not an issue, then why is an alternative
> > minimisation required?  Are you just after the covariance matrix?
> >
> > Regards,
> >
> > Edward
> >
> >
> > On 25 August 2014 11:16, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> 
> > wrote:
> > > Hi Edward.
> > >
> > > In this target function, exist the target function whereby
> > > scipy.optimize.leastsq minimises.
> > >
> > > minfx minimises by minimisation of chi2.
> > >
> > > scipy.optimize.leastsq minimises the return of a function which
> > > calculate the array of difference between
> > > function and measured values.
> > >
> > > It will then (I guess) minimise sum squares.
> > >
> > > You could probably make a call to the target function "exponential" 
> > > instead.
> > >
> > > But I am not so good with C code, and having the code in python makes
> > > it possible for me to bug-fix faster.
> > >
> > > The issue of "fast", does not really come into play here, as a
> > > minimisation takes about 0.02 seconds anyway.
> > >
> > > Best
> > > Troels
> > >
> > >
> > > 2014-08-25 10:33 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
> > > > Hi Troels,
> > > >
> > > > This new target_functions.relax_disp_curve_fit module clashes with the
> > > > C module.  It is repetitive, the code already exists in the
> > > > target_functions.relax_fit module.  Why is the faster C module not
> > > > being used?
> > > >
> > > > Regards,
> > > >
> > > > Edward
> > > >
> > > >
> > > >
> > > > On 25 August 2014 01:08,  <tlinnet@xxxxxxxxxxxxx> wrote:
> > > > > Author: tlinnet
> > > > > Date: Mon Aug 25 01:08:48 2014
> > > > > New Revision: 25230
> > > > >
> > > > > URL: http://svn.gna.org/viewcvs/relax?rev=25230&view=rev
> > > > > Log:
> > > > > Moved the target function for minimisation of exponential fit into the 
> > > > > target functions folder.
> > > > >
> > > > > task #7822(https://gna.org/task/index.php?7822): Implement user 
> > > > > function to estimate R2eff and associated errors for exponential curve 
> > > > > fitting.
> > > > >
> > > > > Added:
> > > > >     trunk/target_functions/relax_disp_curve_fit.py
> > > > >       - copied, changed from r25229, 
> > > > > trunk/test_suite/shared_data/curve_fitting/profiling/relax_fit.py
> > > > > Removed:
> > > > >     trunk/test_suite/shared_data/curve_fitting/profiling/relax_fit.py
> > > > > Modified:
> > > > >     trunk/target_functions/__init__.py
> > > > >     
> > > > > trunk/test_suite/shared_data/curve_fitting/profiling/profiling_relax_fit.py
> > > > >     trunk/test_suite/shared_data/curve_fitting/profiling/verify_error.py
> > > > >
> > > > > Modified: trunk/target_functions/__init__.py
> > > > > URL: 
> > > > > http://svn.gna.org/viewcvs/relax/trunk/target_functions/__init__.py?rev=25230&r1=25229&r2=25230&view=diff
> > > > > ==============================================================================
> > > > > --- trunk/target_functions/__init__.py  (original)
> > > > > +++ trunk/target_functions/__init__.py  Mon Aug 25 01:08:48 2014
> > > > > @@ -32,5 +32,6 @@
> > > > >      'mf',
> > > > >      'n_state_model',
> > > > >      'potential',
> > > > > -    'relax_disp'
> > > > > +    'relax_disp',
> > > > > +    'relax_disp_curve_fit'
> > > > >  ]
> > > > >
> > > > > Copied: trunk/target_functions/relax_disp_curve_fit.py (from r25229, 
> > > > > trunk/test_suite/shared_data/curve_fitting/profiling/relax_fit.py)
> > > > > URL: 
> > > > > http://svn.gna.org/viewcvs/relax/trunk/target_functions/relax_disp_curve_fit.py?p2=trunk/target_functions/relax_disp_curve_fit.py&p1=trunk/test_suite/shared_data/curve_fitting/profiling/relax_fit.py&r1=25229&r2=25230&rev=25230&view=diff
> > > > > ==============================================================================
> > > > >     (empty)
> > > > >
> > > > > Modified: 
> > > > > trunk/test_suite/shared_data/curve_fitting/profiling/profiling_relax_fit.py
> > > > > URL: 
> > > > > http://svn.gna.org/viewcvs/relax/trunk/test_suite/shared_data/curve_fitting/profiling/profiling_relax_fit.py?rev=25230&r1=25229&r2=25230&view=diff
> > > > > ==============================================================================
> > > > > --- 
> > > > > trunk/test_suite/shared_data/curve_fitting/profiling/profiling_relax_fit.py
> > > > >  (original)
> > > > > +++ 
> > > > > trunk/test_suite/shared_data/curve_fitting/profiling/profiling_relax_fit.py
> > > > >  Mon Aug 25 01:08:48 2014
> > > > > @@ -57,7 +57,7 @@
> > > > >  # relax module imports.
> > > > >  from status import Status; status = Status()
> > > > >  from target_functions.relax_fit import setup, func, dfunc, d2func, 
> > > > > back_calc_I
> > > > > -from relax_fit import Exponential
> > > > > +from target_functions.relax_disp_curve_fit import Exponential
> > > > >
> > > > >
> > > > >  # Alter setup.
> > > > >
> > > > > Removed: 
> > > > > trunk/test_suite/shared_data/curve_fitting/profiling/relax_fit.py
> > > > > URL: 
> > > > > http://svn.gna.org/viewcvs/relax/trunk/test_suite/shared_data/curve_fitting/profiling/relax_fit.py?rev=25229&view=auto
> > > > > ==============================================================================
> > > > > --- trunk/test_suite/shared_data/curve_fitting/profiling/relax_fit.py   
> > > > > (original)
> > > > > +++ trunk/test_suite/shared_data/curve_fitting/profiling/relax_fit.py   
> > > > > (removed)
> > > > > @@ -1,135 +0,0 @@
> > > > > -###############################################################################
> > > > > -#                                                                      
> > > > >        #
> > > > > -# Copyright (C) 2013-2014 Edward d'Auvergne                            
> > > > >        #
> > > > > -# Copyright (C) 2014 Troels E. Linnet                                  
> > > > >        #
> > > > > -#                                                                      
> > > > >        #
> > > > > -# This file is part of the program relax (http://www.nmr-relax.com).   
> > > > >        #
> > > > > -#                                                                      
> > > > >        #
> > > > > -# This program is free software: you can redistribute it and/or modify 
> > > > >        #
> > > > > -# it under the terms of the GNU General Public License as published by 
> > > > >        #
> > > > > -# the Free Software Foundation, either version 3 of the License, or    
> > > > >        #
> > > > > -# (at your option) any later version.                                  
> > > > >        #
> > > > > -#                                                                      
> > > > >        #
> > > > > -# This program is distributed in the hope that it will be useful,      
> > > > >        #
> > > > > -# but WITHOUT ANY WARRANTY; without even the implied warranty of       
> > > > >        #
> > > > > -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        
> > > > >        #
> > > > > -# GNU General Public License for more details.                         
> > > > >        #
> > > > > -#                                                                      
> > > > >        #
> > > > > -# You should have received a copy of the GNU General Public License    
> > > > >        #
> > > > > -# along with this program.  If not, see 
> > > > > <http://www.gnu.org/licenses/>.       #
> > > > > -#                                                                      
> > > > >        #
> > > > > -###############################################################################
> > > > > -
> > > > > -# Module docstring.
> > > > > -"""Target functions for relaxation fit."""
> > > > > -
> > > > > -# Python module imports.
> > > > > -from copy import deepcopy
> > > > > -from numpy import exp, multiply, sum
> > > > > -
> > > > > -# relax module imports.
> > > > > -
> > > > > -
> > > > > -class Exponential:
> > > > > -    def __init__(self, num_params=None, num_times=None, values=None, 
> > > > > sd=None, relax_times=None, scaling_matrix=None):
> > > > > -        """Relaxation dispersion target functions for optimisation.
> > > > > -        """
> > > > > -
> > > > > -        # Store variables.
> > > > > -        self.num_params = num_params
> > > > > -        self.num_times = num_times
> > > > > -        self.values = values
> > > > > -        self.errors = sd
> > > > > -        self.relax_times = relax_times
> > > > > -        self.scaling_matrix = scaling_matrix
> > > > > -
> > > > > -        # Create the structure for holding the back-calculated R2eff 
> > > > > values (matching the dimensions of the values structure).
> > > > > -        self.back_calc = deepcopy(self.values)
> > > > > -
> > > > > -        # Define function to minimise.
> > > > > -        self.func = self.func_exp
> > > > > -        self.calc = self.calc_exp
> > > > > -
> > > > > -
> > > > > -    def chi2_rankN(self, data, back_calc_vals, errors):
> > > > > -        """Function to calculate the chi-squared value for multiple 
> > > > > numpy array axis.
> > > > > -
> > > > > -        @param data:            The multi dimensional vectors of yi 
> > > > > values.
> > > > > -        @type data:             numpy multi dimensional array
> > > > > -        @param back_calc_vals:  The multi dimensional vectors of 
> > > > > yi(theta) values.
> > > > > -        @type back_calc_vals:   numpy multi dimensional array
> > > > > -        @param errors:          The multi dimensional vectors of 
> > > > > sigma_i values.
> > > > > -        @type errors:           numpy multi dimensional array
> > > > > -        @return:                The chi-squared value.
> > > > > -        @rtype:                 float
> > > > > -        """
> > > > > -
> > > > > -        # Calculate the chi-squared statistic.
> > > > > -        return sum((1.0 / errors * (data - back_calc_vals))**2)
> > > > > -
> > > > > -
> > > > > -    def calc_exp(self, times=None, r2eff=None, i0=None):
> > > > > -        """Calculate the function values of exponential function.
> > > > > -
> > > > > -        @keyword times: The time points.
> > > > > -        @type times:    float
> > > > > -        @keyword r2eff: The effective transversal relaxation rate.
> > > > > -        @type r2eff:    float
> > > > > -        @keyword i0:    The initial intensity.
> > > > > -        @type i0:       float
> > > > > -        @return:        The function values.
> > > > > -        @rtype:         float
> > > > > -        """
> > > > > -
> > > > > -        # Calculate.
> > > > > -        return i0 * exp( -times * r2eff)
> > > > > -
> > > > > -
> > > > > -    def calc_exp_chi2(self, r2eff=None, i0=None):
> > > > > -        """Calculate the chi-squared value of exponential function.
> > > > > -
> > > > > -
> > > > > -        @keyword r2eff: The effective transversal relaxation rate.
> > > > > -        @type r2eff:    float
> > > > > -        @keyword i0:    The initial intensity.
> > > > > -        @type i0:       float
> > > > > -        @return:        The chi-squared value.
> > > > > -        @rtype:         float
> > > > > -        """
> > > > > -
> > > > > -        # Calculate.
> > > > > -        self.back_calc[:] = self.calc_exp(times=self.relax_times, 
> > > > > r2eff=r2eff, i0=i0)
> > > > > -
> > > > > -        # Return the total chi-squared value.
> > > > > -        return self.chi2_rankN(data=self.values, 
> > > > > back_calc_vals=self.back_calc, errors=self.errors)
> > > > > -
> > > > > -
> > > > > -    def func_exp(self, params):
> > > > > -        """Target function for exponential fit.
> > > > > -
> > > > > -        @param params:  The vector of parameter values.
> > > > > -        @type params:   numpy rank-1 float array
> > > > > -        @return:        The chi-squared value.
> > > > > -        @rtype:         float
> > > > > -        """
> > > > > -
> > > > > -        # Unpack the parameter values.
> > > > > -        r2eff = params[0]
> > > > > -        i0 = params[1]
> > > > > -
> > > > > -        # Calculate and return the chi-squared value.
> > > > > -        return self.calc_exp_chi2(r2eff=r2eff, i0=i0)
> > > > > -
> > > > > -
> > > > > -    def func_exp_general(self, params, xdata, ydata):
> > > > > -        """Target function for minimisation with scipy.optimize.leastsq
> > > > > -        """
> > > > > -
> > > > > -        return self.calc_exp(xdata, *params) - ydata
> > > > > -
> > > > > -
> > > > > -    def func_exp_weighted_general(self, params, xdata, ydata, weights):
> > > > > -        """Target function for minimisation with scipy.optimize.leastsq
> > > > > -        """
> > > > > -
> > > > > -        return weights * (self.calc_exp(xdata, *params) - ydata)
> > > > >
> > > > > Modified: 
> > > > > trunk/test_suite/shared_data/curve_fitting/profiling/verify_error.py
> > > > > URL: 
> > > > > http://svn.gna.org/viewcvs/relax/trunk/test_suite/shared_data/curve_fitting/profiling/verify_error.py?rev=25230&r1=25229&r2=25230&view=diff
> > > > > ==============================================================================
> > > > > --- 
> > > > > trunk/test_suite/shared_data/curve_fitting/profiling/verify_error.py    
> > > > >     (original)
> > > > > +++ 
> > > > > trunk/test_suite/shared_data/curve_fitting/profiling/verify_error.py    
> > > > >     Mon Aug 25 01:08:48 2014
> > > > > @@ -32,7 +32,7 @@
> > > > >  from status import Status; status = Status()
> > > > >
> > > > >  # Initial try for Exponential class.
> > > > > -from relax_fit import Exponential
> > > > > +from target_functions.relax_disp_curve_fit import Exponential
> > > > >
> > > > >  # Define data path.
> > > > >  prev_data_path = status.install_path + 
> > > > > sep+'test_suite'+sep+'shared_data'+sep+'dispersion'+sep+'Kjaergaard_et_al_2013'
> > > > >  +sep+ "check_graphs" +sep+ "mc_2000"  +sep+ "R2eff"
> > > > >
> > > > >
> > > > > _______________________________________________
> > > > > 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