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
