[SciPy-Dev] scipy.optimize.cobyla not consistant in Windows
Casper Skovby
casperskovby@gmail....
Thu Jan 12 15:20:06 CST 2012
Hello again.
I have not got any response on this question yet. Does anyone have a clue
about this?
On Tue, Jan 3, 2012 at 11:05 PM, Casper Skovby <casperskovby@gmail.com>wrote:
> Hello.
>
>
> I have realized that scipy.optimize.cobyla.fmin_cobyla does not always
> give the same result even though the input is exactly the same when running
> in Windows. When running in Linux I do not experience this. Below I have
> added an example. In this case as you can see in the example I am trying
> to represent a curve with a cubic spline. As input to the spline I have 12
> points – two end points and 10 inner points. I am trying to optimize these
> 10 inner points (both in x and y direction) in order to represent the curve
> as good as possible with this spline. But as said it does not always ends
> up with the same result. Sometimes you can run the code several times with
> the same result but suddenly the result differ.
>
> I realized this bug when I used OpenOpt (openopt.org) because it gave
> similar problems. I have a conversation with the Developer of OpenOpt
> Dmitrey (see http://forum.openopt.org/viewtopic.php?id=499). When using
> OpenOpt in Linux some of the solvers do also give inconsistant results.
> Dmitreys conclusion is: Since scipy_cobyla works different in Windows and
> Linux, probably something with f2py is wrong (or something in libraries it
> involves).
>
> He suggested me to contact this group. Have you any ideas where the bug
> can be located?
>
> Kind Regards,
> Casper
>
> from matplotlib.pylab import *
> import numpy as np
> from numpy import linspace
> from scipy import interpolate
> from scipy.optimize import cobyla
>
> def residual(p, r, y):
> tck = interpolate.splrep(np.concatenate(([r[0]], p[0:len(p)/2], [r[-1]])),
> np.concatenate(([y[0]], p[len(p)/2::], [y[-1]])))
> yFit = interpolate.splev(r, tck)
>
> resid = sum((y-yFit)**2)
> return resid
>
> def constraint(x):
> return 1
>
>
> def FitDistribution(r, y):
> rP = linspace(r[0], r[-1], 12)
> yP = np.interp(rP, r, y)
> p0 = np.concatenate((rP[1:-1], yP[1:-1]))
> # form box-bound constraints lb <= x <= ub
> lb = np.concatenate((np.zeros(len(rP)-2), np.ones(len(rP)-2)*(-np.inf))) #
> lower bound
> ub = np.concatenate((np.ones(len(rP)-2)*r[-1],
> np.ones(len(rP)-2)*(np.inf))) # upper bound
> ftol = 10e-5
> f = lambda p: residual(p, r, y)
>
> pvec = cobyla.fmin_cobyla(f, p0, constraint, iprint=1, maxfun=100000)
> #def objective(x):
> # return x[0]*x[1]
> #
> #def constr1(x):
> # return 1 - (x[0]**2 + x[1]**2)
> #
> #def constr2(x):
> # return x[1]
> #
> #x=cobyla.fmin_cobyla(objective, [0.0, 0.1], [constr1, constr2],
> rhoend=1e-7, iprint=0)
>
>
> tck = interpolate.splrep(np.concatenate(([r[0]], pvec[0:len(pvec)/2],
> [r[-1]])), np.concatenate(([y[0]], pvec[len(pvec)/2::], [y[-1]])))
> yFit = interpolate.splev(r, tck)
> return yFit
>
> r = linspace(1,100)
> y = r**2
>
> yFit = FitDistribution(r, y)
>
> plot(r, y, label = 'func')
> plot(r, yFit, label = 'fit')
> legend(loc=0)
> grid()
> show()
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-dev/attachments/20120112/0af99486/attachment.html
More information about the SciPy-Dev
mailing list