[Scipy-tickets] [SciPy] #1588: scipy.optimize.cobyla not consistant in Windows

SciPy Trac scipy-tickets@scipy....
Mon Jan 23 14:11:41 CST 2012

#1588: scipy.optimize.cobyla not consistant in Windows
 Reporter:  casperskovby  |       Owner:  somebody   
     Type:  defect        |      Status:  new        
 Priority:  normal        |   Milestone:  Unscheduled
Component:  Other         |     Version:  0.10.0     
 Keywords:                |  

 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 XP. 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 in Windows XP.
 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 had 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

 He suggested me to contact this group. Have you any ideas where the bug
 can be located?

 Kind Regards,

 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,

         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')

Ticket URL: <http://projects.scipy.org/scipy/ticket/1588>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.

More information about the Scipy-tickets mailing list