Context Navigation

← Previous Changeset
Next Changeset →

Changeset 418

Timestamp:

08/25/10 10:58:37 (6 years ago)

Author:

mmckerns

Message:

centralized support for (modified) scipy helper functions

Location:

branches/alta/mystic-0.2a1

Files:

: 1 added
: 6 edited

Make.mm (modified) (1 diff)
_scipy060optimize.py (modified) (1 diff)
_scipyoptimize.py (added)
scipy_bfgs.py (modified) (1 diff)
scipy_cg.py (modified) (1 diff)
scipy_ncg.py (modified) (1 diff)
scipy_optimize.py (modified) (2 diffs)

Legend:

: Unmodified
: Added
: Removed

branches/alta/mystic-0.2a1/Make.mm

r238	r418
28	28	_genSow.py \
29	29	_scipy060optimize.py \
	30	_scipyoptimize.py \
30	31	abstract_map_solver.py \
31	32	abstract_nested_solver.py \

branches/alta/mystic-0.2a1/_scipy060optimize.py

-                      r124
+                      r418
 """local copy of scipy.optimize"""
 __all__ = ['fmin', 'fmin_powell', 'fmin_ncg',
            'fminbound','brent', 'golden','bracket','rosen','rosen_der',
            'rosen_hess', 'rosen_hess_prod', 'brute', 'approx_fprime',
            'line_search', 'check_grad']
+#__all__ = ['fmin', 'fmin_powell', 'fmin_ncg',
+#           'fminbound','brent', 'golden','bracket','rosen','rosen_der',
+#           'rosen_hess', 'rosen_hess_prod', 'brute', 'approx_fprime',
+#           'line_search', 'check_grad']
 import numpy

branches/alta/mystic-0.2a1/scipy_bfgs.py

-                      r404
+                      r418
 # Helper methods and classes
+_epsilon = sqrt(numpy.finfo(float).eps)
+def zoom(a_lo, a_hi, phi_lo, phi_hi, derphi_lo,
+         phi, derphi, phi0, derphi0, c1, c2):
+    maxiter = 10
+    i = 0
+    delta1 = 0.2  # cubic interpolant check
+    delta2 = 0.1  # quadratic interpolant check
+    phi_rec = phi0
+    a_rec = 0
+    while 1:
+        # interpolate to find a trial step length between a_lo and a_hi
+        # Need to choose interpolation here.  Use cubic interpolation and then if the
+        #  result is within delta * dalpha or outside of the interval bounded by a_lo or a_hi
+        #  then use quadratic interpolation, if the result is still too close, then use bisection
+        dalpha = a_hi-a_lo;
+        if dalpha < 0: a,b = a_hi,a_lo
+        else: a,b = a_lo, a_hi
+        # minimizer of cubic interpolant
+        #    (uses phi_lo, derphi_lo, phi_hi, and the most recent value of phi)
+        #      if the result is too close to the end points (or out of the interval)
+        #         then use quadratic interpolation with phi_lo, derphi_lo and phi_hi
+        #      if the result is stil too close to the end points (or out of the interval)
+        #         then use bisection
+        if (i > 0):
+            cchk = delta1*dalpha
+            a_j = _cubicmin(a_lo, phi_lo, derphi_lo, a_hi, phi_hi, a_rec, phi_rec)
+        if (i==0) or (a_j is None) or (a_j > b-cchk) or (a_j < a+cchk):
+            qchk = delta2*dalpha
+            a_j = _quadmin(a_lo, phi_lo, derphi_lo, a_hi, phi_hi)
+            if (a_j is None) or (a_j > b-qchk) or (a_j < a+qchk):
+                a_j = a_lo + 0.5*dalpha
+#                print "Using bisection."
+#            else: print "Using quadratic."
+#        else: print "Using cubic."
+        # Check new value of a_j
+        phi_aj = phi(a_j)
+        if (phi_aj > phi0 + c1*a_j*derphi0) or (phi_aj >= phi_lo):
+            phi_rec = phi_hi
+            a_rec = a_hi
+            a_hi = a_j
+            phi_hi = phi_aj
+        else:
+            derphi_aj = derphi(a_j)
+            if abs(derphi_aj) <= -c2*derphi0:
+                a_star = a_j
+                val_star = phi_aj
+                valprime_star = derphi_aj
+                break
+            if derphi_aj*(a_hi - a_lo) >= 0:
+                phi_rec = phi_hi
+                a_rec = a_hi
+                a_hi = a_lo
+                phi_hi = phi_lo
+            else:
+                phi_rec = phi_lo
+                a_rec = a_lo
+            a_lo = a_j
+            phi_lo = phi_aj
+            derphi_lo = derphi_aj
+        i += 1
+        if (i > maxiter):
+            a_star = a_j
+            val_star = phi_aj
+            valprime_star = None
+            break
+    return a_star, val_star, valprime_star
+def _cubicmin(a,fa,fpa,b,fb,c,fc):
+    # finds the minimizer for a cubic polynomial that goes through the
+    #  points (a,fa), (b,fb), and (c,fc) with derivative at a of fpa.
+    #
+    # if no minimizer can be found return None
+    #
+    # f(x) = A *(x-a)^3 + B*(x-a)^2 + C*(x-a) + D
+    C = fpa
+    D = fa
+    db = b-a
+    dc = c-a
+    if (db == 0) or (dc == 0) or (b==c): return None
+    denom = (db*dc)**2 * (db-dc)
+    d1 = empty((2,2))
+    d1[0,0] = dc**2
+    d1[0,1] = -db**2
+    d1[1,0] = -dc**3
+    d1[1,1] = db**3
+    [A,B] = numpy.dot(d1,asarray([fb-fa-C*db,fc-fa-C*dc]).flatten())
+    A /= denom
+    B /= denom
+    radical = B*B-3*A*C
+    if radical < 0:  return None
+    if (A == 0): return None
+    xmin = a + (-B + sqrt(radical))/(3*A)
+    return xmin
+def _quadmin(a,fa,fpa,b,fb):
+    # finds the minimizer for a quadratic polynomial that goes through
+    #  the points (a,fa), (b,fb) with derivative at a of fpa
+    # f(x) = B*(x-a)^2 + C*(x-a) + D
+    D = fa
+    C = fpa
+    db = b-a*1.0
+    if (db==0): return None
+    B = (fb-D-C*db)/(db*db)
+    if (B <= 0): return None
+    xmin = a  - C / (2.0*B)
+    return xmin
+def line_search(f, myfprime, xk, pk, gfk, old_fval, old_old_fval,
+                args=(), c1=1e-4, c2=0.9, amax=50, self=None):
+    """Find alpha that satisfies strong Wolfe conditions.
+    :Parameters:
+        f : callable f(x,*args)
+            Objective function.
+        myfprime : callable f'(x,*args)
+            Objective function gradient (can be None).
+        xk : ndarray
+            Starting point.
+        pk : ndarray
+            Search direction.
+        gfk : ndarray
+            Gradient value for x=xk (xk being the current parameter
+            estimate).
+        args : tuple
+            Additional arguments passed to objective function.
+        c1 : float
+            Parameter for Armijo condition rule.
+        c2 : float
+            Parameter for curvature condition rule.
+    :Returns:
+        alpha0 : float
+            Alpha for which ``x_new = x0 + alpha * pk``.
+        fc : int
+            Number of function evaluations made.
+        gc : int
+            Number of gradient evaluations made.
+    :Notes:
+        Uses the line search algorithm to enforce strong Wolfe
+        conditions.  See Wright and Nocedal, 'Numerical Optimization',
+, pg. 59-60.
+        For the zoom phase it uses an algorithm by [...].
+    """
+    global _ls_fc, _ls_gc, _ls_ingfk
+    _ls_fc = 0
+    _ls_gc = 0
+    _ls_ingfk = None
+    def phi(alpha):
+        global _ls_fc
+        _ls_fc += 1
+        return f(xk+alpha*pk,*args)
+    if isinstance(myfprime,type(())):
+        def phiprime(alpha):
+            global _ls_fc, _ls_ingfk
+            _ls_fc += len(xk)+1
+            eps = myfprime[1]
+            fprime = myfprime[0]
+            newargs = (f,eps) + args
+            _ls_ingfk = fprime(xk+alpha*pk,*newargs)  # store for later use
+            return numpy.dot(_ls_ingfk,pk)
+    else:
+        fprime = myfprime
+        def phiprime(alpha):
+            global _ls_gc, _ls_ingfk
+            _ls_gc += 1
+            _ls_ingfk = fprime(xk+alpha*pk,*args)  # store for later use
+            return numpy.dot(_ls_ingfk,pk)
+    alpha0 = 0
+    phi0 = old_fval
+    derphi0 = numpy.dot(gfk,pk)
+    alpha1 = min(1.0,1.01*2*(phi0-old_old_fval)/derphi0)
+    if alpha1 == 0:
+        # This shouldn't happen. Perhaps the increment has slipped below
+        # machine precision?  For now, set the return variables skip the
+        # useless while loop, and raise warnflag=2 due to possible imprecision.
+        alpha_star = None
+        fval_star = old_fval
+        old_fval = old_old_fval
+        fprime_star = None
+    phi_a1 = phi(alpha1)
+    #derphi_a1 = phiprime(alpha1)  evaluated below
+    phi_a0 = phi0
+    derphi_a0 = derphi0
+    i = 1
+    maxiter = 10
+    while 1:         # bracketing phase
+        if alpha1 == 0:
+            break
+        if (phi_a1 > phi0 + c1*alpha1*derphi0) or \
+           ((phi_a1 >= phi_a0) and (i > 1)):
+            alpha_star, fval_star, fprime_star = \
+                        zoom(alpha0, alpha1, phi_a0,
+                             phi_a1, derphi_a0, phi, phiprime,
+                             phi0, derphi0, c1, c2)
+            break
+        derphi_a1 = phiprime(alpha1)
+        if (abs(derphi_a1) <= -c2*derphi0):
+            alpha_star = alpha1
+            fval_star = phi_a1
+            fprime_star = derphi_a1
+            break
+        if (derphi_a1 >= 0):
+            alpha_star, fval_star, fprime_star = \
+                        zoom(alpha1, alpha0, phi_a1,
+                             phi_a0, derphi_a1, phi, phiprime,
+                             phi0, derphi0, c1, c2)
+            break
+        alpha2 = 2 * alpha1   # increase by factor of two on each iteration
+        i = i + 1
+        alpha0 = alpha1
+        alpha1 = alpha2
+        phi_a0 = phi_a1
+        phi_a1 = phi(alpha1)
+        derphi_a0 = derphi_a1
+        # stopping test if lower function not found
+        if (i > maxiter):
+            alpha_star = alpha1
+            fval_star = phi_a1
+            fprime_star = None
+            break
+        if self._EARLYEXIT: # In case it gets stuck in this loop, 'exit' will work
+            alpha_star = alpha1
+            fval_star = phi_a1
+            fprime_star = None
+            break
+    if fprime_star is not None:
+        # fprime_star is a number (derphi) -- so use the most recently
+        # calculated gradient used in computing it derphi = gfk*pk
+        # this is the gradient at the next step no need to compute it
+        # again in the outer loop.
+        fprime_star = _ls_ingfk
+    return alpha_star, _ls_fc, _ls_gc, fval_star, old_fval, fprime_star
+def approx_fprime(xk,f,epsilon,*args):
+    f0 = f(*((xk,)+args))
+    grad = numpy.zeros((len(xk),), float)
+    ei = numpy.zeros((len(xk),), float)
+    for k in range(len(xk)):
+        ei[k] = epsilon
+        grad[k] = (f(*((xk+ei,)+args)) - f0)/epsilon
+        ei[k] = 0.0
+    return grad
+from _scipyoptimize import zoom, _cubicmin, _quadmin, _epsilon
+from _scipyoptimize import line_search, approx_fprime
 ##########################################################################

branches/alta/mystic-0.2a1/scipy_cg.py

-                      r404
+                      r418
 # Helper methods and classes
+_epsilon = sqrt(numpy.finfo(float).eps)
+def zoom(a_lo, a_hi, phi_lo, phi_hi, derphi_lo,
+         phi, derphi, phi0, derphi0, c1, c2):
+    maxiter = 10
+    i = 0
+    delta1 = 0.2  # cubic interpolant check
+    delta2 = 0.1  # quadratic interpolant check
+    phi_rec = phi0
+    a_rec = 0
+    while 1:
+        # interpolate to find a trial step length between a_lo and a_hi
+        # Need to choose interpolation here.  Use cubic interpolation and then if the
+        #  result is within delta * dalpha or outside of the interval bounded by a_lo or a_hi
+        #  then use quadratic interpolation, if the result is still too close, then use bisection
+        dalpha = a_hi-a_lo;
+        if dalpha < 0: a,b = a_hi,a_lo
+        else: a,b = a_lo, a_hi
+        # minimizer of cubic interpolant
+        #    (uses phi_lo, derphi_lo, phi_hi, and the most recent value of phi)
+        #      if the result is too close to the end points (or out of the interval)
+        #         then use quadratic interpolation with phi_lo, derphi_lo and phi_hi
+        #      if the result is stil too close to the end points (or out of the interval)
+        #         then use bisection
+        if (i > 0):
+            cchk = delta1*dalpha
+            a_j = _cubicmin(a_lo, phi_lo, derphi_lo, a_hi, phi_hi, a_rec, phi_rec)
+        if (i==0) or (a_j is None) or (a_j > b-cchk) or (a_j < a+cchk):
+            qchk = delta2*dalpha
+            a_j = _quadmin(a_lo, phi_lo, derphi_lo, a_hi, phi_hi)
+            if (a_j is None) or (a_j > b-qchk) or (a_j < a+qchk):
+                a_j = a_lo + 0.5*dalpha
+#                print "Using bisection."
+#            else: print "Using quadratic."
+#        else: print "Using cubic."
+        # Check new value of a_j
+        phi_aj = phi(a_j)
+        if (phi_aj > phi0 + c1*a_j*derphi0) or (phi_aj >= phi_lo):
+            phi_rec = phi_hi
+            a_rec = a_hi
+            a_hi = a_j
+            phi_hi = phi_aj
+        else:
+            derphi_aj = derphi(a_j)
+            if abs(derphi_aj) <= -c2*derphi0:
+                a_star = a_j
+                val_star = phi_aj
+                valprime_star = derphi_aj
+                break
+            if derphi_aj*(a_hi - a_lo) >= 0:
+                phi_rec = phi_hi
+                a_rec = a_hi
+                a_hi = a_lo
+                phi_hi = phi_lo
+            else:
+                phi_rec = phi_lo
+                a_rec = a_lo
+            a_lo = a_j
+            phi_lo = phi_aj
+            derphi_lo = derphi_aj
+        i += 1
+        if (i > maxiter):
+            a_star = a_j
+            val_star = phi_aj
+            valprime_star = None
+            break
+    return a_star, val_star, valprime_star
+def _cubicmin(a,fa,fpa,b,fb,c,fc):
+    # finds the minimizer for a cubic polynomial that goes through the
+    #  points (a,fa), (b,fb), and (c,fc) with derivative at a of fpa.
+    #
+    # if no minimizer can be found return None
+    #
+    # f(x) = A *(x-a)^3 + B*(x-a)^2 + C*(x-a) + D
+    C = fpa
+    D = fa
+    db = b-a
+    dc = c-a
+    if (db == 0) or (dc == 0) or (b==c): return None
+    denom = (db*dc)**2 * (db-dc)
+    d1 = empty((2,2))
+    d1[0,0] = dc**2
+    d1[0,1] = -db**2
+    d1[1,0] = -dc**3
+    d1[1,1] = db**3
+    [A,B] = numpy.dot(d1,asarray([fb-fa-C*db,fc-fa-C*dc]).flatten())
+    A /= denom
+    B /= denom
+    radical = B*B-3*A*C
+    if radical < 0:  return None
+    if (A == 0): return None
+    xmin = a + (-B + sqrt(radical))/(3*A)
+    return xmin
+def _quadmin(a,fa,fpa,b,fb):
+    # finds the minimizer for a quadratic polynomial that goes through
+    #  the points (a,fa), (b,fb) with derivative at a of fpa
+    # f(x) = B*(x-a)^2 + C*(x-a) + D
+    D = fa
+    C = fpa
+    db = b-a*1.0
+    if (db==0): return None
+    B = (fb-D-C*db)/(db*db)
+    if (B <= 0): return None
+    xmin = a  - C / (2.0*B)
+    return xmin
+def line_search(f, myfprime, xk, pk, gfk, old_fval, old_old_fval,
+                args=(), c1=1e-4, c2=0.9, amax=50, self=None):
+    """Find alpha that satisfies strong Wolfe conditions.
+    :Parameters:
+        f : callable f(x,*args)
+            Objective function.
+        myfprime : callable f'(x,*args)
+            Objective function gradient (can be None).
+        xk : ndarray
+            Starting point.
+        pk : ndarray
+            Search direction.
+        gfk : ndarray
+            Gradient value for x=xk (xk being the current parameter
+            estimate).
+        args : tuple
+            Additional arguments passed to objective function.
+        c1 : float
+            Parameter for Armijo condition rule.
+        c2 : float
+            Parameter for curvature condition rule.
+    :Returns:
+        alpha0 : float
+            Alpha for which ``x_new = x0 + alpha * pk``.
+        fc : int
+            Number of function evaluations made.
+        gc : int
+            Number of gradient evaluations made.
+    :Notes:
+        Uses the line search algorithm to enforce strong Wolfe
+        conditions.  See Wright and Nocedal, 'Numerical Optimization',
+, pg. 59-60.
+        For the zoom phase it uses an algorithm by [...].
+    """
+    global _ls_fc, _ls_gc, _ls_ingfk
+    _ls_fc = 0
+    _ls_gc = 0
+    _ls_ingfk = None
+    def phi(alpha):
+        global _ls_fc
+        _ls_fc += 1
+        return f(xk+alpha*pk,*args)
+    if isinstance(myfprime,type(())):
+        def phiprime(alpha):
+            global _ls_fc, _ls_ingfk
+            _ls_fc += len(xk)+1
+            eps = myfprime[1]
+            fprime = myfprime[0]
+            newargs = (f,eps) + args
+            _ls_ingfk = fprime(xk+alpha*pk,*newargs)  # store for later use
+            return numpy.dot(_ls_ingfk,pk)
+    else:
+        fprime = myfprime
+        def phiprime(alpha):
+            global _ls_gc, _ls_ingfk
+            _ls_gc += 1
+            _ls_ingfk = fprime(xk+alpha*pk,*args)  # store for later use
+            return numpy.dot(_ls_ingfk,pk)
+    alpha0 = 0
+    phi0 = old_fval
+    derphi0 = numpy.dot(gfk,pk)
+    alpha1 = min(1.0,1.01*2*(phi0-old_old_fval)/derphi0)
+    if alpha1 == 0:
+        # This shouldn't happen. Perhaps the increment has slipped below
+        # machine precision?  For now, set the return variables skip the
+        # useless while loop, and raise warnflag=2 due to possible imprecision.
+        alpha_star = None
+        fval_star = old_fval
+        old_fval = old_old_fval
+        fprime_star = None
+    phi_a1 = phi(alpha1)
+    #derphi_a1 = phiprime(alpha1)  evaluated below
+    phi_a0 = phi0
+    derphi_a0 = derphi0
+    i = 1
+    maxiter = 10
+    while 1:         # bracketing phase
+        if alpha1 == 0:
+            break
+        if (phi_a1 > phi0 + c1*alpha1*derphi0) or \
+           ((phi_a1 >= phi_a0) and (i > 1)):
+            alpha_star, fval_star, fprime_star = \
+                        zoom(alpha0, alpha1, phi_a0,
+                             phi_a1, derphi_a0, phi, phiprime,
+                             phi0, derphi0, c1, c2)
+            break
+        derphi_a1 = phiprime(alpha1)
+        if (abs(derphi_a1) <= -c2*derphi0):
+            alpha_star = alpha1
+            fval_star = phi_a1
+            fprime_star = derphi_a1
+            break
+        if (derphi_a1 >= 0):
+            alpha_star, fval_star, fprime_star = \
+                        zoom(alpha1, alpha0, phi_a1,
+                             phi_a0, derphi_a1, phi, phiprime,
+                             phi0, derphi0, c1, c2)
+            break
+        alpha2 = 2 * alpha1   # increase by factor of two on each iteration
+        i = i + 1
+        alpha0 = alpha1
+        alpha1 = alpha2
+        phi_a0 = phi_a1
+        phi_a1 = phi(alpha1)
+        derphi_a0 = derphi_a1
+        # stopping test if lower function not found
+        if (i > maxiter):
+            alpha_star = alpha1
+            fval_star = phi_a1
+            fprime_star = None
+            break
+        if self._EARLYEXIT: # In case it gets stuck in this loop, 'exit' will work
+            alpha_star = alpha1
+            fval_star = phi_a1
+            fprime_star = None
+            break
+    if fprime_star is not None:
+        # fprime_star is a number (derphi) -- so use the most recently
+        # calculated gradient used in computing it derphi = gfk*pk
+        # this is the gradient at the next step no need to compute it
+        # again in the outer loop.
+        fprime_star = _ls_ingfk
+    return alpha_star, _ls_fc, _ls_gc, fval_star, old_fval, fprime_star
+def approx_fprime(xk,f,epsilon,*args):
+    f0 = f(*((xk,)+args))
+    grad = numpy.zeros((len(xk),), float)
+    ei = numpy.zeros((len(xk),), float)
+    for k in range(len(xk)):
+        ei[k] = epsilon
+        grad[k] = (f(*((xk+ei,)+args)) - f0)/epsilon
+        ei[k] = 0.0
+    return grad
+from _scipyoptimize import zoom, _cubicmin, _quadmin, _epsilon
+from _scipyoptimize import line_search, approx_fprime
 ##########################################################################

branches/alta/mystic-0.2a1/scipy_ncg.py

-                      r417
+                      r418
 # Helper methods and classes
+def line_search_BFGS(f, xk, pk, gfk, old_fval, args=(), c1=1e-4, alpha0=1,\
+                     self=None):
+    """Minimize over alpha, the function ``f(xk+alpha pk)``.
+    Uses the interpolation algorithm (Armiijo backtracking) as suggested by
+    Wright and Nocedal in 'Numerical Optimization', 1999, pg. 56-57
+    :Returns: (alpha, fc, gc)
+    """
+    xk = atleast_1d(xk)
+    fc = 0
+    phi0 = old_fval # compute f(xk) -- done in past loop
+    phi_a0 = f(*((xk+alpha0*pk,)+args))
+    fc = fc + 1
+    derphi0 = numpy.dot(gfk,pk)
+    if (phi_a0 <= phi0 + c1*alpha0*derphi0):
+        return alpha0, fc, 0, phi_a0
+    # Otherwise compute the minimizer of a quadratic interpolant:
+    alpha1 = -(derphi0) * alpha0**2 / 2.0 / (phi_a0 - phi0 - derphi0 * alpha0)
+    phi_a1 = f(*((xk+alpha1*pk,)+args))
+    fc = fc + 1
+    if (phi_a1 <= phi0 + c1*alpha1*derphi0):
+        return alpha1, fc, 0, phi_a1
+    # Otherwise loop with cubic interpolation until we find an alpha which
+    # satifies the first Wolfe condition (since we are backtracking, we will
+    # assume that the value of alpha is not too small and satisfies the second
+    # condition.
+    while 1:       # we are assuming pk is a descent direction
+        factor = alpha0**2 * alpha1**2 * (alpha1-alpha0)
+        a = alpha0**2 * (phi_a1 - phi0 - derphi0*alpha1) - \
+            alpha1**2 * (phi_a0 - phi0 - derphi0*alpha0)
+        a = a / factor
+        b = -alpha0**3 * (phi_a1 - phi0 - derphi0*alpha1) + \
+            alpha1**3 * (phi_a0 - phi0 - derphi0*alpha0)
+        b = b / factor
+        alpha2 = (-b + numpy.sqrt(abs(b**2 - 3 * a * derphi0))) / (3.0*a)
+        phi_a2 = f(*((xk+alpha2*pk,)+args))
+        fc = fc + 1
+        if (phi_a2 <= phi0 + c1*alpha2*derphi0):
+            return alpha2, fc, 0, phi_a2
+        if (alpha1 - alpha2) > alpha1 / 2.0 or (1 - alpha2/alpha1) < 0.96:
+            alpha2 = alpha1 / 2.0
+        alpha0 = alpha1
+        alpha1 = alpha2
+        phi_a0 = phi_a1
+        phi_a1 = phi_a2
+        # In case it gets stuck in this loop, to make 'exit' work.
+        if self._EARLYEXIT:
+            return alpha2, fc, 0, phi_a2
+def approx_fhess_p(x0,p,fprime,epsilon,*args):
+    f2 = fprime(*((x0+epsilon*p,)+args))
+    f1 = fprime(*((x0,)+args))
+    return (f2 - f1)/epsilon
+_epsilon = sqrt(numpy.finfo(float).eps)
+# from scipy_bfgs
+def approx_fprime(xk,f,epsilon,*args):
+    xk = xk.copy()
+    f0 = f(*((xk,)+args))
+    grad = numpy.zeros((len(xk),), float)
+    ei = numpy.zeros((len(xk),), float)
+    for k in range(len(xk)):
+        ei[k] = epsilon
+        grad[k] = (f(*((xk+ei,)+args)) - f0)/epsilon
+        ei[k] = 0.0
+    return grad
+from _scipyoptimize import line_search_BFGS, _epsilon
+from _scipyoptimize import approx_fprime, approx_fhess_p
 #########################################################################

branches/alta/mystic-0.2a1/scipy_optimize.py

-                      r403
+                      r418
 from numpy import clip, squeeze
+abs = absolute
+from _scipy060optimize import brent #XXX: local copy to avoid dependency!
+from _scipyoptimize import _linesearch_powell
 from abstract_solver import AbstractSolver
 …
 ############################################################################
-def _linesearch_powell(func, p, xi, tol=1e-3):
-    # line-search algorithm using fminbound
-    #  find the minimium of the function
-    #  func(x0+ alpha*direc)
-    def myfunc(alpha):
-        return func(p + alpha * xi)
-    alpha_min, fret, iter, num = brent(myfunc, full_output=1, tol=tol)
-    xi = alpha_min*xi
-    return squeeze(fret), p+xi, xi
 class PowellDirectionalSolver(AbstractSolver):
     """

Note: See TracChangeset for help on using the changeset viewer.