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]