Changeset 87

Timestamp:

02/05/09 20:28:35 (7 years ago)

Author:

mmckerns

Message:

added tools.wrap_bounds; added range bounds to fmin (modified from park.simplex)

Location:

python

Files:

• nmtools.py (modified) (1 diff)
• scipy_optimize_fmin.py (modified) (10 diffs)
• tools.py (modified) (1 diff)

python/nmtools.py

@@ -16 +16 @@
          sim = inst.population
          fsim = inst.popEnergy
-         return (max(numpy.ravel(abs(sim[1:]-sim[0]))) <= xtol \
-                and max(abs(fsim[0]-fsim[1:])) <= ftol)
+         #FIXME: abs(inf - inf) will raise a warning...
+         errdict = numpy.seterr(invalid='ignore') #FIXME: turn off warning 
+         answer = (max(numpy.ravel(abs(sim[1:]-sim[0]))) <= xtol \
+                  and max(abs(fsim[0]-fsim[1:])) <= ftol)
+         numpy.seterr(invalid=errdict['invalid']) #FIXME: turn on warnings
+         return answer
     return _
 

python/scipy_optimize_fmin.py

@@ -5 +5 @@
 # adapted from scipy.optimize, scipy version 0.4.8
 # by Patrick Hung, Caltech.
-# adapted to class by mmckerns@caltech.edu
+# adapted to class (& added bounds)
+# by mmckerns@caltech.edu
 
 """
@@ -18 +19 @@
 
 from mystic.tools import Null, wrap_function
+from mystic.tools import wrap_bounds
 
 import numpy
 from numpy import atleast_1d, eye, zeros, shape, \
      asarray, absolute, sqrt, Inf, asfarray
+from numpy import clip
 
 abs = absolute
@@ -59 +62 @@
 
     def SetStrictRanges(self, min, max):
-        raise NotImplementedError, "bound constraints not implemented"
+        self._useStrictRange = True
+        self._strictMin = min
+        self._strictMax = max
+        return
+
+    def setSimplexWithinRangeBoundary(self, x0, radius):
+        """ensure that initial simplex is set within bounds"""
+        #code modified from park-1.2/park/simplex.py (version 1257)
+        if self._useStrictRange:
+            lo = asarray(self._strictMin)
+            hi = asarray(self._strictMax)
+            # crop x0 at bounds
+            x0[x0<lo] = lo[x0<lo]
+            x0[x0>hi] = hi[x0>hi]
+
+        val = x0*(1+radius)
+        val[val==0] = radius
+        if not self._useStrictRange:
+            return x0, val
+
+        radius = clip(radius,0,0.5)
+        # rescale val by bounded range...
+        # (increases fit for tight bounds; makes worse[?] for large bounds)
+        bounded = ~numpy.isinf(lo) & ~numpy.isinf(hi)
+        val[bounded] = x0[bounded] + (hi[bounded]-lo[bounded])*radius
+        # crop val at bounds
+        val[val<lo] = lo[val<lo]
+        val[val>hi] = hi[val>hi]
+        # handle collisions (when val[i] == x0[i])
+        collision = val==x0
+        if numpy.any(collision):
+            rval = x0*(1-radius)
+            rval[rval==0] = -radius
+            rval[bounded] = x0[bounded] - (hi[bounded]-lo[bounded])*radius
+            val[collision] = rval[collision]
+        # make tolerance relative for bounded parameters
+     #  tol = numpy.ones(x0.shape)*xtol
+     #  tol[bounded] = (hi[bounded]-lo[bounded])*xtol
+     #  xtol = tol
+        return x0, val
 
     def UpdateGenealogyRecords(self, id, newchild):
@@ -128 +170 @@
         disp=0         #non-zero to print convergence messages.
         retall=0
+        radius=0.05    #percentage change for initial simplex values
         if kwds.has_key('disp'): disp = kwds['disp']
+        if kwds.has_key('radius'): radius = kwds['radius']
         #-------------------------------------------------------------
 
@@ -136 +180 @@
 
         fcalls, func = wrap_function(func, args, EvaluationMonitor)
+        if self._useStrictRange:
+            func = wrap_bounds(func, self._strictMin, self._strictMax)
 
         def handler(signum, frame):
@@ -195 +241 @@
             allvecs = [sim[0]]
         fsim[0] = func(x0)
-        nonzdelt = 0.05
-        zdelt = 0.00025
+
+        #--- ensure initial simplex is within bounds ---
+        x0,val = self.setSimplexWithinRangeBoundary(x0,radius)
+        #--- end bounds code ---
         for k in range(0,N):
             y = numpy.array(x0,copy=True)
-            if y[k] != 0:
-                y[k] = (1+nonzdelt)*y[k]
-            else:
-                y[k] = zdelt
-
+            y[k] = val[k]
             sim[k+1] = y
             f = func(y)
@@ -209 +253 @@
     
         ind = numpy.argsort(fsim)
-        fsim = numpy.take(fsim,ind)  # sort so sim[0,:] has the lowest function value
+        fsim = numpy.take(fsim,ind)
+        # sort so sim[0,:] has the lowest function value
         sim = numpy.take(sim,ind,0)
-        self.bestSolution = sim[0] #XXX: goes here?
-        self.bestEnergy = min(fsim) #XXX: goes here?
-        self.population = sim #XXX: OK?
-        self.popEnergy = fsim #XXX: OK?
+        self.bestSolution = sim[0]
+        self.bestEnergy = min(fsim)
+        self.population = sim
+        self.popEnergy = fsim
 
         iterations = 1
@@ -276 +321 @@
                 allvecs.append(sim[0])
 
-            self.bestSolution = sim[0] #XXX: goes here?
-            self.bestEnergy = min(fsim) #XXX: goes here?
-            self.population = sim #XXX: OK?
-            self.popEnergy = fsim #XXX: OK?
+            self.bestSolution = sim[0]
+            self.bestEnergy = min(fsim)
+            self.population = sim
+            self.popEnergy = fsim
             self.energy_history.append(self.bestEnergy)
 
@@ -326 +371 @@
     stepmon = Sow()
     evalmon = Sow()
-    from nmtools import IterationRelativeError as IRE
+    from mystic.nmtools import IterationRelativeError as IRE
 
     solver = NelderMeadSimplexSolver(len(x0))
@@ -375 +420 @@
     algor = []
     x0 = [0.8,1.2,0.7]
+   #x0 = [0.8,1.2,1.7]                  #... better when using "bad" range
+    min = [-0.999, -0.999, 0.999]       #XXX: behaves badly when large range
+    max = [200.001, 100.001, numpy.inf] #... for >=1 x0 out of bounds; (up xtol)
+  # min = [-0.999, -0.999, -0.999]
+  # max = [200.001, 100.001, numpy.inf]
+ #  min = [-0.999, -0.999, 0.999]
+ #  max = [2.001, 1.001, 1.001]
     print "Nelder-Mead Simplex"
     print "==================="
     start = time.time()
-    from tools import VerboseSow
+    from mystic.tools import VerboseSow
     stepmon = VerboseSow(10)
-    from nmtools import IterationRelativeError as IRE
+    from mystic.nmtools import IterationRelativeError as IRE
 
    #print fmin(__rosen,x0,retall=0,full_output=0)
     solver = NelderMeadSimplexSolver(len(x0))
     solver.SetInitialPoints(x0)
+    solver.SetStrictRanges(min,max)
     solver.enable_signal_handler()
-    solver.Solve(__rosen,termination=IRE(),StepMonitor=stepmon)
+    solver.Solve(__rosen,termination=IRE(xtol=1e-5),StepMonitor=stepmon)
     print solver.Solution()
 

python/tools.py

@@ -189 +189 @@
     return ncalls, function_wrapper
 
+def wrap_bounds(function, min=None, max=None):
+    from numpy import asarray, any, inf
+    bounds = True
+    if min is not None and max is not None: #has upper & lower bound
+        min = asarray(min)
+        max = asarray(max)
+    elif min is not None: #has lower bound
+        min = asarray(min)
+        max = asarray([inf for i in min])
+    elif max is not None: #has upper bound
+        max = asarray(max)
+        min = asarray([-inf for i in max])
+    else: #not bounded
+        bounds = False
+    if bounds:
+        def function_wrapper(x):
+            if any((x<min)|(x>max)): #if violate bounds, evaluate as inf
+                return inf
+            return function(x)
+    else:
+        def function_wrapper(x):
+            return function(x)
+    return function_wrapper
+
 def wrap_cf(CF, REG=None, cfmult = 1.0, regmult = 0.0):
     def _(*args, **kwargs):