Changeset 176

Timestamp:

08/09/09 13:45:49 (7 years ago)

Author:

altafang

Message:

Fixing termination, gradient call counter

Location:

branches/alta/mystic-0.1a2/mystic

Files:

• scipy_bfgs.py (modified) (12 diffs)
• scipy_ncg.py (modified) (12 diffs)

branches/alta/mystic-0.1a2/mystic/scipy_bfgs.py

@@ -367 +367 @@
       epsilon : int or ndarray
           If fprime is approximated, use this value for the step size.
-      callback : callable
-          An optional user-supplied function to call after each
-          iteration.  Called as callback(xk), where xk is the
-          current parameter vector.
 
         """
@@ -376 +372 @@
         args = ExtraArgs
         x0 = self.population[0]
+
         fprime=None
         epsilon=_epsilon
-        disp=1
-        callback=None
+        self.disp = 0
+        self.callback = None
 
         if kwds.has_key('fprime'): fprime = kwds['fprime']
         if kwds.has_key('epsilon'): epsilon = kwds['epsilon']
-        if kwds.has_key('callback'): callback = kwds['callback']
-        if kwds.has_key('disp'): disp = kwds['disp']
+        if kwds.has_key('disp'): self.disp = kwds['disp']
+        if kwds.has_key('callback'): self.callback = kwds['callback']
+
         #-------------------------------------------------------------
 
@@ -407 +405 @@
             self._maxiter = len(x0)*200
 
-        # Wrap gradient monitor        
-        # ok to set monitor to Null here, or is there another way to do this?
+        # Wrap gradient monitor
         if fprime is None:
             gcalls, myfprime = wrap_function(approx_fprime, (func, epsilon), Null)
@@ -424 +421 @@
         xk = x0
         #sk = [2*gtol]
-        sk = [2*1e-5]
         warnflag = 0
         while k < self._maxiter:
@@ -455 +451 @@
             yk = gfkp1 - gfk
             gfk = gfkp1
-            if callback is not None:
-                callback(xk)
+            if self.callback is not None:
+                self.callback(xk)
             k += 1
 
@@ -465 +461 @@
             # Pass the gradient value to self.gfk for termination
             self.gfk = gfk
-    
+
+            # Calculate the inverse Hessian    
             try: # this was handled in numeric, let it remaines for more safety
                 rhok = 1.0 / (numpy.dot(yk,sk))
@@ -485 +482 @@
 
         self.generations = k 
-        signal.signal(signal.SIGINT,signal.default_int_handler)    
-
-        if disp:
+        signal.signal(signal.SIGINT,signal.default_int_handler) 
+
+        # Fix me?
+        self.gcalls = gcalls[0]
+
+        # Print messages if self.disp is non-zero
+        if self.disp:
             fval = old_fval
         if warnflag == 2:
-            if disp:
+            if self.disp:
                 print "Warning: Desired error not necessarily achieved" \
                       "due to precision loss"
@@ -499 +500 @@
     
         elif k >= self._maxiter:
-            if disp:
+            if self.disp:
                 print "Warning: Maximum number of iterations has been exceeded"
                 print "         Current function value: %f" % fval
@@ -506 +507 @@
                 print "         Gradient evaluations: %d" % gcalls[0]
         else:
-            if disp:
+            if self.disp:
                 print "Optimization terminated successfully."
                 print "         Current function value: %f" % fval
@@ -583 +584 @@
 
     solver = BFGSSolver(len(x0))
+
+    solver.disp = disp
+    solver.callback = callback
+
     solver.SetInitialPoints(x0)
     solver.enable_signal_handler()
@@ -588 +593 @@
     solver.Solve(func,termination=GradientTermination(gtol, norm),\
                  EvaluationMonitor=evalmon,StepMonitor=stepmon,\
-                 disp=disp, ExtraArgs=args, callback=callback,\
-                 epsilon=epsilon, fprime=fprime)
+                 ExtraArgs=args, epsilon=epsilon, fprime=fprime)
     solution = solver.Solution()
 
@@ -598 +602 @@
     warnflag = 0
     fcalls = len(evalmon.x)
+    iterations = len(stepmon.x)
+
+    # fix me?
     #gcalls = len(gradmon.x)
-    iterations = len(stepmon.x)
-    gcalls = iterations + 1 # fix me?
+    gcalls = solver.gcalls
+
     allvecs = []
     for i in range(iterations):

branches/alta/mystic-0.1a2/mystic/scipy_ncg.py

@@ -18 +18 @@
    - StepMonitor = CustomSow() with 4 columns
    - enable_signal_handler()
+   - termination = SolutionImprovement(tolerance)
 
 Usage
@@ -130 +131 @@
 
 
-    def Solve(self, func, fprime, termination, sigint_callback=None,
-              EvaluationMonitor=Null, StepMonitor=Null, GradMonitor=Null,
-              HessianMonitor=Null, ExtraArgs=(), **kwds):
+    def Solve(self, func, termination, sigint_callback=None,
+              EvaluationMonitor=Null, StepMonitor=Null, #GradientMonitor=Null,
+              ExtraArgs=(), **kwds):
         """Minimize a function using NCG.
 
@@ -143 +144 @@
 
     func -- the Python function or method to be minimized.
-    fprime -- callable f'(x,*args)
-            Gradient of f.
     termination -- callable object providing termination conditions.
 
@@ -158 +157 @@
 Further Inputs:
 
+    fprime -- callable f'(x,*args)
+            Gradient of f.
     fhess_p : callable fhess_p(x,p,*args)
         Function which computes the Hessian of f times an
@@ -186 +187 @@
         x0 = self.population[0]
         x0 = asarray(x0).flatten()
+
         epsilon = _epsilon
-        disp = 1
-        callback = None
+        self.disp = 1
+        self.callback = None
         fhess_p = None
         fhess = None
 
         if kwds.has_key('epsilon'): epsilon = kwds['epsilon']
-        if kwds.has_key('callback'): callback = kwds['callback']
-        if kwds.has_key('disp'): disp = kwds['disp']
+        if kwds.has_key('callback'): self.callback = kwds['callback']
+        if kwds.has_key('disp'): self.disp = kwds['disp']
         if kwds.has_key('fhess'): fhess = kwds['fhess']
         if kwds.has_key('fhess_p'): fhess_p = kwds['fhess_p']
+
+        # fprime is actually required. Temporary fix?:
+        if kwds.has_key('fprime'): fprime = kwds['fprime']
         #-------------------------------------------------------------
 
@@ -216 +221 @@
             self._maxiter = len(x0)*200
 
-        # Wrap gradient function?        
+        # Wrap gradient function? 
+        # gcalls, fprime = wrap_function(fprime, args, GradientMonitor)   
         gcalls, fprime = wrap_function(fprime, args, Null)
 
@@ -286 +292 @@
             update = alphak * pk
 
+            # Put last solution in trialSolution for termination()
             self.trialSolution = xk
 
             xk = xk + update        # upcast if necessary
-            if callback is not None:
-                callback(xk)
+            if self.callback is not None:
+                self.callback(xk)
             k += 1
 
@@ -304 +311 @@
         self.generations = k
 
+        # Fix me?
         self.hcalls = hcalls
+        self.gcalls = gcalls[0]
 
         signal.signal(signal.SIGINT,signal.default_int_handler)
     
-        if disp:
+        if self.disp:
             fval = old_fval
         if k >= self._maxiter:
-            if disp:
+            if self.disp:
                 print "Warning: Maximum number of iterations has been exceeded"
                 print "         Current function value: %f" % fval
@@ -319 +328 @@
                 print "         Hessian evaluations: %d" % hcalls
         else:
-            if disp:
+            if self.disp:
                 print "Optimization terminated successfully."
                 print "         Current function value: %f" % fval
@@ -402 +411 @@
  
     from mystic.tools import Sow, CustomSow
-    from mystic.termination import XTermination
+    from mystic.termination import SolutionImprovement
     #stepmon = Sow()
     stepmon = CustomSow('x','y','g','h', x='x', y='fval', \
@@ -413 +422 @@
     solver.SetEvaluationLimits(maxiter,None)
     # Does requiring fprime break abstract_solver interface?
-    solver.Solve(func, fprime, termination=XTermination(xtol),\
+    solver.Solve(func, SolutionImprovement(tolerance=xtol),\
                  EvaluationMonitor=evalmon,StepMonitor=stepmon,\
                  disp=disp, ExtraArgs=args, callback=callback,\
                  epsilon=epsilon, fhess_p=fhess_p,\
-                 fhess=fhess)
+                 fhess=fhess, fprime=fprime)
     solution = solver.Solution()
 
@@ -427 +436 @@
     fcalls = len(evalmon.x)
     iterations = len(stepmon.x)
-    gcalls = iterations   # Fix me?
+
+    # Fix me?
+    gcalls = solver.gcalls  
     hcalls = solver.hcalls
+
     allvecs = []
     for i in range(iterations):