Anyway, before minfx can handle constraints in for example BFGS, this is just a waste of time. I think there will be a 10 x speed up, just for the Jacobian. And when you have the Jacobian, estimating the errors are trivial. std(q) = sqrt ( (dq/dx std(x))*2 + (dq/dz std(z))*2 ) where q is the function. x and z are R1 and R1rho_prime. So, until then, implementing the Jacobian is only for testing the error estimation compared to Monte-Carlo simulations. Best Troels 2014-09-01 12:11 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:
Hi, I forgot to mention, for these symbol derivations, you should apply trigonometric simplifications. I don't know how sympy does this, but I always use these simplification functions in wxMaxima (or FullSimplify[] in Mathematica). It will significantly simplify the equations and make them faster to calculate. And when implementing these, you should always check against a numeric example. Here Numdifftools is very useful. You can pass dfunc_DPL94 into Numdifftools and then hardcode the numeric gradient and Hessian into a unit test. The Jacobian requires a different input function. Regards, Edward On 1 September 2014 12:07, Edward d'Auvergne <edward@xxxxxxxxxxxxx> wrote:Hi, For the DPL model, the symbolic gradient, Hessian, and Jacobian matrices are very basic. These could be added to lib.dispersion.dpl94. And then maybe as the target function class methods dfunc_DPL94(), d2func_DPL94(), and jacobian_DPL94(). For a permanent reference, it would be great to have these added to the relax manual. Maybe as a new section in Chapter 14, "Optimisation of relaxation data - values, gradients, and Hessians". The symbolic derivatives for all other analytic models should also not be too complicated, I hope. Anyway, if you are interested in having this functionality we can incorporate Numdifftools into relax to provide slow numerical estimates for the other dispersion models (https://code.google.com/p/numdifftools/). One day I might incorporate this into minfx as well, so that minfx can use numerical gradients and Hessians automatically for optimisation when the user does not provide them themselves. Cheers, Edward On 1 September 2014 11:49, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:Hi Edward. I think you dont randomize for R1. This should be a bug. Ugh. Do you submit this? If R1 is fitted, then one can just take the fitted values. I have just made the Jacobian and Hessian for DPL94. wiki.nmr-relax.com/DPL94_derivatives When the Jacobians are defined like this, the only thing necessary is: ------------------------------------------------- def func_chi2_grad(self, params=None, times=None, values=None, errors=None): """Target function for the gradient (Jacobian matrix) to minfx, for exponential fit . @param params: The vector of parameter values. @type params: numpy rank-1 float array @keyword times: The time points. @type times: numpy array @param values: The measured values. @type values: numpy array @param errors: The standard deviation of the measured intensity values per time point. @type errors: numpy array @return: The Jacobian matrix with 'm' rows of function derivatives per 'n' columns of parameters, which have been summed together. @rtype: numpy array """ # Get the back calc. back_calc = self.func(params=params) # Get the Jacobian, with partial derivative, with respect to r2eff and i0. grad = self.func_grad(params=params) # Transpose back, to get rows. grad_t = transpose(grad) # n is number of fitted parameters. n = len(params) # Define array to update parameters in. jacobian_chi2_minfx = zeros([n]) # Update value elements. dchi2(dchi2=jacobian_chi2_minfx, M=n, data=values, back_calc_vals=back_calc, back_calc_grad=grad_t, errors=errors) # Return Jacobian matrix. return jacobian_chi2_minfx ---------------------------------------------------- 2014-09-01 10:51 GMT+02:00 Edward d'Auvergne <edward@xxxxxxxxxxxxx>:Hi, Please see below: On 1 September 2014 10:20, Troels Emtekær Linnet <tlinnet@xxxxxxxxxxxxx> wrote:Yes. That was a seriously hard bug to find. Especially when you consider the MC simulations as the "Golden Standard". And then the "Golden Standard" is wrong... Ouch.True! But there are bugs everywhere and you should never assume that parts of relax, or any software in general is bug free. Never trust black-boxes! This is a good lesson ;) All software has bugs, and this will not be the last in relax.Should we set the GUI to have exp_sim = -1? There is no assumption, that 100000 simulations of exponential fits are better than the co-variance.Just in case someone has much harder to optimise peak intensity data, for example if the data is extremely noisy or if it is not exactly mono-exponential, then this is not a good idea. It is better to spend time to obtain the best result rather than obtaining a quick result which in most cases works, but is known to theoretically fail. You don't want to be in that theoretical failure group and not know about it. So the user can set it themselves but they should compare it to MC simulations anyway to be sure. Note that the curvature of the optimisation space for the dispersion models is far more complicated than the 2-parameter exponential curves. For the dispersion models, the covariance matrix approach will not be anywhere near as good. For these models, most big names in the field would never consider the covariance matrix approach. Many people are wary of the edge-case failures of this technique. For the best results you would always pick the best technique which, by statistical theory, is Monte Carlo simulations by far.Btw. Can we check Monte-Carlo simulations for the dispersion models?That's a great idea! This is probably also untested in the test suite. The covariance matrix approach is perfect for checking that the Monte Carlo results are reasonable. However you do require the Jacobian matrix which is not derived for any dispersion model. There are no gradients derived, though it could be done numerically in the test_suite/shared_data/dispersion directories using the very useful Numdifftools package (https://pypi.python.org/pypi/Numdifftools). Or an even better way would be to create the error_analysis.covariance_matrix user function which, like the pipe_control.monte_carlo module, uses the specific API to obtain, in this case the Jacobian and weighting matrix via one new method, calls lib.statistics.multifit_covar() to create the covariance matrix, and then calls the API again to set the errors via the already existing api.set_error() API method. Then you can use the covariance matrix approach for all the dispersion models. Due to the licencing of Numdifftools, we could even bundle that with relax in the extern package and use numerical Jacobian integration so that even the numeric dispersion models can have a covariance matrix.Where is that performed?The specific analysis API. See the functions in the pipe_control.monte_carlo module. The API object is obtained as: # The specific analysis API object. api = return_api() Then you can see the methods called in specific_analyses.relax_disp.api as, for example in the pipe_control.monte_carlo module: # Create the Monte Carlo data. if method == 'back_calc': data = api.create_mc_data(data_index) You will therefore find the create_mc_data() method in the dispersion API module. If you search for all of the api.*() calls, then you'll find all those methods in the API object (or the default with a RelaxImplementError in the API base class). It's rather simple.Do you randomize only R1rho' or do you also include randomize for R1?This is performed in the pipe_control.monte_carlo.create_data() function. See if you can trace the API method calls back and find the source of this! It would be good if you check as well. Cheers, Edward