Changeset 422

Timestamp:

09/08/10 22:54:59 (6 years ago)

Author:

mmckerns

Message:

added several new termination conditions;
added tools.wrap_nested function;
uncommented local copies of scipy gradient solvers
minor bugfix in abstract_solver warning

Location:

mystic

Files:

• examples/test_rosenbrock2.py (modified) (1 diff)
• examples/test_rosenbrock3.py (modified) (1 diff)
• mystic/Make.mm (modified) (2 diffs)
• mystic/_scipy060optimize.py (modified) (6 diffs)
• mystic/_scipyoptimize.py (added)
• mystic/abstract_solver.py (modified) (2 diffs)
• mystic/linesearch.py (added)
• mystic/termination.py (modified) (2 diffs)
• mystic/tools.py (modified) (2 diffs)

mystic/examples/test_rosenbrock2.py

@@ -29 +29 @@
     from mystic.termination import CandidateRelativeTolerance as CRT
 
-    from scipy.optimize import fmin
+   #from scipy.optimize import fmin
     from mystic.scipy_optimize import fmin, NelderMeadSimplexSolver
-    print fmin(rosen,x0,retall=0,full_output=0,maxiter=121)
+   #print fmin(rosen,x0,retall=0,full_output=0,maxiter=121)
     solver = NelderMeadSimplexSolver(len(x0))
     solver.SetInitialPoints(x0)

mystic/examples/test_rosenbrock3.py

@@ -29 +29 @@
     from mystic.termination import NormalizedChangeOverGeneration as NCOG
 
-    from scipy.optimize import fmin_powell
+   #from scipy.optimize import fmin_powell
     from mystic.scipy_optimize import fmin_powell, PowellDirectionalSolver
-    print fmin_powell(rosen,x0,retall=0,full_output=0,maxiter=14)
+   #print fmin_powell(rosen,x0,retall=0,full_output=0)#,maxiter=14)
     solver = PowellDirectionalSolver(len(x0))
     solver.SetInitialPoints(x0)
     solver.SetStrictRanges(min,max)
-    solver.SetEvaluationLimits(maxiter=13)
+   #solver.SetEvaluationLimits(maxiter=13)
     solver.enable_signal_handler()
     solver.Solve(rosen,termination=NCOG(tolerance=1e-4),StepMonitor=stepmon,disp=1)

mystic/mystic/Make.mm

@@ -28 +28 @@
     _genSow.py \
     _scipy060optimize.py \
+    _scipyoptimize.py \
     abstract_map_solver.py \
     abstract_nested_solver.py \
@@ -35 +36 @@
     forward_model.py \
     helputil.py \
+    linesearch.py \
     mystic_math.py \
     nested.py \

mystic/mystic/_scipy060optimize.py

@@ -6 +6 @@
 # guarantee implied provided you keep this notice in all copies.
 # *****END NOTICE************
-
-# Minimization routines
-# (removed: fmin_bfgs, fmin_cg)
 """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', 'fmin_cg', 'fmin_bfgs',
+#           'fminbound','brent', 'golden','bracket','rosen','rosen_der',
+#           'rosen_hess', 'rosen_hess_prod', 'brute', 'approx_fprime',
+#           'line_search', 'check_grad']
 
 import numpy
 from numpy import atleast_1d, eye, mgrid, argmin, zeros, shape, empty, \
      squeeze, isscalar, vectorize, asarray, absolute, sqrt, Inf, asfarray, isinf
-#import linesearch
+import linesearch
 
 # These have been copied from Numeric's MLab.py
@@ -622 +619 @@
     return (f2 - f1)/epsilon
 
-'''
 def fmin_bfgs(f, x0, fprime=None, args=(), gtol=1e-5, norm=Inf,
               epsilon=_epsilon, maxiter=None, full_output=0, disp=1,
@@ -808 +804 @@
 
     return retlist
-'''
-
-'''
+
 def fmin_cg(f, x0, fprime=None, args=(), gtol=1e-5, norm=Inf, epsilon=_epsilon,
               maxiter=None, full_output=0, disp=1, retall=0, callback=None):
@@ -979 +973 @@
 
     return retlist
-'''
 
 def fmin_ncg(f, x0, fprime, fhess_p=None, fhess=None, args=(), avextol=1e-5,
@@ -2042 +2035 @@
     algor.append('Powell Direction Set Method.')
 
-#   print
-#   print "Nonlinear CG"
-#   print "============"
-#   start = time.time()
-#   x = fmin_cg(rosen, x0, fprime=rosen_der, maxiter=200)
-#   print x
-#   times.append(time.time() - start)
-#   algor.append('Nonlinear CG     \t')
-
-#   print
-#   print "BFGS Quasi-Newton"
-#   print "================="
-#   start = time.time()
-#   x = fmin_bfgs(rosen, x0, fprime=rosen_der, maxiter=80)
-#   print x
-#   times.append(time.time() - start)
-#   algor.append('BFGS Quasi-Newton\t')
-
-#   print
-#   print "BFGS approximate gradient"
-#   print "========================="
-#   start = time.time()
-#   x = fmin_bfgs(rosen, x0, gtol=1e-4, maxiter=100)
-#   print x
-#   times.append(time.time() - start)
-#   algor.append('BFGS without gradient\t')
-
+    print
+    print "Nonlinear CG"
+    print "============"
+    start = time.time()
+    x = fmin_cg(rosen, x0, fprime=rosen_der, maxiter=200)
+    print x
+    times.append(time.time() - start)
+    algor.append('Nonlinear CG     \t')
+
+    print
+    print "BFGS Quasi-Newton"
+    print "================="
+    start = time.time()
+    x = fmin_bfgs(rosen, x0, fprime=rosen_der, maxiter=80)
+    print x
+    times.append(time.time() - start)
+    algor.append('BFGS Quasi-Newton\t')
+
+    print
+    print "BFGS approximate gradient"
+    print "========================="
+    start = time.time()
+    x = fmin_bfgs(rosen, x0, gtol=1e-4, maxiter=100)
+    print x
+    times.append(time.time() - start)
+    algor.append('BFGS without gradient\t')
 
     print
@@ -2079 +2071 @@
     algor.append('Newton-CG with hessian product')
 
-
     print
     print "Newton-CG with full Hessian"

mystic/mystic/abstract_solver.py

@@ -131 +131 @@
         import mystic.termination
         self._EARLYEXIT       = mystic.termination.EARLYEXIT
-
+        return
 
     def Solution(self):
@@ -283 +283 @@
                     return
                 else:
-                    print "unknown option : %s ", s
+                    print "unknown option : %s" % s
             return
         self.signal_handler = handler

mystic/mystic/termination.py

@@ -8 +8 @@
 from numpy import absolute
 abs = absolute
+Inf = numpy.Inf
 
 # a module level singleton.
@@ -66 +67 @@
     return _CandidateRelativeTolerance
 
+def SolutionImprovement(tolerance = 1e-5):  
+    """sum of change in each parameter is < tolerance:
+
+sum(abs(last_params - current_params)) <= tolerance"""
+    def _SolutionImprovement(inst):
+        update = inst.bestSolution - inst.trialSolution #XXX: if inf - inf ?
+        answer = numpy.add.reduce(abs(update)) <= tolerance
+        return answer
+    return _SolutionImprovement
+
+def NormalizedCostTarget(fval = None, tolerance = 1e-6, generations = 30):
+    """normalized absolute difference from given cost value is < tolerance:
+(if fval is not provided, then terminate when no improvement over g iterations)
+
+abs(cost[-1] - fval)/fval <= tolerance *or* (cost[-1] - cost[-g]) = 0 """
+    #NOTE: modified from original behavior
+    #  original --> if generations: then return cost[-g] - cost[-1] < 0
+    #           --> else: return fval != 0 and abs((best - fval)/fval) < tol
+    def _NormalizedCostTarget(inst):
+         if generations and fval == None:
+             hist = inst.energy_history
+             lg = len(hist)
+             #XXX: throws error when hist is shorter than generations ?
+             return lg > generations and (hist[-generations]-hist[-1]) < 0
+         if not generations and fval == None: return True
+         return abs(inst.bestEnergy-fval) <= abs(tolerance * fval)
+    return _NormalizedCostTarget
+
+def VTRChangeOverGenerations(ftol = 0.005, gtol = 1e-6, generations = 30):
+    """change in cost is < gtol over a number of generations,
+or cost of last iteration is < ftol:
+
+cost[-g] - cost[-1] < gtol, where g=generations *or* cost[-1] < ftol."""
+    def _VTRChangeOverGenerations(inst):
+         hist = inst.energy_history
+         lg = len(hist)
+         #XXX: throws error when hist is shorter than generations ?
+         return (lg > generations and (hist[-generations]-hist[-1]) < gtol)\
+                or ( hist[-1] < ftol )
+    return _VTRChangeOverGenerations
+
+def PopulationSpread(tolerance=1e-6):
+    """normalized absolute deviation from best candidate is < tolerance:
+
+abs(params - params[0]) < tolerance"""
+    def _PopulationSpread(inst):
+         sim = numpy.array(inst.population)
+         #if not len(sim[1:]):
+         #    print "Warning: Invalid termination condition (nPop < 2)"
+         #    return True
+         return all(abs(sim - sim[0]) <= abs(tolerance * sim[0]))
+    return _PopulationSpread
+
+def GradientNormTolerance(tolerance=1e-5, norm=Inf): 
+    """gradient norm is < tolerance, given user-supplied norm:
+
+sum( abs(gradient)**norm )**(1.0/norm) < tolerance"""
+    def _GradientNormTolerance(inst):
+        try:
+            gfk = inst.gfk #XXX: need to ensure that gfk is an array ?
+        except:
+            print "Warning: Invalid termination condition (no gradient)"
+            return True
+        if norm == Inf:
+            gnorm = numpy.amax(abs(gfk))
+        elif norm == -Inf:
+            gnorm = numpy.amin(abs(gfk))
+        else: #XXX: throws error when norm = 0.0
+           #XXX: as norm > large, gnorm approaches amax(abs(gfk)) --> then inf
+           #XXX: as norm < -large, gnorm approaches amin(abs(gfk)) --> then -inf
+            gnorm = numpy.sum(abs(gfk)**norm,axis=0)**(1.0/norm)
+        return gnorm <= tolerance
+    return _GradientNormTolerance
+
 # end of file

mystic/mystic/tools.py

@@ -21 +21 @@
     - CustomSow: A customizable 'n-variable' version of the basic Sow
     - random_seed: sets the seed for calls to 'random()'
+    - wrap_nested: nest a function call within a function object
     - wrap_function: bind an EvaluationMonitor and an evaluation counter
         to a function object
@@ -291 +292 @@
     return
 
+def wrap_nested(function, inner_function):
+    """nest a function call within a function object
+
+This is useful for nesting a constraints function in a cost function;
+thus, the constraints will be enforced at every cost function evaluation.
+    """
+    def function_wrapper(x):
+        _x = x[:] #XXX: trouble if x not a list or ndarray... maybe "deepcopy"?
+        return function(inner_function(_x))
+    return function_wrapper
+
 def wrap_function(function, args, EvaluationMonitor):
     """bind an EvaluationMonitor and an evaluation counter