Changeset 104

Timestamp:

02/20/09 18:36:22 (7 years ago)

Author:

mmckerns

Message:

moved fmin_powell to scipy_optimize;
removed scipy dependency;
cleanup of Solver documentation

Files:

• DEV_NOTES (modified) (1 diff)
• mystic/examples/test_rosenbrock3.py (modified) (1 diff)
• mystic/mystic/Make.mm (modified) (1 diff)
• mystic/mystic/__init__.py (modified) (1 diff)
• mystic/mystic/_scipy060optimize.py (added)
• mystic/mystic/differential_evolution.py (modified) (1 diff)
• mystic/mystic/fmin_powell.py (deleted)
• mystic/mystic/scipy_optimize.py (modified) (6 diffs)

DEV_NOTES

@@ -5 +5 @@
 */*/*py:
  - move "circle" to mystic.models, and move more of "mogi" to mystic.models
-mystic/mystic/fmin_powell.py:
+mystic/mystic/scipy_optimize.py:
  - derive an AbstractSolver class for all solvers to inherit from
- - add scipy solver "Brent", to _remove_ the scipy dependency (*)
- - migrate all contents into scipy_optimize.py
- - clean up & edit inline documentation
- - better way to deal with callback, direc, maxfun, maxiter (?)
-mystic/mystic/scipy_optimize.py:
- - add the "__all__" scipy trick
+ - better way to deal with direc, maxfun, maxiter (?)
 mystic/mystic/differential_evolution.py:
  - add ability to choose initial point

mystic/examples/test_rosenbrock3.py

@@ -29 +29 @@
 
    #from scipy.optimize import fmin_powell
-    from mystic.fmin_powell import fmin_powell, PowellDirectionalSolver
+    from mystic.scipy_optimize import fmin_powell, PowellDirectionalSolver
    #print fmin_powell(rosen,x0,retall=0,full_output=0)
     solver = PowellDirectionalSolver(len(x0))

mystic/mystic/Make.mm

@@ -29 +29 @@
     differential_evolution.py \
     scipy_optimize.py \
-    fmin_powell.py \
+    _scipy060optimize.py \
     termination.py \
     strategy.py \

mystic/mystic/__init__.py

@@ -16 +16 @@
 # solvers
 import differential_evolution, scipy_optimize
-import fmin_powell #FIXME: move to scipy_optimize after remove scipy dependency
 
 # strategies, termination conditions

mystic/mystic/differential_evolution.py

@@ -42 +42 @@
 
 """
+__all__ = ['DifferentialEvolutionSolver','DifferentialEvolutionSolver2']
 
 from mystic.tools import Null, wrap_function

mystic/mystic/scipy_optimize.py

@@ -4 +4 @@
 # (derives from optimize.py module by Travis E. Oliphant)
 #
-# adapted from scipy.optimize.fmin, scipy version 0.4.8
+# adapted scipy.optimize.fmin (from scipy version 0.4.8)
 # by Patrick Hung, Caltech.
 #
-# adapted to class (& added bounds)
-# updated to scipy version 0.6.0
+# adapted from function to class (& added bounds)
+# adapted scipy.optimize.fmin_powell
+# updated solvers to scipy version 0.6.0
 # by mmckerns@caltech.edu
 
 """
-Algorithm adapted from scipy.optimize, 
+Algorithms adapted from scipy.optimize, 
 
 fmin        ---      Nelder-Mead Simplex algorithm.
@@ -19 +20 @@
                         -- StepMonitor
                      
+fmin_powell ---      Powell's Directional Search method.
+                     Takes two additional input args
+                        -- EvaluationMonitor
+                        -- StepMonitor
+                     
 """
+__all__ = ['NelderMeadSimplexSolver','PowellDirectionalSolver',
+           'fmin','fmin_powell']
+
 
 from mystic.tools import Null, wrap_function
@@ -27 +36 @@
 from numpy import atleast_1d, eye, zeros, shape, \
      asarray, absolute, sqrt, Inf, asfarray
-from numpy import clip
+from numpy import clip, squeeze
 
 abs = absolute
+
+from _scipy060optimize import brent #XXX: local copy to avoid dependency!
+
 
 class NelderMeadSimplexSolver(object):
@@ -176 +188 @@
       ExtraArgs -- extra arguments for func.
 
-          """
+    Further Inputs:
+
+      callback -- an optional user-supplied function to call after each
+                  iteration.  It is called as callback(xk), where xk is the
+                  current parameter vector
+      disp -- non-zero to print convergence messages.
+      radius -- percentage change for initial simplex values
+
+"""
         # set arg names to scipy.optimize.fmin names; set fixed inputs
         x0 = self.population[0]
@@ -185 +205 @@
         callback=None  #user-supplied function, called after each step
         radius=0.05    #percentage change for initial simplex values
+        if kwds.has_key('callback'): callback = kwds['callback']
         if kwds.has_key('disp'): disp = kwds['disp']
         if kwds.has_key('radius'): radius = kwds['radius']
-        if kwds.has_key('callback'): callback = kwds['callback']
         #-------------------------------------------------------------
 
@@ -431 +451 @@
     return retlist
 
+############################################################################
+
+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(NelderMeadSimplexSolver): #FIXME: not a simplex solver
+    """
+    Powell Direction Search optimization adapted from scipy.optimize.fmin_powell.
+    """
+    
+    def __init__(self, dim):
+        """
+ Takes one initial input: 
+   dim      -- dimensionality of the problem
+        """
+        NelderMeadSimplexSolver.__init__(self,dim)
+        self._direc = None #FIXME: this is the easy way to return 'direc'...
+       #FIXME: NO SIMPLEX, so maybe best not to inherit __init__?
+
+
+    def Solve(self, func, termination,
+              maxiter=None, maxfun=None, sigint_callback=None,
+              EvaluationMonitor=Null, StepMonitor=Null, ExtraArgs=(), **kwds):
+        """Minimize a function using modified Powell's method.
+
+    Description:
+
+      Uses a modified Powell Directional Search algorithm to find
+      the minimum of function of one or more variables.
+
+    Inputs:
+
+      func -- the Python function or method to be minimized.
+      termination -- callable object providing termination conditions.
+
+    Additional Inputs:
+
+      maxiter -- the maximum number of iterations to perform.
+      maxfun -- the maximum number of function evaluations.
+      sigint_callback -- callback function for signal handler.
+      EvaluationMonitor -- a callable object that will be passed x, fval
+           whenever the cost function is evaluated.
+      StepMonitor -- a callable object that will be passed x, fval
+           after the end of a simplex iteration.
+      ExtraArgs -- extra arguments for func.
+
+    Further Inputs:
+
+      callback -- an optional user-supplied function to call after each
+                  iteration.  It is called as callback(xk), where xk is the
+                  current parameter vector
+      direc -- initial direction set
+      xtol -- line-search error tolerance.
+      disp -- non-zero to print convergence messages.
+
+"""
+
+        # set arg names to scipy.optimize.fmin_powell names; set fixed inputs
+        x0 = self.population[0]
+        args = ExtraArgs
+        full_output=1  #non-zero if fval and warnflag outputs are desired.
+        disp=0         #non-zero to print convergence messages.
+        retall=0       #non-zero to return all steps
+        direc=None
+        callback=None  #user-supplied function, called after each step
+        xtol=1e-4      #line-search error tolerance
+        if kwds.has_key('callback'): callback = kwds['callback']
+        if kwds.has_key('direc'): direc = kwds['direc']  #XXX: best interface?
+        if kwds.has_key('xtol'): xtol = kwds['xtol']
+        if kwds.has_key('disp'): disp = kwds['disp']
+        #-------------------------------------------------------------
+
+        import signal
+        import mystic.termination as detools
+        detools.EARLYEXIT = False
+
+        fcalls, func = wrap_function(func, args, EvaluationMonitor)
+        if self._useStrictRange:
+            x0 = self._setGuessWithinRangeBoundary(x0)
+            func = wrap_bounds(func, self._strictMin, self._strictMax)
+
+        def handler(signum, frame):
+            import inspect
+            print inspect.getframeinfo(frame)
+            print inspect.trace()
+            while 1:
+                s = raw_input(\
+"""
+ 
+ Enter sense switch.
+
+   sol: Write current best solution.
+   cont: Continue calculation.
+   call: Executes sigint_callback [%s].
+   exit: Exits with current best solution.
+
+ >>> """ % sigint_callback)
+                if s.lower() == 'sol': 
+                    print "sw1."
+                    print self.bestSolution
+                elif s.lower() == 'cont': 
+                    return
+                elif s.lower() == 'call': 
+                    # sigint call_back
+                    if sigint_callback is not None:
+                        sigint_callback(self.bestSolution)
+                elif s.lower() == 'exit': 
+                    detools.EARLYEXIT = True
+                    return
+                else:
+                    print "unknown option : %s ", s
+
+        self.signal_handler = handler
+
+        if self._handle_sigint: signal.signal(signal.SIGINT, self.signal_handler)
+        #-------------------------------------------------------------
+
+        x = asfarray(x0).flatten()
+        if retall:
+            allvecs = [x]
+        N = len(x) #XXX: this should be equal to self.nDim
+        rank = len(x.shape)
+        if not -1 < rank < 2:
+            raise ValueError, "Initial guess must be a scalar or rank-1 sequence."
+        if maxiter is None:
+            maxiter = N * 1000
+        if maxfun is None:
+            maxfun = N * 1000
+        self._maxiter = maxiter #XXX: better to just copy the code?
+        self._maxfun = maxfun   #XXX: better to just copy the code?
+
+        if direc is None:
+            direc = eye(N, dtype=float)
+        else:
+            direc = asarray(direc, dtype=float)
+        fval = squeeze(func(x))
+        x1 = x.copy()
+
+        self._direc = direc #XXX: instead, use a monitor?
+        self.bestSolution = x
+        self.bestEnergy = fval
+        self.population = [x]    #XXX: pointless, if simplex not used
+        self.popEnergy = [fval]  #XXX: pointless, if simplex not used
+        self.energy_history.append(self.bestEnergy)
+        StepMonitor(x,fval) # get initial values
+
+        iter = 0;
+        ilist = range(N)
+
+        CONTINUE = True
+        while CONTINUE:
+            fx = fval
+            bigind = 0
+            delta = 0.0
+            for i in ilist:
+                direc1 = direc[i]
+                fx2 = fval
+                fval, x, direc1 = _linesearch_powell(func, x, direc1, tol=xtol*100)
+                if (fx2 - fval) > delta:
+                    delta = fx2 - fval
+                    bigind = i
+
+            iter += 1
+            if callback is not None:
+                callback(x)
+            if retall:
+                allvecs.append(x)
+
+            self.energy_history.append(fval) #XXX: the 'best' for now...
+            if detools.EARLYEXIT or termination(self): CONTINUE = False #break
+            elif fcalls[0] >= maxfun: CONTINUE = False #break
+            elif iter >= maxiter: CONTINUE = False #break
+
+            else: # Construct the extrapolated point
+                direc1 = x - x1
+                x2 = 2*x - x1
+                x1 = x.copy()
+                fx2 = squeeze(func(x2))
+    
+                if (fx > fx2):
+                    t = 2.0*(fx+fx2-2.0*fval)
+                    temp = (fx-fval-delta)
+                    t *= temp*temp
+                    temp = fx-fx2
+                    t -= delta*temp*temp
+                    if t < 0.0:
+                        fval, x, direc1 = _linesearch_powell(func, x, direc1, tol=xtol*100)
+                        direc[bigind] = direc[-1]
+                        direc[-1] = direc1
+
+                self.energy_history[-1] = fval #...update to 'best' energy
+
+            self._direc = direc #XXX: instead, use a monitor?
+            self.bestSolution = x
+            self.bestEnergy = fval
+            self.population = [x]    #XXX: pointless, if simplex not used
+            self.popEnergy = [fval]  #XXX: pointless, if simplex not used
+            StepMonitor(x,fval) # get ith values; #XXX: should be [x],[fval] ?
+    
+        self.generations = iter
+        signal.signal(signal.SIGINT,signal.default_int_handler)
+
+        # code below here is dead, unless disp!=0
+        warnflag = 0
+
+        if fcalls[0] >= maxfun:
+            warnflag = 1
+            if disp:
+                print "Warning: Maximum number of function evaluations has "\
+                      "been exceeded."
+        elif iter >= maxiter:
+            warnflag = 2
+            if disp:
+                print "Warning: Maximum number of iterations has been exceeded"
+        else:
+            if disp:
+                print "Optimization terminated successfully."
+                print "         Current function value: %f" % fval
+                print "         Iterations: %d" % iter
+                print "         Function evaluations: %d" % fcalls[0]
+    
+        x = squeeze(x)
+
+        if full_output:
+            retlist = x, fval, direc, iter, fcalls[0], warnflag
+            if retall:
+                retlist += (allvecs,)
+        else:
+            retlist = x
+            if retall:
+                retlist = (x, allvecs)
+
+        return #retlist
+
+
+def fmin_powell(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None,
+                maxfun=None, full_output=0, disp=1, retall=0, callback=None,
+                direc=None):
+    """fmin_powell using the 'original' scipy.optimize.fmin_powell interface"""
+
+    #FIXME: need to resolve "direc"
+    #        - should just pass 'direc', and then hands-off ?  How return it ?
+
+    from mystic.tools import Sow
+    stepmon = Sow()
+    evalmon = Sow()
+    from mystic.termination import NormalizedChangeOverGeneration as NCOG
+
+    solver = PowellDirectionalSolver(len(x0))
+    solver.SetInitialPoints(x0)
+   #solver.enable_signal_handler()
+    solver.Solve(func,termination=NCOG(ftol),\
+                 maxiter=maxiter,maxfun=maxfun,\
+                 EvaluationMonitor=evalmon,StepMonitor=stepmon,\
+                 xtol=xtol, callback=callback, \
+                 disp=disp, direc=direc)   #XXX: last two lines use **kwds
+    solution = solver.Solution()
+
+    # code below here pushes output to scipy.optimize.fmin_powell interface
+   #x = list(solver.bestSolution)
+    x = solver.bestSolution
+    fval = solver.bestEnergy
+    warnflag = 0
+    fcalls = len(evalmon.x)
+    iterations = len(stepmon.x) - 1
+    allvecs = []
+    for i in range(len(stepmon.x)):
+       #allvecs.append(list(stepmon.x[i]))
+        allvecs.append(stepmon.x[i])
+    direc = solver._direc #FIXME: better way to get direc from Solve() ?
+
+    if fcalls >= solver._maxfun:
+        warnflag = 1
+    elif iterations >= solver._maxiter:
+        warnflag = 2
+
+    x = squeeze(x) #FIXME: write squeezed x to stepmon instead?
+
+    if full_output:
+        retlist = x, fval, direc, iterations, fcalls, warnflag
+        if retall:
+            retlist += (allvecs,)
+    else:
+        retlist = x
+        if retall:
+            retlist = (x, allvecs)
+
+    return retlist
+
 
 if __name__=='__main__':