r25255 - /trunk/specific_analyses/relax_disp/estimate_r2eff.py -- August 25, 2014 - 20:57

Author: tlinnet
Date: Mon Aug 25 20:57:15 2014
New Revision: 25255

URL: http://svn.gna.org/viewcvs/relax?rev=25255&view=rev
Log:
Tried algorithms:

newton
Newton-CG
Steepest descent
Fletcher-Reeves

None of them work.

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=25255&r1=25254&r2=25255&view=diff
==============================================================================
--- trunk/specific_analyses/relax_disp/estimate_r2eff.py        (original)
+++ trunk/specific_analyses/relax_disp/estimate_r2eff.py        Mon Aug 25 
20:57:15 2014
@@ -26,6 +26,7 @@
 from copy import deepcopy
 from numpy import asarray, array, diag, dot, exp, inf, log, sqrt, sum, 
transpose, zeros
 from minfx.generic import generic_minimise
+from re import match, search
 import sys
 from warnings import warn
 
@@ -113,7 +114,19 @@
 
         else:
             #self.min_algor = 'simplex'
-            self.min_algor = 'newton'
+
+            # Newton does not work.
+            #self.min_algor = 'newton'
+
+            # Newton-CG does not work.
+            self.min_algor = 'Newton-CG'
+
+            # Also not work.
+            #self.min_algor = 'Steepest descent'
+
+            # Also not work.#
+            #self.min_algor = 'Fletcher-Reeves'
+
             self.min_options = ()
             self.A = None
             self.b = None
@@ -122,7 +135,11 @@
         self.back_calc = deepcopy(self.values)
 
         # Define function to minimise for minfx.
-        self.func = self.func_exp
+        if match('^[Ss]implex$', self.min_algor):
+            self.func = self.func_exp
+        else:
+            self.func = self.func_exp_chi2
+
         self.dfunc = self.func_exp_grad
         self.d2func = self.func_exp_hess
 
@@ -213,6 +230,29 @@
         return self.calc_exp_chi2(r2eff=r2eff, i0=i0)
 
 
+    def func_exp_chi2(self, params):
+        """Target function for exponential fit in minfx.
+
+        @param params:  The vector of parameter values.
+        @type params:   numpy rank-1 float array
+        @return:        The chi-squared value.
+        @rtype:         float
+        """
+
+        # Scaling.
+        if self.scaling_flag:
+            params = dot(params, self.scaling_matrix)
+
+        # Unpack the parameter values.
+        r2eff = params[0]
+        i0 = params[1]
+
+        chi2 = 1.0 * ( -i0 * exp( -r2eff * self.relax_times) + 
self.values)**2 / self.errors**2
+
+        # Calculate and return the chi-squared value.
+        return chi2
+
+
     def func_exp_grad(self, params):
         """Target function for the gradient (Jacobian matrix) for 
exponential fit in minfx.
 
@@ -244,6 +284,8 @@
         # 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] ) )
 
+        #print jacobian_matrix
+
         # Return Jacobian matrix.
         return jacobian_matrix
 
@@ -283,6 +325,8 @@
 
         # 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] ) )
+
+        #print hessian_matrix
 
         # Return Jacobian matrix.
         return hessian_matrix