Changeset 219

Timestamp:

04/29/10 14:23:44 (6 years ago)

Author:

mmckerns

Message:

bugfix for #102, DE uses initial population;
SOW now captures initial pop,E for each solver;
consistent maxiter termination per solver & w/ scipy

Location:

mystic

Files:

• examples/test_rosenbrock.py (modified) (6 diffs)
• examples/test_rosenbrock2.py (modified) (1 diff)
• examples/test_rosenbrock3.py (modified) (1 diff)
• mystic/abstract_solver.py (modified) (3 diffs)
• mystic/differential_evolution.py (modified) (8 diffs)
• mystic/scipy_optimize.py (modified) (3 diffs)
• mystic/tools.py (modified) (2 diffs)

mystic/examples/test_rosenbrock.py

@@ -16 +16 @@
 from mystic.termination import ChangeOverGeneration, VTR
 from mystic.models import rosen
-from mystic import Sow
+from mystic import Sow, VerboseSow
 
 import random
@@ -23 +23 @@
 ND = 3
 NP = 30
+#MAX_GENERATIONS = 29
 MAX_GENERATIONS = 99999
 
@@ -32 +33 @@
 
     solver.Solve(rosen, termination = VTR(0.0001), \
-                 CrossProbability=0.5, ScalingFactor=0.6)
+                 CrossProbability=0.5, ScalingFactor=0.6, disp=1)
 
     solution = solver.Solution()
@@ -44 +45 @@
     print "without bounds..."
     from timeit import Timer
+    print "Differential Evolution"
+    print "======================"
     t = Timer("main()", "from __main__ import main")
     timetaken =  t.timeit(number=1)
@@ -58 +61 @@
     esow= Sow()
     ssow= Sow()
+   #ssow= VerboseSow(1)
 
  #  xinit = [random.random() for j in range(ND)]
@@ -75 +79 @@
     solver.SetEvaluationLimits(maxiter=MAX_GENERATIONS)
     solver.Solve(rosen, VTR(0.0001), EvaluationMonitor=esow, StepMonitor=ssow,\
-                 CrossProbability=0.5, ScalingFactor=0.6)
+                 CrossProbability=0.5, ScalingFactor=0.6, disp=1)
     sol = solver.Solution()
     print sol

mystic/examples/test_rosenbrock2.py

@@ -24 +24 @@
     print "==================="
     start = time.time()
-    from mystic.tools import Sow #VerboseSow
+    from mystic.tools import Sow, VerboseSow
+   #stepmon = VerboseSow(1)
     stepmon = Sow() #VerboseSow(10)
     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)
+    print fmin(rosen,x0,retall=0,full_output=0,maxiter=121)
     solver = NelderMeadSimplexSolver(len(x0))
     solver.SetInitialPoints(x0)
     solver.SetStrictRanges(min,max)
+    solver.SetEvaluationLimits(maxiter=146)
     solver.enable_signal_handler()
     solver.Solve(rosen,termination=CRT(xtol=4e-5),StepMonitor=stepmon,disp=1)

mystic/examples/test_rosenbrock3.py

@@ -24 +24 @@
     print "==========================="
     start = time.time()
-    from mystic.tools import Sow #VerboseSow
-    stepmon = Sow() #VerboseSow(1)
+    from mystic.tools import Sow, VerboseSow
+   #stepmon = VerboseSow(1)
+    stepmon = Sow() #VerboseSow(10)
     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)
+    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.enable_signal_handler()
     solver.Solve(rosen,termination=NCOG(tolerance=1e-4),StepMonitor=stepmon,disp=1)

mystic/mystic/abstract_solver.py

@@ -116 +116 @@
 
         self._init_popEnergy  = 1.0E20 #XXX: or numpy.inf?
-        self.popEnergy        = [0.0] * NP
+        self.popEnergy        = [self._init_popEnergy] * NP
         self.population       = [[0.0 for i in range(dim)] for j in range(NP)]
         self.energy_history   = []
@@ -211 +211 @@
             for j in range(self.nDim):
                 self.population[i][j] = random.uniform(min[j],max[j])
-            self.popEnergy[i] = self._init_popEnergy
 
     def SetMultinormalInitialPoints(self, mean, var = None):
@@ -236 +235 @@
         for i in range(self.nPop):
             self.population[i] = multivariate_normal(mean, var).tolist()
-            self.popEnergy[i] = self._init_popEnergy
         return
 

mystic/mystic/differential_evolution.py

@@ -249 +249 @@
         self.scale = ScalingFactor
 
-        self.bestEnergy = self._init_popEnergy
+        #set initial solution and energy by running a single iteration
+        self.bestSolution = self.population[0][:]
+        self.bestEnergy = self.popEnergy[0]
+        for candidate in range(self.nPop):
+            # generate trialSolution (within valid range)
+            self.trialSolution[:] = self.population[candidate][:]
+            trialEnergy = costfunction(self.trialSolution)
+
+            if trialEnergy < self.popEnergy[candidate]:
+                # New low for this candidate
+                self.popEnergy[candidate] = trialEnergy
+                self.population[candidate][:] = self.trialSolution[:]
+                self.UpdateGenealogyRecords(candidate, self.trialSolution[:])
+
+                # Check if all-time low
+                if trialEnergy < self.bestEnergy:
+                    self.bestEnergy = trialEnergy
+                    self.bestSolution[:] = self.trialSolution[:]
+                            
+        self.energy_history.append(self.bestEnergy)
+        StepMonitor(self.bestSolution[:], self.bestEnergy)
+        self.generations = 0  #XXX: above currently *not* counted as an iteration
+        if callback is not None:
+            callback(self.bestSolution)
          
         if self._maxiter is None:
@@ -256 +279 @@
             self._maxfun = self.nDim * self.nPop * 1000 #XXX: set better defaults?
 
-        generation = 0
-        for generation in range(self._maxiter):
-            StepMonitor(self.bestSolution[:], self.bestEnergy)
+        #run for generations <= maxiter
+        for generation in range(self._maxiter - self.generations):
             if fcalls[0] >= self._maxfun: break
             for candidate in range(self.nPop):
@@ -277 +299 @@
                             
             self.energy_history.append(self.bestEnergy)
-
+            StepMonitor(self.bestSolution[:], self.bestEnergy)
+            self.generations += 1
             if callback is not None:
                 callback(self.bestSolution)
@@ -283 +306 @@
             if self._EARLYEXIT or termination(self):
                 break
-
-        self.generations = generation + 1
 
         signal.signal(signal.SIGINT,signal.default_int_handler)
@@ -394 +415 @@
         self.scale = ScalingFactor
 
-        self.bestEnergy = self._init_popEnergy
+        #set initial solution and energy by running a single iteration
+        self.bestSolution = self.population[0][:]
+        self.bestEnergy = self.popEnergy[0]
+        trialPop = [[0.0 for i in range(self.nDim)] for j in range(self.nPop)]
+        for candidate in range(self.nPop):
+            # generate trialSolution (within valid range)
+            self.trialSolution[:] = self.population[candidate][:]
+            #XXX:[HACK] apply constraints
+            self.trialSolution[:] = constraints(self.trialSolution[:])
+            #XXX:[HACK] -end-
+            trialPop[candidate][:] = self.trialSolution[:]
+
+        trialEnergy = map(costfunction, trialPop)
+
+        for candidate in range(self.nPop):
+            if trialEnergy[candidate] < self.popEnergy[candidate]:
+                # New low for this candidate
+                self.popEnergy[candidate] = trialEnergy[candidate]
+                self.population[candidate][:] = trialPop[candidate][:]
+                self.UpdateGenealogyRecords(candidate, trialPop[candidate][:])
+               #XXX: was self.UpdateGenealogyRecords(candidate, self.trialSolution[:])
+
+                # Check if all-time low
+                if trialEnergy[candidate] < self.bestEnergy:
+                    self.bestEnergy = trialEnergy[candidate]
+                    self.bestSolution[:] = trialPop[candidate][:]
+                        
+        self.energy_history.append(self.bestEnergy)
+        StepMonitor(self.bestSolution[:], self.bestEnergy)
+        self.generations = 0  #XXX: above currently *not* counted as an iteration
+        if callback is not None:
+            callback(self.bestSolution)
          
         if self._maxiter is None:
@@ -400 +452 @@
         if self._maxfun is None:
             self._maxfun = self.nDim * self.nPop * 1000 #XXX: set better defaults?
-        trialPop = [[0.0 for i in range(self.nDim)] for j in range(self.nPop)]
-
-        generation = 0
-        for generation in range(self._maxiter):
-            StepMonitor(self.bestSolution[:], self.bestEnergy)
+
+        #run for generations <= maxiter
+        for generation in range(self._maxiter - self.generations):
             if fcalls[0] >= self._maxfun: break
             for candidate in range(self.nPop):
@@ -435 +485 @@
                             
             self.energy_history.append(self.bestEnergy)
-
+            StepMonitor(self.bestSolution[:], self.bestEnergy)
+            self.generations += 1
             if callback is not None:
                 callback(self.bestSolution)
@@ -441 +492 @@
             if self._EARLYEXIT or termination(self):
                 break
-
-        self.generations = generation + 1
 
         signal.signal(signal.SIGINT,signal.default_int_handler)

mystic/mystic/scipy_optimize.py

@@ -207 +207 @@
             allvecs = [sim[0]]
         fsim[0] = func(x0)
+        StepMonitor(sim[0], fsim[0]) # sim = all values; "best" is sim[0]
 
         #--- ensure initial simplex is within bounds ---
@@ -227 +228 @@
         self.popEnergy = fsim
         self.energy_history.append(self.bestEnergy)
+        StepMonitor(sim[0], fsim[0]) # sim = all values; "best" is sim[0]
 
         iterations = 1
 
         while (fcalls[0] < self._maxfun and iterations < self._maxiter):
-            StepMonitor(sim[0], fsim[0]) # sim = all values; "best" is sim[0]
             if self._EARLYEXIT or termination(self):
                 break
@@ -295 +296 @@
             self.popEnergy = fsim
             self.energy_history.append(self.bestEnergy)
+            StepMonitor(sim[0], fsim[0]) # sim = all values; "best" is sim[0]
 
         self.generations = iterations

mystic/mystic/tools.py

@@ -161 +161 @@
         from numpy import ndarray
         Sow.__call__(self, x, y)
-        self._step += 1
         if isinstance(y,(list,ndarray)):
             y = y[0] #XXX: get the "best" fit... which should be in y[0]
@@ -171 +170 @@
         if int(self._step % self._xinterval) == 0:
             print "Generation %d has best fit parameters:\n %s" % (self._step, x)
+        self._step += 1
         return
     pass