mailr25330 - /trunk/specific_analyses/relax_disp/estimate_r2eff.py


Others Months | Index by Date | Thread Index
>>   [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Header


Content

Posted by tlinnet on August 27, 2014 - 11:29:
Author: tlinnet
Date: Wed Aug 27 11:29:24 2014
New Revision: 25330

URL: http://svn.gna.org/viewcvs/relax?rev=25330&view=rev
Log:
Tried to implement the Jacobian from C-code.

This though also report errors which are to small.

Maybe some scaling is wrong.

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=25330&r1=25329&r2=25330&view=diff
==============================================================================
--- trunk/specific_analyses/relax_disp/estimate_r2eff.py        (original)
+++ trunk/specific_analyses/relax_disp/estimate_r2eff.py        Wed Aug 27 
11:29:24 2014
@@ -42,7 +42,7 @@
 from specific_analyses.relax_disp.variables import MODEL_R2EFF
 from specific_analyses.relax_fit.optimisation import func_wrapper, 
dfunc_wrapper, d2func_wrapper
 from target_functions.chi2 import chi2_rankN
-from target_functions.relax_fit import setup
+from target_functions.relax_fit import jacobian, setup
 
 
 # Scipy installed.
@@ -734,7 +734,7 @@
     E.set_settings_minfx(min_algor=min_algor)
 
     # Do C code
-    do_C = False
+    do_C = True
 
     if do_C:
         # Initialise the function to minimise.
@@ -766,19 +766,27 @@
     param_vector, chi2, iter_count, f_count, g_count, h_count, warning = 
results_minfx
 
     # Get the Jacobian.
-    # First make a call to the Jacobian function, which store it in the 
class.
-    E.func_exp_grad(params=param_vector)
-    jacobian_matrix = deepcopy(E.jacobian_matrix)
-
+    if do_C:
+        # First make a call to the Jacobian function, which store it in the 
class.
+        jacobian_matrix = transpose(asarray( jacobian(param_vector) ) )
+
+        # Compare with python code.
+        #E.func_exp_grad(params=param_vector)
+        #jacobian_matrix2 = deepcopy(E.jacobian_matrix)
+        #print jacobian_matrix
+        #print " "
+        #print jacobian_matrix2
+    else:
+        jacobian_matrix = deepcopy(E.jacobian_matrix)
+
+    # Get the co-variance
+    pcov = E.multifit_covar(J=jacobian_matrix)
+
+    # To compute one standard deviation errors on the parameters, take the 
square root of the diagonal covariance.
+    param_vector_error = sqrt(diag(pcov))
     # Set error to inf.
     #param_vector_error = [inf, inf]
 
-    # Get the co-variance
-    pcov = E.multifit_covar(J=jacobian_matrix)
-
-    # To compute one standard deviation errors on the parameters, take the 
square root of the diagonal covariance.
-    param_vector_error = sqrt(diag(pcov))
-
     # Pack to list.
     results = [param_vector, param_vector_error, chi2, iter_count, f_count, 
g_count, h_count, warning]
 




Related Messages


Powered by MHonArc, Updated Wed Aug 27 12:00:03 2014