leastsq_bounds, a wrapper for scipy.optimize.leastsq with box bounds

Leastsq_bounds.py

# leastsq_bounds.py
# see also test_leastsq_bounds.py on gist.github.com/denis-bz

from __future__ import division
import numpy as np
from scipy.optimize import leastsq

__version__ = "2015-01-10 jan  denis"  # orig 2012


#...............................................................................
def leastsq_bounds( func, x0, bounds, boundsweight=10, **kwargs ):
    """ leastsq with bound conatraints lo <= p <= hi
    run leastsq with additional constraints to minimize the sum of squares of
        [func(p) ...]
        + boundsweight * [max( lo_i - p_i, 0, p_i - hi_i ) ...]

    Parameters
    ----------
    func() : a list of function of parameters `p`, [err0 err1 ...]
    bounds : an n x 2 list or array `[[lo_0,hi_0], [lo_1, hi_1] ...]`.
        Use e.g. [0, inf]; do not use NaNs.
        A bound e.g. [2,2] pins that x_j == 2.
    boundsweight : weights the bounds constraints
    kwargs : keyword args passed on to leastsq

    Returns
    -------
    exactly as for leastsq,
http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html

    Notes
    -----
    The bounds may not be met if boundsweight is too small;
    check that with e.g. check_bounds( p, bounds ) below.

    To access `x` in `func(p)`, `def func( p, x=xouter )`
    or make it global, or `self.x` in a class.

    There are quite a few methods for box constraints;
    you'll maybe sing a longer song ...
    Comments are welcome, test cases most welcome.

"""
    # Example: test_leastsq_bounds.py

    if bounds is not None  and  boundsweight > 0:
        check_bounds( x0, bounds )
        if "args" in kwargs:  # 8jan 2015
            args = kwargs["args"]
            del kwargs["args"]
        else:
            args = ()
#...............................................................................
        funcbox = lambda p: \
            np.hstack(( func( p, *args ),
                        _inbox( p, bounds, boundsweight ))) 
    else:
        funcbox = func
    return leastsq( funcbox, x0, **kwargs )


def _inbox( X, box, weight=1 ):
    """ -> [tub( Xj, loj, hij ) ... ]
        all 0  <=>  X in box, lo <= X <= hi
    """
    assert len(X) == len(box), \
        "len X %d != len box %d" % (len(X), len(box))
    return weight * np.array([
        np.fmax( lo - x, 0 ) + np.fmax( 0, x - hi )
            for x, (lo,hi) in zip( X, box )])

# def tub( x, lo, hi ):
#     """ \___/  down to lo, 0 lo .. hi, up from hi """
#     return np.fmax( lo - x, 0 ) + np.fmax( 0, x - hi )

#...............................................................................
def check_bounds( X, box ):
    """ print Xj not in box, loj <= Xj <= hij
        return nr not in
    """
    nX, nbox = len(X), len(box)
    assert nX == nbox, \
        "len X %d != len box %d" % (nX, nbox)
    nnotin = 0
    for j, x, (lo,hi) in zip( range(nX), X, box ):
        if not (lo <= x <= hi):
            print "check_bounds: x[%d] %g is not in box %g .. %g" % (j, x, lo, hi)
            nnotin += 1
    return nnotin

test_leastsq_bounds.py

# test_leastsq_bounds.py 10jan
# from http://stackoverflow.com/questions/6949370/scipy-leastsq-dfun-usage

from __future__ import division
import sys
import numpy as np
from scipy.optimize import leastsq  # $scipy/optimize/minpack.py

import leastsq_bounds as lsqb

#...............................................................................
pmin = np.r_[3,5,1]
pstart = np.r_[10,10,10]
box = np.c_[ pmin + 1, pmin + 10 ]
wt = 10
nx = 10
plot = 0
seed = 0

    # to change these params in sh or ipython, run this.py  a=1  b=None  c=[3] ...
for arg in sys.argv[1:]:
    exec( arg )
np.set_printoptions( 1, threshold=100, edgeitems=10, linewidth=100,
    formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g
np.random.seed(seed)

if plot:
    from matplotlib import pyplot as pl
    try:
        import seaborn  # http://stanford.edu/~mwaskom/software/seaborn
    except ImportError:
        pass
    pl.rcParams['figure.figsize'] = 10, 5
    pl.subplots_adjust( top=.9 )

#...............................................................................
def func(p, x):
    return p[0]*np.exp(-p[1]*x) + p[2]

def dfunc(p, x, y):     #Derivative
    return [np.exp(-p[1]*x), - x*p[0]*np.exp(-p[1]*x), np.ones(len(x))]

def residuals(p, x, y):
    return func(p, x) - y  # care

#Generate messy data
x = np.linspace( 0, 5, 5*nx | 1 )
y0 = func( pmin, x )
y = y0 + np.random.normal(size=len(y0))

    # lsq = leastsq( residuals, pstart, args=(x, y), full_output=True )
            # Dfun=dfunc, col_deriv=1,
    # p, cov, infodict, mesg, ier = lsq
    # print "leastsq p:", p  # pmin

print "pmin:", pmin
print "box:\n", box

#...............................................................................
lsq = lsqb.leastsq_bounds( residuals, pstart, bounds=box, boundsweight=wt,
        args=(x, y), full_output=True )
p, cov, infodict, mesg, ier = lsq

print "leastsq_bounds weight %.2g -> %s" % (wt, p)
if not (1 <= ier <= 4):
    print mesg
funcp = func( p, x )
res = funcp - y
print "residuals func( p, x ) - y:  av %.2g  |av| %.2g \n%s" % (
        res.mean(), np.fabs(res).mean(), res )

if plot:
    pl.title( "leastsq_bounds fit p0 e^(-p1 x) + p2\n"
            "in box outside pmin %s  -> p %s" % (pmin, p) )
    pl.ylim( -1, 5 )
    pl.plot( x, funcp, label="func( p, x )" )
    pl.plot( x, y,     label="y" )
    pl.legend()
    pl.show()