Changeset 171

Timestamp:

08/07/09 11:27:13 (7 years ago)

Author:

altafang

Message:

Adding uncertainty and fixing SetInitialPoints?. More revisions to come

File:

• branches/alta/mystic-0.1a2/mystic/snobfit_solver.py (modified) (27 diffs)

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

@@ -27 +27 @@
     - StepMonitor = Sow()
     - enable_signal_handler()
-    - termination = SnobfitTermination(rtol, ftol, nstop)
+    - termination = SnobfitTermination(rtol, ftol, generations)
 
 Usage
@@ -72 +72 @@
 # Import Mystic tools
 from mystic.tools import Null, wrap_function
-from mystic.tools import wrap_bounds, random_seed
+from mystic.tools import wrap_bounds
 
 from abstract_solver import AbstractSolver
@@ -123 +123 @@
     p -- probability of generating a evaluation point of Class 4.
     fac     -- Factor for multiplicative perturbation of the data.
-    disp   -- non-zero if fval and warnflag outputs are desired.
-    callback  -- an optional user-supplied function to call after each
-                 iteration. It's called as callback(xbest), where xbest
-                 is the current best point.
     fglob -- the user specified global function value.
 
@@ -134 +130 @@
         self.p = 0.5
         self.fac = 0.0
-        disp = 0
-        callback = None
         self.fglob = 0
 
@@ -141 +135 @@
         if kwds.has_key('p'): self.p = kwds['p']
         if kwds.has_key('fac'): self.fac = kwds['fac']
-        if kwds.has_key('disp'): disp = kwds['disp']
-        if kwds.has_key('callback'): callback = kwds['callback']
         if kwds.has_key('fglob'): self.fglob = kwds['fglob']
 
@@ -160 +152 @@
         #-------------------------------------------------------------
 
+        self.func = func
+
         x0 = self.population[0]
         self.x = self.population
@@ -182 +176 @@
 
         # Function call counter
-        ncall = self.nPop
+        #ncall = self.nPop   # <-- original snobfit's function call counter
+        initfcalls = fcalls[0]
 
         # iteration counter
@@ -195 +190 @@
 
         # Repeat until the limit on function calls is reached
-        while ncall < self._maxiter:
+        #while ncall < self._maxiter:
+        while fcalls[0] < self._maxiter:
              StepMonitor(self.bestSolution, self.bestEnergy)
-             
-             if ncall == self.nPop:
+
+             #if ncall == self.nPop:
+             if fcalls[0] == initfcalls:
                  # initial call
                  self.fit(func, initialCall=True)
@@ -213 +210 @@
 
              # update function call counter
-             ncall += self.nPop
+             #ncall += self.nPop
 
              # best of the new function values
@@ -223 +220 @@
                  xopt = x[jbest,:]
 
-             if callback is not None:
-                 callback(x[jbest])
+             if self.callback is not None:
+                 self.callback(x[jbest])
 
              # Monitor progress
@@ -245 +242 @@
 
         if fcalls[0] >= self._maxfun:
-            if disp:
+            if self.disp:
                 print "Warning: Maximum number of function evaluations has "\
                       "been exceeded."
         elif i >= self._maxiter:
-            if disp:
+            if self.disp:
                 print "Warning: Maximum number of iterations has been exceeded"
         else:
-            if disp:
+            if self.disp:
                 print "Optimization terminated successfully."
                 print "         Current function value: %f" % fval
@@ -259 +256 @@
         return 
 
-    # Overwrite abstract_solver's SetInitialPoints()
-    def SetInitialPoints(self, x0):
-        """Set initial points with guess x0."""
-        self.population = self._setInitRecommendedEvalPoints(x0)
-
+    # Overwrite abstract_solver's SetRandomInitialPoints()
+    def SetRandomInitialPoints(self, min=None, max=None):
+        """Set random initial points."""
+        self.population = self._setInitRecommendedEvalPoints()
+
+  #  def SetInitialPoints(self, x0):
+  #      self.SetRandomInitialPoints()
+  #      self.population[0] = x0
+        
+        
 #################################################################################
 
 # Helper methods for SnobfitSolver
 
-    def _setInitRecommendedEvalPoints(self, x0):
+    def _setInitRecommendedEvalPoints(self):
         """
         From the initial user-supplied guess, we set up the first set of
@@ -277 +279 @@
         Make sure the initial guess x0 is in recommended evaluation points.
         """
-        x = numpy.dot( rand(self.nPop-1, self.nDim),
+        x = numpy.dot( rand(self.nPop, self.nDim),
                        numpy.diag(self._strictMax-self._strictMin)
                        ) + \
-            numpy.array( [self._strictMin] ).repeat(self.nPop-1,axis=0)
-
-        x = numpy.append( [x0], x, axis=0 )
+            numpy.array( [self._strictMin] ).repeat(self.nPop,axis=0)
+
+        #x = numpy.append( [x0], x, axis=0 )
         return x
 
@@ -2163 +2165 @@
             self._warning()
 
+###########################################################################
+# Helper methods for SnobfitSolver to calculate the covariance 
+# matrix, uncertainty:
+
+    def _calcCovarAtSolution(self, x, f, df=None):
+        """ calc the Correlation Matrix or Uncertainty at the solution
+
+        OutPut:
+        -------
+        covar   the Correlation Matrix.
+        """
+        if x  is not None: self.xnew  = x.copy()
+        if f  is not None: self.fnew  = f.copy()
+        if df is not None:
+            self.dfnew = df.copy()
+        else:
+            self.dfnew = numpy.ones(len(f))*max(3*self.fac,numpy.sqrt(eps))
+
+        self.updateWorkspace(self.func)
+
+        # Current best solution
+        [fSolution, j] = min_( self.f )
+        xSolution = self.x[j,:]
+
+        K = min( self.nx-1, self.nDim*(self.nDim+3) )
+        distance=numpy.sum((self.x - \
+                           numpy.array([xSolution]).repeat(self.nx,axis=0))**2,
+                           axis=1
+                           )
+        ind = numpy.argsort(distance)[1:K+1]
+
+        # S^k = x^k - xSolution
+        S = self.x[ind,:] - numpy.array([xSolution]).repeat(K,axis=0)
+
+        #---------- Fixed typo -------------------
+        #R = numpy.linalg.qr(a, mode='r')
+        R = numpy.linalg.qr(S, mode='r')
+        L = inv(R).T
+
+        # Scale matrix
+        sc = numpy.sum( numpy.dot(S,L.T)**2, axis=1) ** (3.0/2.0 )
+        b = (self.f[ind]-fSolution)/sc
+
+        xx  = self.x[ind,:] - numpy.array([xSolution]).repeat(K,axis=0)
+        xx2 = 0.5*xx*xx
+        A   = numpy.bmat( 'xx xx2')
+
+        for i in xrange(self.nDim-1):
+            xx =( self.x[ind,i] - xSolution[i])
+            B  = numpy.array( toCol(xx) ).repeat( self.nDim-i-1, axis=1)
+            xx = B*(self.x[ind,i+1:self.nDim] - \
+                    numpy.array([ xSolution[i+1:self.nDim] ]).repeat(K,axis=0) )
+            A  = numpy.bmat('A xx')
+
+        xx = numpy.array( toCol(sc) ).repeat(A.shape[1], axis=1)
+        A = A/xx
+
+        # In q, it includes the components of g and the upper triangle of G
+        q = mldiv(A, toCol(b))
+
+        # Jacobian matrix
+        J = numpy.diag( q[self.nDim:2*self.nDim] )/2
+        k = 2*self.nDim
+        for i in xrange( self.nDim-1 ):
+           for j in xrange(self.nDim-i-1):
+               J[i, j+i+1] = q[k]/2
+               J[j+i+1, i] = q[k]/2
+               k += 1
+
+        # The covariance matrix is computed by,
+        # covar = (J^T J)^{-1}
+        covar = inv(J.T * J)
+
+        # If the minimisation uses the weighted least-squares function.
+        # the covariance matrix should be multiplied by the variance of
+        # the residuals about the best-solution least-squares.
+        # For example, the fit of the reflectometry problem.
+        if self.isLeastSqrt:
+            # Calc the squared residuals(i.e., the goodness of fit)
+            # at the solution.
+            sumsq = abs(fSolution)**2/self.nDim
+            covar *= sumsq
+
+
+        return covar
+
+
+    def _covarianceMatrix(self):
+        """
+        Calculate the Covariance Matrix at the solution
+        """
+        # Compute function values at the suggested points
+        nrequest = self.request.shape[0]
+        x = numpy.zeros( (nrequest,self.nDim) )
+        f = vector( nrequest )
+        for j in xrange( nrequest ):
+            x[j,:] = self.request[j,0:self.nDim]
+            f[j]   = self.func(x[j,:])
+
+        return self. _calcCovarAtSolution(x, f)
+
+    # The alias
+    covar = _covarianceMatrix
+
+
+    def uncertainty(self, isLeastSqrt=False):
+        """
+        Calculate the uncertainty of the fitting parameter at the solution.
+
+        Optional input: isLeastSqrt -- boolean -- True if the minimisation
+                                    uses the least squares function. 
+        """
+        self.isLeastSqrt=isLeastSqrt
+        covar = self._covarianceMatrix()
+        print covar
+        print numpy.diag(covar)
+
+        # calculate the uncertainty of fit parameters.
+        uct = numpy.sqrt( numpy.diag(covar) )
+
+        return uct
+
+
+
 ###############################################################################
-# If using an instance of SoftConstraints, use SnobfitSoftSolver instead of SnobfitSolver
+# If using an instance of SoftConstraints (constraint is not None),
+# use SnobfitSoftSolver instead of SnobfitSolver
 
 class SnobfitSoftSolver(SnobfitSolver):
-    """Inherits from SnobfitSolver. Use if there are soft constraints."""
+    """Inherits from SnobfitSolver. Use this solver only if 
+there are soft constraints."""
 
     def Solve(self, func, termination, sigint_callback=None,
@@ -2191 +2319 @@
     p -- probability of generating a evaluation point of Class 4.
     fac     -- Factor for multiplicative perturbation of the data.
-    disp   -- non-zero if fval and warnflag outputs are desired.
-    callback  -- an optional user-supplied function to call after each
-                 iteration. It's called as callback(xbest), where xbest
-                 is the current best point.
     fglob -- the user specified global function value.
 
@@ -2202 +2326 @@
         self.p = 0.5
         self.fac = 0.0
-        disp = 0
-        callback = None
         self.fglob = 0
 
@@ -2209 +2331 @@
         if kwds.has_key('p'): self.p = kwds['p']
         if kwds.has_key('fac'): self.fac = kwds['fac']
-        if kwds.has_key('disp'): disp = kwds['disp']
-        if kwds.has_key('callback'): callback = kwds['callback']
+        if kwds.has_key('fglob'): self.fglob = kwds['fglob']
+
         # SoftConstraints requires a constraint passed to Solve()....
         if kwds.has_key('constraint'): self.constraint = kwds['constraint']
-        if kwds.has_key('fglob'): self.fglob = kwds['fglob']
 
         #-------------------------------------------------------------
@@ -2253 +2374 @@
         ncall = self.nPop
         xopt  = inf
-        improvement = 0
         change = 0
 
@@ -2342 +2462 @@
                  xopt = x[jbest,:]
 
-             if  callback is not None :
-                 callback(i, x[jbest], fbestn, improvement==0)
+             if  self.callback is not None :
+                 self.callback(x[jbest])
 
              if  fopt < 0 and change == 0:
@@ -2399 +2519 @@
                 print "         Function evaluations: %d" % fcalls[0]
         return
-
-    # Overwrite abstractsolver's SetInitialPoints()
-    def SetInitialPoints(self, x0):
-        """Set initial points with guess x0."""
-        self.population = self._setInitRecommendedEvalPoints(x0)
 
 #######################################################################
@@ -2488 +2603 @@
             retall = 0,
             callback    = None,
-            seed        = None,
             constraint = None,
-            full_output = 0):
+            full_output = 0,
+            retuct = False):
     """Minimize a function using Snobfit - Stable Noisy Optimization by Branch and Fit.
 
@@ -2519 +2634 @@
                  iteration. It's called as callback(xbest), where xbest
                  is the current best point.
-    seed     -- The seed for numpy.random. If the user set this number,it makes
-                sure that it returns the same solution from repeated tests
     constraint -- an instance of SoftConstraints. If constraint is not None, then use 
                   SnobfitSoftSolver instead of SnobfitSolver.
     full_output -- non-zero if fval and warnflag outputs are
                   desired.
-
-
-Returns: (xopt, {fopt, iter, funcalls, warnflag}, {allvecs})
+    retuct  -- True to return the uncertainty of the fitting parameters
+
+
+Returns: (xopt, {fopt, iter, funcalls, warnflag}, {allvecs}, {uncertainty})
 
     xopt -- ndarray - minimizer of function
@@ -2537 +2651 @@
         2 : 'Maximum number of iterations.'
     allvecs -- list - a list of solutions at each iteration
+    uncertainty -- list - Returned if retuct is True. The uncertainty of 
+                   the fitting parameters.
 
     """
@@ -2545 +2661 @@
     from mystic.termination import SnobfitTermination
 
-    # Set the random seed
-    if seed:
-        random_seed(seed)
-
-    # Messy?:
-    # Using npop instead of dn may not be the best idea?
     if constraint is None:
         if npop:
@@ -2562 +2672 @@
             solver = SnobfitSoftSolver(len(x0))
 
+    solver.disp = disp
+    solver.callback = callback
+
     # Problem?: SetStrictRanges must be called before SetInitialPoints
-    # In fact, SetStrictRanges must always be called. Perhaps
-    # integrate bounds into __init__?
+    # In fact, _strictMin and _strictMax must always be set. Perhaps
+    # integrate bounds into __init__? Make a default?
     solver.SetStrictRanges(bounds[0],bounds[1])
     solver.SetInitialPoints(x0)
@@ -2572 +2685 @@
                  EvaluationMonitor=evalmon,StepMonitor=stepmon,\
                  ExtraArgs=args,\
-                 p=p, fac=fac, disp=disp, fglob=fglob,\
-                 callback=callback,
-                 constraint=constraint) 
+                 p=p, fac=fac, fglob=fglob,constraint=constraint) 
     solution = solver.Solution()
 
-    # code below here pushes output to scipy.optimize interface
     x = solver.bestSolution
     fval = solver.bestEnergy
@@ -2601 +2711 @@
             retlist = (x, allvecs)
 
+    if retuct:
+        retlist += (solver.uncertainty(),)
+
     return retlist