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 <tlin...@nmr-relax.com> 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 <edw...@nmr-relax.com>: >> Hi, >> >> Please see below: >> >> On 1 September 2014 10:20, Troels Emtekær Linnet <tlin...@nmr-relax.com> >> 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 _______________________________________________ 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
