It seems you form: self.xk = x0 self.I = identity(len(self.xk))
p = self.I So self.I will be 2 x 2 matrix. And a is 4x4. Best Troels 2014-08-25 19:55 GMT+02:00 Troels Emtekær Linnet <[email protected]>: > Hi Ed. > > I get: > results = newton(func=func, dfunc=dfunc, d2func=d2func, args=args, > x0=x0, min_options=min_options, func_tol=func_tol, grad_tol=grad_tol, > maxiter=maxiter, full_output=full_output, print_flag=print_flag, > print_prefix=print_prefix) > File > "/sbinlab2/software/python-enthought-dis/canopy-1.4.0-full-rh5-64/Canopy_64bit/User/lib/python2.7/site-packages/minfx/newton.py", > line 51, in newton > results = min.minimise() > File > "/sbinlab2/software/python-enthought-dis/canopy-1.4.0-full-rh5-64/Canopy_64bit/User/lib/python2.7/site-packages/minfx/base_classes.py", > line 262, in minimise > self.new_param_func() > File > "/sbinlab2/software/python-enthought-dis/canopy-1.4.0-full-rh5-64/Canopy_64bit/User/lib/python2.7/site-packages/minfx/newton.py", > line 169, in new_param_func > self.pk = self.get_pk() > File > "/sbinlab2/software/python-enthought-dis/canopy-1.4.0-full-rh5-64/Canopy_64bit/User/lib/python2.7/site-packages/minfx/base_classes.py", > line 668, in gmw > return gmw(self.dfk, self.d2fk, self.I, self.n, self.mach_acc, > self.print_prefix, self.print_flag, return_matrix) > File > "/sbinlab2/software/python-enthought-dis/canopy-1.4.0-full-rh5-64/Canopy_64bit/User/lib/python2.7/site-packages/minfx/hessian_mods/gmw81.py", > line 101, in gmw > a = dot(p, dot(a, p)) > ValueError: matrices are not aligned > > > print p.shape = (2, 2) > print a.shape = (4, 4) > > Do you know what is wrong? > > Best > Troels > > 2014-08-25 19:15 GMT+02:00 Edward d'Auvergne <[email protected]>: >> Hi, >> >> Maybe to help work out what is happing both in Scipy and in relax, a >> simple test data set could be created. For example if you set I0 to >> 1000, R to 1, and then pick the time points 0, 1, 2, 3, 4, you would >> then have the intensity data: >> >> 0 1000.0 >> 1 367.879441171 >> 2 135.335283237 >> 3 49.7870683679 >> 4 18.3156388887 >> >> It might take some Bootstrapping simulations to obtain intensity >> errors. Just randomise I0 and R with a Gaussian distribution a large >> number of times, then work out the intensity SD values for each time >> point. >> >> Then you could pass the data and errors into the Scipy and relax code >> paths and see what you get for a covariance matrix. That might show >> if the double R2eff error is a bug or a failure of the covariance >> matrix error estimate. As it is synthetic data, you will know the >> exact parameter values. And you will know the exact error values to >> search for in the covariance matrix. >> >> Regards, >> >> Edward >> >> >> >> On 25 August 2014 18:02, Edward d'Auvergne <[email protected]> wrote: >>> Hi, >>> >>> I just found a useful description of Jacobian construction with an example: >>> >>> http://www.mathworks.de/de/help/optim/ug/writing-objective-functions.html#brkjt9_ >>> >>> In our case we have a vector function. Given a time, the peak >>> intensity will be returned as >>> >>> I = I0.exp(-R.t). >>> >>> Regards, >>> >>> Edward >>> >>> >>> On 25 August 2014 17:51, Edward d'Auvergne <[email protected]> wrote: >>>> They should be the same ordering as the parameter vector, as each >>>> column is a partial derivative with respect to that parameter. The >>>> rows should be the x data or the time points. Well, I hope this is >>>> correct, it's been many years since I coded the model-free Jacobian >>>> matrix. Rather than using the new C code, maybe a test data set >>>> should be calculated by hand to determine if the covariance matrix >>>> calculation, and Jacobian calculation, is correct. >>>> >>>> Regards, >>>> >>>> Edward >>>> >>>> >>>> >>>> On 25 August 2014 17:47, Troels Emtekær Linnet <[email protected]> >>>> wrote: >>>>> Is the order of the columns in the Jacobian matrix important? >>>>> >>>>> Best >>>>> Troels >>>>> >>>>> 2014-08-25 17:42 GMT+02:00 Edward d'Auvergne <[email protected]>: >>>>>> That looks correct. If you calculate: >>>>>> >>>>>> linalg.inv(dot(transpose(mat), mat)) >>>>>> >>>>>> Do you get the covariance matrix? >>>>>> >>>>>> Regards, >>>>>> >>>>>> Edward >>>>>> >>>>>> >>>>>> >>>>>> On 25 August 2014 17:35, Troels Emtekær Linnet <[email protected]> >>>>>> wrote: >>>>>>> Let me exemplify: >>>>>>> >>>>>>> There is 4 time points. >>>>>>> >>>>>>> df_d_i0 = - ( 2. * ( self.values - i0 / exp(r2eff * self.relax_times) >>>>>>> ) ) / ( exp(r2eff * self.relax_times) * self.errors**2) >>>>>>> df_d_r2eff = (2 * i0 * self.relax_times * (self.values - i0 / >>>>>>> exp(r2eff * self.relax_times) ) ) / ( exp(r2eff * self.relax_times) * >>>>>>> self.errors**2) >>>>>>> >>>>>>> Should the return then be: >>>>>>> >>>>>>> print df_d_i0.shape, df_d_i0 >>>>>>> (4,) [-0.004826723918314 -0.00033019968656 0.002366308749814 >>>>>>> -0.000232558176186] >>>>>>> >>>>>>> print df_d_r2eff.shape, df_d_r2eff >>>>>>> (4,) [ 0. 2.66126225080615 -47.678483702132965 >>>>>>> 9.371576058231405] >>>>>>> >>>>>>> mat = transpose(array( [df_d_i0, df_d_r2eff ]) ) >>>>>>> print mat.shape, mat >>>>>>> >>>>>>> >>>>>>> (4, 2) [[ -4.826723918313830e-03 0.000000000000000e+00] >>>>>>> [ -3.301996865596296e-04 2.661262250806150e+00] >>>>>>> [ 2.366308749814298e-03 -4.767848370213297e+01] >>>>>>> [ -2.325581761857821e-04 9.371576058231405e+00]] >>>>>>> >>>>>>> >>>>>>> Best >>>>>>> Troels >>>>>>> >>>>>>> 2014-08-25 17:27 GMT+02:00 Edward d'Auvergne <[email protected]>: >>>>>>>> Hi, >>>>>>>> >>>>>>>> I may have explained this incorrectly earlier: >>>>>>>> >>>>>>>> - The vector of partial derivatives with respect to the chi-squared >>>>>>>> equation is the gradient. >>>>>>>> - The vector of partial derivatives with respect to the exponential >>>>>>>> function is the Jacobian. >>>>>>>> >>>>>>>> The equations and code for the exponential partial derivatives are the >>>>>>>> same in both. It's just that they are used differently. Does this >>>>>>>> help? >>>>>>>> >>>>>>>> Regards, >>>>>>>> >>>>>>>> Edward >>>>>>>> >>>>>>>> On 25 August 2014 17:26, Troels Emtekær Linnet <[email protected]> >>>>>>>> wrote: >>>>>>>>> Hi Edward. >>>>>>>>> >>>>>>>>> When writing the Jacobian, do you then derivative according to ( i0 * >>>>>>>>> exp( -times * r2eff) ) >>>>>>>>> or do you do the derivative according to chi2 function? >>>>>>>>> >>>>>>>>> I have a little hard time to figure out the code text. >>>>>>>>> >>>>>>>>> From minfx: >>>>>>>>> @keyword func: The function which returns the value. >>>>>>>>> @type func: func >>>>>>>>> >>>>>>>>> So, this is the chi2 function. >>>>>>>>> >>>>>>>>> @keyword dfunc: The function which returns the gradient. >>>>>>>>> @type dfunc: func >>>>>>>>> >>>>>>>>> So, this must be the derivative of the chi2 function? >>>>>>>>> >>>>>>>>> So in essence. >>>>>>>>> >>>>>>>>> Does minfx expect a "dfunc" function which calculate the: >>>>>>>>> >>>>>>>>> one gradient chi2 value, subject to the input parameters? >>>>>>>>> >>>>>>>>> Or >>>>>>>>> A jacobian matrix of the form: >>>>>>>>> m X n matrix, where m is the number of time elements and n is number >>>>>>>>> of parameters = 2. >>>>>>>>> >>>>>>>>> Best >>>>>>>>> Troels >>>>>>>>> >>>>>>>>> 2014-08-25 15:52 GMT+02:00 Edward d'Auvergne <[email protected]>: >>>>>>>>>> Hi Troels, >>>>>>>>>> >>>>>>>>>> Please see below: >>>>>>>>>> >>>>>>>>>> On 25 August 2014 13:01, Troels Emtekær Linnet >>>>>>>>>> <[email protected]> wrote: >>>>>>>>>>> 2014-08-25 11:19 GMT+02:00 Edward d'Auvergne <[email protected]>: >>>>>>>>>>>> Hi Troels, >>>>>>>>>>>> >>>>>>>>>>>> Unfortunately you have gone ahead an implemented a solution without >>>>>>>>>>>> first discussing or planning it. Hence the current solution has a >>>>>>>>>>>> number of issues: >>>>>>>>>>>> >>>>>>>>>>>> 1) Target function replication. The solution should have reused >>>>>>>>>>>> the >>>>>>>>>>>> C modules. The original Python code for fitting exponential curves >>>>>>>>>>>> was converted to C code for speed >>>>>>>>>>>> (http://gna.org/forum/forum.php?forum_id=1043). Note that two >>>>>>>>>>>> point >>>>>>>>>>>> exponentials that decay to zero is not the only way that data can >>>>>>>>>>>> be >>>>>>>>>>>> collected, and that is the reason for Sebastien Morin's >>>>>>>>>>>> inversion-recovery branch (which was never completed). Anyway, the >>>>>>>>>>>> code duplication is not acceptable. If the C module is extended >>>>>>>>>>>> with >>>>>>>>>>>> new features, such as having the true gradient and Hessian >>>>>>>>>>>> functions, >>>>>>>>>>>> then the Python module will then be out of sync. And vice-versa. >>>>>>>>>>>> If >>>>>>>>>>>> a bug is found in one module and fixed, it may still be present in >>>>>>>>>>>> the >>>>>>>>>>>> second. This is a very non-ideal situation for relax to be in, >>>>>>>>>>>> and is >>>>>>>>>>>> the exact reason why I did not allow the cst branch to be merged >>>>>>>>>>>> back >>>>>>>>>>>> to trunk. >>>>>>>>>>> >>>>>>>>>>> Hi Edward. >>>>>>>>>>> >>>>>>>>>>> I prefer not to make this target function dependent on C-code >>>>>>>>>>> compilation. >>>>>>>>>>> >>>>>>>>>>> Compilation of code on windows can be quite a hairy thing. >>>>>>>>>>> >>>>>>>>>>> For example see: >>>>>>>>>>> http://wiki.nmr-relax.com/Installation_windows_Python_x86-32_Visual_Studio_Express_for_Windows_Desktop#Install_Visual_Studio_Express_2012_for_Windows_Desktop >>>>>>>>>>> >>>>>>>>>>> Visua Studio Express is several hundreds of megabyte installation, >>>>>>>>>>> for >>>>>>>>>>> just compiling an exponential curve. ? >>>>>>>>>>> This is way, way overkill for this situation. >>>>>>>>>> >>>>>>>>>> The C code compilation has been a requirement in relax since 2006. >>>>>>>>>> This was added not only for speed, but as a framework to copy for >>>>>>>>>> other analysis types in the future. Once a Python target function >>>>>>>>>> has >>>>>>>>>> been fully optimised, for the last speed up the code can be converted >>>>>>>>>> to C. This is the future plan for a number of the relax analyses. >>>>>>>>>> But first the Python code is used for prototyping and for finding the >>>>>>>>>> fastest implementation/algorithm. >>>>>>>>>> >>>>>>>>>> The C compilation will become an even greater requirement once I >>>>>>>>>> write >>>>>>>>>> C wrapper code for QUADPACK to eliminate the last dependencies on >>>>>>>>>> Scipy. And the C compilation framework allows for external C and >>>>>>>>>> FORTRAN libraries to be added to the 'extern' package in the future, >>>>>>>>>> as there are plenty of open source libraries out there with >>>>>>>>>> compatible >>>>>>>>>> licences which could be very useful to use within relax. >>>>>>>>>> >>>>>>>>>> >>>>>>>>>>>> 2) Scipy is now a dependency for the dispersion analysis! Why was >>>>>>>>>>>> this not discussed? Coding a function for calculating the >>>>>>>>>>>> covariance >>>>>>>>>>>> matrix is basic. Deriving and coding the real gradient function is >>>>>>>>>>>> also basic. I do not understand why Scipy is now a dependency. I >>>>>>>>>>>> have been actively trying to remove Scipy as a relax dependency and >>>>>>>>>>>> only had a single call for numeric quadratic intergration via >>>>>>>>>>>> QUADPACK >>>>>>>>>>>> wrappers left to remove for the frame order analysis. Now Scipy is >>>>>>>>>>>> back :( >>>>>>>>>>> >>>>>>>>>>> Hi Edward. >>>>>>>>>>> >>>>>>>>>>> Scipy is a dependency for trying calculation with >>>>>>>>>>> scipy.optimize.leastsq. >>>>>>>>>>> >>>>>>>>>>> How could it be anymore different? >>>>>>>>>>> >>>>>>>>>>> What you are aiming at, is to add yet another feature for >>>>>>>>>>> estimating the errors. >>>>>>>>>>> A third solution. >>>>>>>>>>> >>>>>>>>>>> What ever the third solution would come up with of dependency, would >>>>>>>>>>> depend on the method implemented. >>>>>>>>>>> One could also possible imagine to extend this procedures in R, >>>>>>>>>>> Matlab >>>>>>>>>>> or whatever. >>>>>>>>>>> >>>>>>>>>>> Byt they would also need to meet some dependencies. >>>>>>>>>>> >>>>>>>>>>> Of course the best solution would always try to make relax most >>>>>>>>>>> independent. >>>>>>>>>>> >>>>>>>>>>> But if the desire is to try with scipy.optimize.leastsq, then you >>>>>>>>>>> are >>>>>>>>>>> bound with this dependency. >>>>>>>>>> >>>>>>>>>> That's why I asked if only the covariance matrix is required. Then >>>>>>>>>> we >>>>>>>>>> can replace the use of scipy.optimize.leastsq() with a single >>>>>>>>>> function >>>>>>>>>> for calculating the covariance matrix. >>>>>>>>>> >>>>>>>>>> >>>>>>>>>>>> 3) If the covariance function was coded, then the specific >>>>>>>>>>>> analysis >>>>>>>>>>>> API could be extended with a new covariance method and the >>>>>>>>>>>> relax_disp.r2eff_estimate user function could have simply been >>>>>>>>>>>> called >>>>>>>>>>>> error_estimate.covariance_matrix, or something like that. Then >>>>>>>>>>>> this >>>>>>>>>>>> new error_estimate.covariance_matrix user function could replace >>>>>>>>>>>> the >>>>>>>>>>>> monte_carlo user functions for all analyses, as a rough error >>>>>>>>>>>> estimator. >>>>>>>>>>> >>>>>>>>>>> That would be the third possibility. >>>>>>>>>> >>>>>>>>>> ..., that would give the same result, save the same amount of time, >>>>>>>>>> but would avoid the new Scipy dependency and be compatible with all >>>>>>>>>> analysis types ;) >>>>>>>>>> >>>>>>>>>> >>>>>>>>>>>> 4) For the speed of optimisation part of the new >>>>>>>>>>>> relax_disp.r2eff_estimate user function, this is not because scipy >>>>>>>>>>>> is >>>>>>>>>>>> faster than minfx!!! It is the choice of algorithms, the numerical >>>>>>>>>>>> gradient estimate, etc. >>>>>>>>>>>> (http://thread.gmane.org/gmane.science.nmr.relax.scm/22979/focus=6812). >>>>>>>>>>> >>>>>>>>>>> This sound good. >>>>>>>>>>> >>>>>>>>>>> But I can only say, that as I user I meet a "big wall of time >>>>>>>>>>> consumption", for the error >>>>>>>>>>> estimation of R2eff via Monte-Carlo. >>>>>>>>>>> >>>>>>>>>>> As a user, I needed more options to try out. >>>>>>>>>> >>>>>>>>>> The idea of adding the covariance matrix error estimate to relax is a >>>>>>>>>> great idea. Despite its lower quality, it is hugely faster than >>>>>>>>>> Monte >>>>>>>>>> Carlo simulations. It has been considered it before, see >>>>>>>>>> http://thread.gmane.org/gmane.science.nmr.relax.user/602/focus=629 >>>>>>>>>> and >>>>>>>>>> the discussions in that thread. But the time required for Monte >>>>>>>>>> Carlo >>>>>>>>>> simulations was never an issue so the higher quality estimate >>>>>>>>>> remained >>>>>>>>>> the only implementation. >>>>>>>>>> >>>>>>>>>> What I'm trying to do, is to direct your solution to be general and >>>>>>>>>> reusable. I'm also thinking of other techniques at the same time, >>>>>>>>>> Jackknife simulations for example, which could be added in the future >>>>>>>>>> by developers with completely different interests. >>>>>>>>>> >>>>>>>>>> >>>>>>>>>>>> 5) Back to Scipy. Scipy optimisation is buggy full stop. The >>>>>>>>>>>> developers ignored my feedback back in 2003. I assumed that the >>>>>>>>>>>> original developers had just permanently disappeared, and they >>>>>>>>>>>> really >>>>>>>>>>>> never came back. The Scipy optimisation code did not change for >>>>>>>>>>>> many, >>>>>>>>>>>> many years. While it looks like optimisation works, in some cases >>>>>>>>>>>> it >>>>>>>>>>>> does fails hard, stopping in a position in the space where there >>>>>>>>>>>> is no >>>>>>>>>>>> minimum! I added the dx.map user function to relax to understand >>>>>>>>>>>> these Scipy rubbish results. And I created minfx to work around >>>>>>>>>>>> these >>>>>>>>>>>> nasty hidden failures. I guess such failures are due to them not >>>>>>>>>>>> testing the functions as part of a test suite. Maybe they have >>>>>>>>>>>> fixed >>>>>>>>>>>> the bugs now, but I really can no longer trust Scipy optimisation. >>>>>>>>>>>> >>>>>>>>>>> >>>>>>>>>>> I am sorry to hear about this. >>>>>>>>>>> >>>>>>>>>>> And I am totally convinced that minfx is better for minimising the >>>>>>>>>>> dispersion models. >>>>>>>>>>> You have proven that quite well in your papers. >>>>>>>>>>> >>>>>>>>>>> I do though have a hard time believing that minimisation of an >>>>>>>>>>> exponential function should be >>>>>>>>>>> subject to erroneous results. >>>>>>>>>>> >>>>>>>>>>> Anyway, this is still left to "freedom of choice" for the user. >>>>>>>>>> >>>>>>>>>> The error in the original Scipy optimisation code was causing quite >>>>>>>>>> different results. The 3 algorithms, now that I look back at my >>>>>>>>>> emails from 2003, are: >>>>>>>>>> >>>>>>>>>> - Nelder-Mead simplex, >>>>>>>>>> - Levenberg-Marquardt, >>>>>>>>>> - NCG. >>>>>>>>>> >>>>>>>>>> These are still all present in Scipy, though I don't know if the code >>>>>>>>>> is different from back in 2003. The error in the Levenberg-Marquardt >>>>>>>>>> algorithm was similar to the Modelfree4 problem, in that a lamba >>>>>>>>>> matrix updating condition was incorrectly checked for. When the >>>>>>>>>> gradient was positive, i.e. up hill, the matrix should update and the >>>>>>>>>> algorithm continue to try to find a downhill step. If the conditions >>>>>>>>>> are not correctly checked for, the algorithm thinks that the up hill >>>>>>>>>> step means that it is at the minimum. But this is not the case, it >>>>>>>>>> is >>>>>>>>>> just pointing in the wrong direction. I don't remember what the NCG >>>>>>>>>> bug was, but that one was much more severe and the results were >>>>>>>>>> strange. >>>>>>>>>> >>>>>>>>>> Failures of optimisation algorithms due to bugs can be quite random. >>>>>>>>>> And you often don't see them, as you don't know what the true result >>>>>>>>>> really is. But such bugs will affect exponential functions, despite >>>>>>>>>> their simplicity. >>>>>>>>>> >>>>>>>>>> Regards, >>>>>>>>>> >>>>>>>>>> Edward _______________________________________________ relax (http://www.nmr-relax.com) This is the relax-devel mailing list [email protected] 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
