Hi Edward. I is just target_functions/c_chi2.c where you dont sum the elements, but return as array.
Best Troels 2014-08-28 15:30 GMT+02:00 Edward d'Auvergne <edw...@nmr-relax.com>: > Hi Troels, > > Could you derive the chi-squared Jacobian? Maybe the Jacobian I have > been using is not correct - this is the one required for the > Levenberg-Marquardt optimisation algorithm. Because the chi-squared > is squared, its derivative will have a factor of 2 out the front, just > like the gradient: > > http://www.nmr-relax.com/manual/chi_squared_gradient.html > > It might be useful to add a Jacobian section to this part of the > manual with the equations. > > Cheers, > > Edward > > > > On 28 August 2014 15:14, <tlin...@nmr-relax.com> wrote: >> Author: tlinnet >> Date: Thu Aug 28 15:14:16 2014 >> New Revision: 25379 >> >> URL: http://svn.gna.org/viewcvs/relax?rev=25379&view=rev >> Log: >> Modified systemtest test Relax_disp.test_estimate_r2eff_err_methods() to >> show the difference between using the direct function Jacobian, or the chi2 >> function Jacobian. >> >> Added also the functionality to the estimate R2eff module, to switch between >> using the different Jacobians. >> >> The results show, that R2eff can be estimated better. >> >> ---------------------- >> The results are: >> >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 431.0. >> r2eff=8.646/8.646 r2eff_err=0.0348/0.0692 i0=202664.191/202664.191 >> i0_err=699.6443/712.4201 >> >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 651.2. >> r2eff=10.377/10.377 r2eff_err=0.0403/0.0810 i0=206049.558/206049.558 >> i0_err=776.4215/782.1833 >> >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 800.5. >> r2eff=10.506/10.506 r2eff_err=0.0440/0.0853 i0=202586.332/202586.332 >> i0_err=763.9678/758.7052 >> >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 984.0. >> r2eff=10.903/10.903 r2eff_err=0.0476/0.0922 i0=203455.021/203455.021 >> i0_err=837.8779/828.7280 >> >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1341.1. >> r2eff=10.684/10.684 r2eff_err=0.0446/0.0853 i0=218670.412/218670.412 >> i0_err=850.0210/830.9558 >> >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1648.5. >> r2eff=10.501/10.501 r2eff_err=0.0371/0.0742 i0=206502.512/206502.512 >> i0_err=794.0523/772.9843 >> >> R1rho at 799.8 MHz, for offset=124.247 ppm and dispersion point 1341.1. >> r2eff=11.118/11.118 r2eff_err=0.0413/0.0827 i0=216447.241/216447.241 >> i0_err=784.6562/788.0384 >> >> R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 800.5. >> r2eff=7.866/7.866 r2eff_err=0.0347/0.0695 i0=211869.715/211869.715 >> i0_err=749.2776/763.6930 >> >> R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1341.1. >> r2eff=9.259/9.259 r2eff_err=0.0331/0.0661 i0=217703.151/217703.151 >> i0_err=682.2137/685.5838 >> >> R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1648.5. >> r2eff=9.565/9.565 r2eff_err=0.0373/0.0745 i0=211988.939/211988.939 >> i0_err=839.0313/827.0373 >> >> R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 800.5. >> r2eff=3.240/3.240 r2eff_err=0.0127/0.0253 i0=214417.382/214417.382 >> i0_err=595.8865/613.4378 >> >> R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 1341.1. >> r2eff=5.084/5.084 r2eff_err=0.0177/0.0352 i0=226358.691/226358.691 >> i0_err=660.5314/655.7670 >> >> R1rho at 799.8 MHz, for offset=179.768 ppm and dispersion point 1341.1. >> r2eff=2.208/2.208 r2eff_err=0.0091/0.0178 i0=228620.553/228620.553 >> i0_err=564.8353/560.0873 >> >> R1rho at 799.8 MHz, for offset=241.459 ppm and dispersion point 1341.1. >> r2eff=1.711/1.711 r2eff_err=0.0077/0.0155 i0=224087.486/224087.486 >> i0_err=539.4300/546.4217 >> >> Fitting with minfx to: 52V @N >> ----------------------------- >> >> min_algor='Newton', c_code=True, constraints=False, chi2_jacobian?=False >> ------------------------------------------------------------------------ >> >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 431.0, with >> 4 time points. r2eff=8.646 r2eff_err=0.0692, i0=202664.2, i0_err=712.4201, >> chi2=3.758. >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 651.2, with >> 5 time points. r2eff=10.377 r2eff_err=0.0810, i0=206049.6, i0_err=782.1833, >> chi2=27.291. >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 800.5, with >> 5 time points. r2eff=10.506 r2eff_err=0.0853, i0=202586.3, i0_err=758.7052, >> chi2=13.357. >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 984.0, with >> 5 time points. r2eff=10.903 r2eff_err=0.0922, i0=203455.0, i0_err=828.7280, >> chi2=33.632. >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1341.1, with >> 5 time points. r2eff=10.684 r2eff_err=0.0853, i0=218670.4, i0_err=830.9558, >> chi2=35.818. >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1648.5, with >> 5 time points. r2eff=10.501 r2eff_err=0.0742, i0=206502.5, i0_err=772.9843, >> chi2=7.356. >> R1rho at 799.8 MHz, for offset=124.247 ppm and dispersion point 1341.1, with >> 5 time points. r2eff=11.118 r2eff_err=0.0827, i0=216447.2, i0_err=788.0384, >> chi2=15.587. >> R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 800.5, with >> 5 time points. r2eff=7.866 r2eff_err=0.0695, i0=211869.7, i0_err=763.6930, >> chi2=14.585. >> R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1341.1, with >> 5 time points. r2eff=9.259 r2eff_err=0.0661, i0=217703.2, i0_err=685.5838, >> chi2=79.498. >> R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1648.5, with >> 5 time points. r2eff=9.565 r2eff_err=0.0745, i0=211988.9, i0_err=827.0373, >> chi2=0.447. >> R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 800.5, with >> 5 time points. r2eff=3.240 r2eff_err=0.0253, i0=214417.4, i0_err=613.4378, >> chi2=1.681. >> R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 1341.1, with >> 5 time points. r2eff=5.084 r2eff_err=0.0352, i0=226358.7, i0_err=655.7670, >> chi2=23.170. >> R1rho at 799.8 MHz, for offset=179.768 ppm and dispersion point 1341.1, with >> 5 time points. r2eff=2.208 r2eff_err=0.0178, i0=228620.6, i0_err=560.0873, >> chi2=7.794. >> R1rho at 799.8 MHz, for offset=241.459 ppm and dispersion point 1341.1, with >> 5 time points. r2eff=1.711 r2eff_err=0.0155, i0=224087.5, i0_err=546.4217, >> chi2=21.230. >> >> Fitting with minfx to: 52V @N >> ----------------------------- >> >> min_algor='BFGS', c_code=False, constraints=False, chi2_jacobian?=True >> ---------------------------------------------------------------------- >> >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 431.0, with >> 4 time points. r2eff=8.646 r2eff_err=0.0524, i0=202664.2, i0_err=1239.0827, >> chi2=3.758. >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 651.2, with >> 5 time points. r2eff=10.377 r2eff_err=0.0228, i0=206049.6, i0_err=178.1907, >> chi2=27.291. >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 800.5, with >> 5 time points. r2eff=10.506 r2eff_err=0.0345, i0=202586.3, i0_err=705.7630, >> chi2=13.357. >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 984.0, with >> 5 time points. r2eff=10.903 r2eff_err=0.0206, i0=203455.0, i0_err=186.0857, >> chi2=33.632. >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1341.1, with >> 5 time points. r2eff=10.684 r2eff_err=0.0198, i0=218670.4, i0_err=165.0420, >> chi2=35.818. >> R1rho at 799.8 MHz, for offset=118.078 ppm and dispersion point 1648.5, with >> 5 time points. r2eff=10.501 r2eff_err=0.0407, i0=206502.5, i0_err=321.3685, >> chi2=7.356. >> R1rho at 799.8 MHz, for offset=124.247 ppm and dispersion point 1341.1, with >> 5 time points. r2eff=11.118 r2eff_err=0.0301, i0=216447.2, i0_err=248.9394, >> chi2=15.587. >> R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 800.5, with >> 5 time points. r2eff=7.866 r2eff_err=0.0280, i0=211869.7, i0_err=259.8845, >> chi2=14.585. >> R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1341.1, with >> 5 time points. r2eff=9.259 r2eff_err=0.0108, i0=217703.2, i0_err=88.1514, >> chi2=79.498. >> R1rho at 799.8 MHz, for offset=130.416 ppm and dispersion point 1648.5, with >> 5 time points. r2eff=9.565 r2eff_err=0.1630, i0=211988.9, i0_err=2054.6615, >> chi2=0.447. >> R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 800.5, with >> 5 time points. r2eff=3.240 r2eff_err=0.0485, i0=214417.4, i0_err=611.7573, >> chi2=1.681. >> R1rho at 799.8 MHz, for offset=142.754 ppm and dispersion point 1341.1, with >> 5 time points. r2eff=5.084 r2eff_err=0.0124, i0=226358.7, i0_err=122.7341, >> chi2=23.170. >> R1rho at 799.8 MHz, for offset=179.768 ppm and dispersion point 1341.1, with >> 5 time points. r2eff=2.208 r2eff_err=0.0086, i0=228620.6, i0_err=219.4208, >> chi2=7.794. >> R1rho at 799.8 MHz, for offset=241.459 ppm and dispersion point 1341.1, with >> 5 time points. r2eff=1.711 r2eff_err=0.0101, i0=224087.5, i0_err=166.9081, >> chi2=21.230. >> >> 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 >> trunk/test_suite/system_tests/relax_disp.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=25379&r1=25378&r2=25379&view=diff >> ============================================================================== >> --- trunk/specific_analyses/relax_disp/estimate_r2eff.py (original) >> +++ trunk/specific_analyses/relax_disp/estimate_r2eff.py Thu Aug 28 >> 15:14:16 2014 >> @@ -175,7 +175,7 @@ >> print(print_string), >> >> >> -def multifit_covar(J=None, epsrel=0.0, errors=None): >> +def multifit_covar(J=None, epsrel=0.0, errors=None, use_weights=True): >> """This is the implementation of the multifit covariance. >> >> This is inspired from GNU Scientific Library (GSL). >> @@ -184,9 +184,15 @@ >> >> The parameter 'epsrel' is used to remove linear-dependent columns when >> J is rank deficient. >> >> + The weighting matrix 'W', is a square symmetric matrix. For independent >> measurements, this is a diagonal matrix. Larger values indicate greater >> significance. It is formed by multiplying the supplied errors as >> 1./errors^2 with an Identity matrix:: >> + >> + W = I.(1/errors^2) >> + >> + If 'use_weights' is set to 'False', the errors are set to 1.0. >> + >> The covariance matrix is given by:: >> >> - covar = (J^T J)^{-1} , >> + covar = (J^T.W.J)^{-1} , >> >> and is computed by QR decomposition of J with column-pivoting. Any >> columns of R which satisfy:: >> >> @@ -224,6 +230,8 @@ >> @type epsrel: float >> @keyword errors: The standard deviation of the measured >> intensity values per time point. >> @type errors: numpy array >> + @keyword use_weights: If the supplied weights should be used. >> + @type use_weights: bool >> @return: The co-variance matrix >> @rtype: square numpy array >> """ >> @@ -237,6 +245,10 @@ >> # Now form the error matrix, with errors down the diagonal. >> weights = 1. / errors**2 >> >> + if use_weights == False: >> + weights[:] = 1.0 >> + >> + # Form weight matrix. >> W = multiply(weights, eye_mat) >> >> # The covariance matrix (sometimes referred to as the >> variance-covariance matrix), Qxx, is defined as: >> @@ -344,7 +356,7 @@ >> self.factor = factor >> >> >> - def set_settings_minfx(self, scaling_matrix=None, min_algor='simplex', >> c_code=True, constraints=False, func_tol=1e-25, grad_tol=None, >> max_iterations=10000000): >> + def set_settings_minfx(self, scaling_matrix=None, min_algor='simplex', >> c_code=True, constraints=False, chi2_jacobian=False, func_tol=1e-25, >> grad_tol=None, max_iterations=10000000): >> """Setup options to minfx. >> >> @keyword scaling_matrix: The square and diagonal scaling matrix. >> @@ -355,6 +367,8 @@ >> @type c_code: bool >> @keyword constraints: If constraints should be used. >> @type constraints: bool >> + @keyword chi2_jacobian: If the chi2 Jacobian should be used. >> + @type chi2_jacobian: bool >> @keyword func_tol: The function tolerance which, when >> reached, terminates optimisation. Setting this to None turns of the check. >> @type func_tol: None or float >> @keyword grad_tol: The gradient tolerance which, when >> reached, terminates optimisation. Setting this to None turns of the check. >> @@ -366,6 +380,7 @@ >> # Store variables. >> self.scaling_matrix = scaling_matrix >> self.c_code = c_code >> + self.chi2_jacobian = chi2_jacobian >> >> # Scaling initialisation. >> self.scaling_flag = False >> @@ -561,7 +576,7 @@ >> return 1. / self.errors * (self.func_exp(self.times, *params) - >> self.values) >> >> >> -def estimate_r2eff(method='minfx', min_algor='simplex', c_code=True, >> constraints=False, spin_id=None, ftol=1e-15, xtol=1e-15, maxfev=10000000, >> factor=100.0, verbosity=1): >> +def estimate_r2eff(method='minfx', min_algor='simplex', c_code=True, >> constraints=False, chi2_jacobian=False, spin_id=None, ftol=1e-15, >> xtol=1e-15, maxfev=10000000, factor=100.0, verbosity=1): >> """Estimate r2eff and errors by exponential curve fitting with >> scipy.optimize.leastsq or minfx. >> >> THIS IS ONLY FOR TESTING. >> @@ -583,10 +598,12 @@ >> @type method: string >> @keyword min_algor: The minimisation algorithm >> @type min_algor: string >> + @keyword c_code: If optimise with C code. >> + @type c_code: bool >> @keyword constraints: If constraints should be used. >> @type constraints: bool >> - @keyword c_code: If optimise with C code. >> - @type c_code: bool >> + @keyword chi2_jacobian: If the chi2 Jacobian should be used. >> + @type chi2_jacobian: bool >> @keyword spin_id: The spin identification string. >> @type spin_id: str >> @keyword ftol: The function tolerance for the relative >> error desired in the sum of squares, parsed to leastsq. >> @@ -661,7 +678,7 @@ >> top += 2 >> subsection(file=sys.stdout, text="Fitting with %s to: >> %s"%(method, spin_string), prespace=top) >> if method == 'minfx': >> - subsection(file=sys.stdout, text="min_algor='%s', >> c_code=%s, constraints=%s"%(min_algor, c_code, constraints), prespace=0) >> + subsection(file=sys.stdout, text="min_algor='%s', >> c_code=%s, constraints=%s, chi2_jacobian?=%s"%(min_algor, c_code, >> constraints, chi2_jacobian), prespace=0) >> >> # Loop over each spectrometer frequency and dispersion point. >> for exp_type, frq, offset, point, ei, mi, oi, di in >> loop_exp_frq_offset_point(return_indices=True): >> @@ -692,7 +709,7 @@ >> >> elif method == 'minfx': >> # Set settings. >> - E.set_settings_minfx(min_algor=min_algor, c_code=c_code, >> constraints=constraints) >> + E.set_settings_minfx(min_algor=min_algor, c_code=c_code, >> chi2_jacobian=chi2_jacobian, constraints=constraints) >> >> # Acquire results. >> results = minimise_minfx(E=E) >> @@ -737,7 +754,7 @@ >> point_info = "%s at %3.1f MHz, for offset=%3.3f ppm and >> dispersion point %-5.1f, with %i time points." % (exp_type, frq/1E6, offset, >> point, len(times)) >> print_strings.append(point_info) >> >> - par_info = "r2eff=%3.3f r2eff_err=%3.3f, i0=%6.1f, >> i0_err=%3.3f, chi2=%3.3f.\n" % ( r2eff, r2eff_err, i0, i0_err, chi2) >> + par_info = "r2eff=%3.3f r2eff_err=%3.4f, i0=%6.1f, >> i0_err=%3.4f, chi2=%3.3f.\n" % ( r2eff, r2eff_err, i0, i0_err, chi2) >> print_strings.append(par_info) >> >> if E.verbosity >= 2: >> @@ -912,14 +929,24 @@ >> #jacobian_matrix_exp2 = E.jacobian_matrix_exp >> #print jacobian_matrix_exp - jacobian_matrix_exp2 >> else: >> - # Call class, to store value. >> - E.func_exp_grad(param_vector) >> - jacobian_matrix_exp = E.jacobian_matrix_exp >> - #E.func_exp_chi2_grad(param_vector) >> - #jacobian_matrix_exp = E.jacobian_matrix_exp_chi2 >> + if E.chi2_jacobian: >> + # Call class, to store value. >> + E.func_exp_chi2_grad(param_vector) >> + jacobian_matrix_exp = E.jacobian_matrix_exp_chi2 >> + else: >> + # Call class, to store value. >> + E.func_exp_grad(param_vector) >> + jacobian_matrix_exp = E.jacobian_matrix_exp >> + #E.func_exp_chi2_grad(param_vector) >> + #jacobian_matrix_exp = E.jacobian_matrix_exp_chi2 >> >> # Get the co-variance >> - pcov = multifit_covar(J=jacobian_matrix_exp, errors=E.errors) >> + if E.chi2_jacobian: >> + use_weights = False >> + else: >> + use_weights = True >> + >> + pcov = multifit_covar(J=jacobian_matrix_exp, errors=E.errors, >> use_weights=use_weights) >> >> # To compute one standard deviation errors on the parameters, take the >> square root of the diagonal covariance. >> param_vector_error = sqrt(diag(pcov)) >> >> Modified: trunk/test_suite/system_tests/relax_disp.py >> URL: >> http://svn.gna.org/viewcvs/relax/trunk/test_suite/system_tests/relax_disp.py?rev=25379&r1=25378&r2=25379&view=diff >> ============================================================================== >> --- trunk/test_suite/system_tests/relax_disp.py (original) >> +++ trunk/test_suite/system_tests/relax_disp.py Thu Aug 28 15:14:16 2014 >> @@ -2946,12 +2946,13 @@ >> >> >> # Now do it manually. >> - estimate_r2eff(method='scipy.optimize.leastsq') >> - estimate_r2eff(method='minfx', min_algor='simplex', c_code=True, >> constraints=False) >> - estimate_r2eff(method='minfx', min_algor='simplex', c_code=False, >> constraints=False) >> - estimate_r2eff(method='minfx', min_algor='BFGS', c_code=True, >> constraints=False) >> - estimate_r2eff(method='minfx', min_algor='BFGS', c_code=False, >> constraints=False) >> + #estimate_r2eff(method='scipy.optimize.leastsq') >> + #estimate_r2eff(method='minfx', min_algor='simplex', c_code=True, >> constraints=False) >> + #estimate_r2eff(method='minfx', min_algor='simplex', c_code=False, >> constraints=False) >> + #estimate_r2eff(method='minfx', min_algor='BFGS', c_code=True, >> constraints=False) >> + #estimate_r2eff(method='minfx', min_algor='BFGS', c_code=False, >> constraints=False) >> estimate_r2eff(method='minfx', min_algor='Newton', c_code=True, >> constraints=False) >> + estimate_r2eff(method='minfx', min_algor='BFGS', c_code=False, >> constraints=False, chi2_jacobian=True) >> >> >> def test_exp_fit(self): >> >> >> _______________________________________________ >> relax (http://www.nmr-relax.com) >> >> This is the relax-commits mailing list >> relax-comm...@gna.org >> >> 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@gna.org > > 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 _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-devel mailing list relax-devel@gna.org 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
