Changeset 171


Ignore:
Timestamp:
08/07/09 11:27:13 (7 years ago)
Author:
altafang
Message:

Adding uncertainty and fixing SetInitialPoints?. More revisions to come

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/alta/mystic-0.1a2/mystic/snobfit_solver.py

    r168 r171  
    2727    - StepMonitor = Sow() 
    2828    - enable_signal_handler() 
    29     - termination = SnobfitTermination(rtol, ftol, nstop) 
     29    - termination = SnobfitTermination(rtol, ftol, generations) 
    3030 
    3131Usage 
     
    7272# Import Mystic tools 
    7373from mystic.tools import Null, wrap_function 
    74 from mystic.tools import wrap_bounds, random_seed 
     74from mystic.tools import wrap_bounds 
    7575 
    7676from abstract_solver import AbstractSolver 
     
    123123    p -- probability of generating a evaluation point of Class 4. 
    124124    fac     -- Factor for multiplicative perturbation of the data. 
    125     disp   -- non-zero if fval and warnflag outputs are desired. 
    126     callback  -- an optional user-supplied function to call after each 
    127                  iteration. It's called as callback(xbest), where xbest 
    128                  is the current best point. 
    129125    fglob -- the user specified global function value. 
    130126 
     
    134130        self.p = 0.5 
    135131        self.fac = 0.0 
    136         disp = 0 
    137         callback = None 
    138132        self.fglob = 0 
    139133 
     
    141135        if kwds.has_key('p'): self.p = kwds['p'] 
    142136        if kwds.has_key('fac'): self.fac = kwds['fac'] 
    143         if kwds.has_key('disp'): disp = kwds['disp'] 
    144         if kwds.has_key('callback'): callback = kwds['callback'] 
    145137        if kwds.has_key('fglob'): self.fglob = kwds['fglob'] 
    146138 
     
    160152        #------------------------------------------------------------- 
    161153 
     154        self.func = func 
     155 
    162156        x0 = self.population[0] 
    163157        self.x = self.population 
     
    182176 
    183177        # Function call counter 
    184         ncall = self.nPop 
     178        #ncall = self.nPop   # <-- original snobfit's function call counter 
     179        initfcalls = fcalls[0] 
    185180 
    186181        # iteration counter 
     
    195190 
    196191        # Repeat until the limit on function calls is reached 
    197         while ncall < self._maxiter: 
     192        #while ncall < self._maxiter: 
     193        while fcalls[0] < self._maxiter: 
    198194             StepMonitor(self.bestSolution, self.bestEnergy) 
    199               
    200              if ncall == self.nPop: 
     195 
     196             #if ncall == self.nPop: 
     197             if fcalls[0] == initfcalls: 
    201198                 # initial call 
    202199                 self.fit(func, initialCall=True) 
     
    213210 
    214211             # update function call counter 
    215              ncall += self.nPop 
     212             #ncall += self.nPop 
    216213 
    217214             # best of the new function values 
     
    223220                 xopt = x[jbest,:] 
    224221 
    225              if callback is not None: 
    226                  callback(x[jbest]) 
     222             if self.callback is not None: 
     223                 self.callback(x[jbest]) 
    227224 
    228225             # Monitor progress 
     
    245242 
    246243        if fcalls[0] >= self._maxfun: 
    247             if disp: 
     244            if self.disp: 
    248245                print "Warning: Maximum number of function evaluations has "\ 
    249246                      "been exceeded." 
    250247        elif i >= self._maxiter: 
    251             if disp: 
     248            if self.disp: 
    252249                print "Warning: Maximum number of iterations has been exceeded" 
    253250        else: 
    254             if disp: 
     251            if self.disp: 
    255252                print "Optimization terminated successfully." 
    256253                print "         Current function value: %f" % fval 
     
    259256        return  
    260257 
    261     # Overwrite abstract_solver's SetInitialPoints() 
    262     def SetInitialPoints(self, x0): 
    263         """Set initial points with guess x0.""" 
    264         self.population = self._setInitRecommendedEvalPoints(x0) 
    265  
     258    # Overwrite abstract_solver's SetRandomInitialPoints() 
     259    def SetRandomInitialPoints(self, min=None, max=None): 
     260        """Set random initial points.""" 
     261        self.population = self._setInitRecommendedEvalPoints() 
     262 
     263  #  def SetInitialPoints(self, x0): 
     264  #      self.SetRandomInitialPoints() 
     265  #      self.population[0] = x0 
     266         
     267         
    266268################################################################################# 
    267269 
    268270# Helper methods for SnobfitSolver 
    269271 
    270     def _setInitRecommendedEvalPoints(self, x0): 
     272    def _setInitRecommendedEvalPoints(self): 
    271273        """ 
    272274        From the initial user-supplied guess, we set up the first set of 
     
    277279        Make sure the initial guess x0 is in recommended evaluation points. 
    278280        """ 
    279         x = numpy.dot( rand(self.nPop-1, self.nDim), 
     281        x = numpy.dot( rand(self.nPop, self.nDim), 
    280282                       numpy.diag(self._strictMax-self._strictMin) 
    281283                       ) + \ 
    282             numpy.array( [self._strictMin] ).repeat(self.nPop-1,axis=0) 
    283  
    284         x = numpy.append( [x0], x, axis=0 ) 
     284            numpy.array( [self._strictMin] ).repeat(self.nPop,axis=0) 
     285 
     286        #x = numpy.append( [x0], x, axis=0 ) 
    285287        return x 
    286288 
     
    21632165            self._warning() 
    21642166 
     2167########################################################################### 
     2168# Helper methods for SnobfitSolver to calculate the covariance  
     2169# matrix, uncertainty: 
     2170 
     2171    def _calcCovarAtSolution(self, x, f, df=None): 
     2172        """ calc the Correlation Matrix or Uncertainty at the solution 
     2173 
     2174        OutPut: 
     2175        ------- 
     2176        covar   the Correlation Matrix. 
     2177        """ 
     2178        if x  is not None: self.xnew  = x.copy() 
     2179        if f  is not None: self.fnew  = f.copy() 
     2180        if df is not None: 
     2181            self.dfnew = df.copy() 
     2182        else: 
     2183            self.dfnew = numpy.ones(len(f))*max(3*self.fac,numpy.sqrt(eps)) 
     2184 
     2185        self.updateWorkspace(self.func) 
     2186 
     2187        # Current best solution 
     2188        [fSolution, j] = min_( self.f ) 
     2189        xSolution = self.x[j,:] 
     2190 
     2191        K = min( self.nx-1, self.nDim*(self.nDim+3) ) 
     2192        distance=numpy.sum((self.x - \ 
     2193                           numpy.array([xSolution]).repeat(self.nx,axis=0))**2, 
     2194                           axis=1 
     2195                           ) 
     2196        ind = numpy.argsort(distance)[1:K+1] 
     2197 
     2198        # S^k = x^k - xSolution 
     2199        S = self.x[ind,:] - numpy.array([xSolution]).repeat(K,axis=0) 
     2200 
     2201        #---------- Fixed typo ------------------- 
     2202        #R = numpy.linalg.qr(a, mode='r') 
     2203        R = numpy.linalg.qr(S, mode='r') 
     2204        L = inv(R).T 
     2205 
     2206        # Scale matrix 
     2207        sc = numpy.sum( numpy.dot(S,L.T)**2, axis=1) ** (3.0/2.0 ) 
     2208        b = (self.f[ind]-fSolution)/sc 
     2209 
     2210        xx  = self.x[ind,:] - numpy.array([xSolution]).repeat(K,axis=0) 
     2211        xx2 = 0.5*xx*xx 
     2212        A   = numpy.bmat( 'xx xx2') 
     2213 
     2214        for i in xrange(self.nDim-1): 
     2215            xx =( self.x[ind,i] - xSolution[i]) 
     2216            B  = numpy.array( toCol(xx) ).repeat( self.nDim-i-1, axis=1) 
     2217            xx = B*(self.x[ind,i+1:self.nDim] - \ 
     2218                    numpy.array([ xSolution[i+1:self.nDim] ]).repeat(K,axis=0) ) 
     2219            A  = numpy.bmat('A xx') 
     2220 
     2221        xx = numpy.array( toCol(sc) ).repeat(A.shape[1], axis=1) 
     2222        A = A/xx 
     2223 
     2224        # In q, it includes the components of g and the upper triangle of G 
     2225        q = mldiv(A, toCol(b)) 
     2226 
     2227        # Jacobian matrix 
     2228        J = numpy.diag( q[self.nDim:2*self.nDim] )/2 
     2229        k = 2*self.nDim 
     2230        for i in xrange( self.nDim-1 ): 
     2231           for j in xrange(self.nDim-i-1): 
     2232               J[i, j+i+1] = q[k]/2 
     2233               J[j+i+1, i] = q[k]/2 
     2234               k += 1 
     2235 
     2236        # The covariance matrix is computed by, 
     2237        # covar = (J^T J)^{-1} 
     2238        covar = inv(J.T * J) 
     2239 
     2240        # If the minimisation uses the weighted least-squares function. 
     2241        # the covariance matrix should be multiplied by the variance of 
     2242        # the residuals about the best-solution least-squares. 
     2243        # For example, the fit of the reflectometry problem. 
     2244        if self.isLeastSqrt: 
     2245            # Calc the squared residuals(i.e., the goodness of fit) 
     2246            # at the solution. 
     2247            sumsq = abs(fSolution)**2/self.nDim 
     2248            covar *= sumsq 
     2249 
     2250 
     2251        return covar 
     2252 
     2253 
     2254    def _covarianceMatrix(self): 
     2255        """ 
     2256        Calculate the Covariance Matrix at the solution 
     2257        """ 
     2258        # Compute function values at the suggested points 
     2259        nrequest = self.request.shape[0] 
     2260        x = numpy.zeros( (nrequest,self.nDim) ) 
     2261        f = vector( nrequest ) 
     2262        for j in xrange( nrequest ): 
     2263            x[j,:] = self.request[j,0:self.nDim] 
     2264            f[j]   = self.func(x[j,:]) 
     2265 
     2266        return self. _calcCovarAtSolution(x, f) 
     2267 
     2268    # The alias 
     2269    covar = _covarianceMatrix 
     2270 
     2271 
     2272    def uncertainty(self, isLeastSqrt=False): 
     2273        """ 
     2274        Calculate the uncertainty of the fitting parameter at the solution. 
     2275 
     2276        Optional input: isLeastSqrt -- boolean -- True if the minimisation 
     2277                                    uses the least squares function.  
     2278        """ 
     2279        self.isLeastSqrt=isLeastSqrt 
     2280        covar = self._covarianceMatrix() 
     2281        print covar 
     2282        print numpy.diag(covar) 
     2283 
     2284        # calculate the uncertainty of fit parameters. 
     2285        uct = numpy.sqrt( numpy.diag(covar) ) 
     2286 
     2287        return uct 
     2288 
     2289 
     2290 
    21652291############################################################################### 
    2166 # If using an instance of SoftConstraints, use SnobfitSoftSolver instead of SnobfitSolver 
     2292# If using an instance of SoftConstraints (constraint is not None), 
     2293# use SnobfitSoftSolver instead of SnobfitSolver 
    21672294 
    21682295class SnobfitSoftSolver(SnobfitSolver): 
    2169     """Inherits from SnobfitSolver. Use if there are soft constraints.""" 
     2296    """Inherits from SnobfitSolver. Use this solver only if  
     2297there are soft constraints.""" 
    21702298 
    21712299    def Solve(self, func, termination, sigint_callback=None, 
     
    21912319    p -- probability of generating a evaluation point of Class 4. 
    21922320    fac     -- Factor for multiplicative perturbation of the data. 
    2193     disp   -- non-zero if fval and warnflag outputs are desired. 
    2194     callback  -- an optional user-supplied function to call after each 
    2195                  iteration. It's called as callback(xbest), where xbest 
    2196                  is the current best point. 
    21972321    fglob -- the user specified global function value. 
    21982322 
     
    22022326        self.p = 0.5 
    22032327        self.fac = 0.0 
    2204         disp = 0 
    2205         callback = None 
    22062328        self.fglob = 0 
    22072329 
     
    22092331        if kwds.has_key('p'): self.p = kwds['p'] 
    22102332        if kwds.has_key('fac'): self.fac = kwds['fac'] 
    2211         if kwds.has_key('disp'): disp = kwds['disp'] 
    2212         if kwds.has_key('callback'): callback = kwds['callback'] 
     2333        if kwds.has_key('fglob'): self.fglob = kwds['fglob'] 
     2334 
    22132335        # SoftConstraints requires a constraint passed to Solve().... 
    22142336        if kwds.has_key('constraint'): self.constraint = kwds['constraint'] 
    2215         if kwds.has_key('fglob'): self.fglob = kwds['fglob'] 
    22162337 
    22172338        #------------------------------------------------------------- 
     
    22532374        ncall = self.nPop 
    22542375        xopt  = inf 
    2255         improvement = 0 
    22562376        change = 0 
    22572377 
     
    23422462                 xopt = x[jbest,:] 
    23432463 
    2344              if  callback is not None : 
    2345                  callback(i, x[jbest], fbestn, improvement==0) 
     2464             if  self.callback is not None : 
     2465                 self.callback(x[jbest]) 
    23462466 
    23472467             if  fopt < 0 and change == 0: 
     
    23992519                print "         Function evaluations: %d" % fcalls[0] 
    24002520        return 
    2401  
    2402     # Overwrite abstractsolver's SetInitialPoints() 
    2403     def SetInitialPoints(self, x0): 
    2404         """Set initial points with guess x0.""" 
    2405         self.population = self._setInitRecommendedEvalPoints(x0) 
    24062521 
    24072522####################################################################### 
     
    24882603            retall = 0, 
    24892604            callback    = None, 
    2490             seed        = None, 
    24912605            constraint = None, 
    2492             full_output = 0): 
     2606            full_output = 0, 
     2607            retuct = False): 
    24932608    """Minimize a function using Snobfit - Stable Noisy Optimization by Branch and Fit. 
    24942609 
     
    25192634                 iteration. It's called as callback(xbest), where xbest 
    25202635                 is the current best point. 
    2521     seed     -- The seed for numpy.random. If the user set this number,it makes 
    2522                 sure that it returns the same solution from repeated tests 
    25232636    constraint -- an instance of SoftConstraints. If constraint is not None, then use  
    25242637                  SnobfitSoftSolver instead of SnobfitSolver. 
    25252638    full_output -- non-zero if fval and warnflag outputs are 
    25262639                  desired. 
    2527  
    2528  
    2529 Returns: (xopt, {fopt, iter, funcalls, warnflag}, {allvecs}) 
     2640    retuct  -- True to return the uncertainty of the fitting parameters 
     2641 
     2642 
     2643Returns: (xopt, {fopt, iter, funcalls, warnflag}, {allvecs}, {uncertainty}) 
    25302644 
    25312645    xopt -- ndarray - minimizer of function 
     
    25372651        2 : 'Maximum number of iterations.' 
    25382652    allvecs -- list - a list of solutions at each iteration 
     2653    uncertainty -- list - Returned if retuct is True. The uncertainty of  
     2654                   the fitting parameters. 
    25392655 
    25402656    """ 
     
    25452661    from mystic.termination import SnobfitTermination 
    25462662 
    2547     # Set the random seed 
    2548     if seed: 
    2549         random_seed(seed) 
    2550  
    2551     # Messy?: 
    2552     # Using npop instead of dn may not be the best idea? 
    25532663    if constraint is None: 
    25542664        if npop: 
     
    25622672            solver = SnobfitSoftSolver(len(x0)) 
    25632673 
     2674    solver.disp = disp 
     2675    solver.callback = callback 
     2676 
    25642677    # Problem?: SetStrictRanges must be called before SetInitialPoints 
    2565     # In fact, SetStrictRanges must always be called. Perhaps 
    2566     # integrate bounds into __init__? 
     2678    # In fact, _strictMin and _strictMax must always be set. Perhaps 
     2679    # integrate bounds into __init__? Make a default? 
    25672680    solver.SetStrictRanges(bounds[0],bounds[1]) 
    25682681    solver.SetInitialPoints(x0) 
     
    25722685                 EvaluationMonitor=evalmon,StepMonitor=stepmon,\ 
    25732686                 ExtraArgs=args,\ 
    2574                  p=p, fac=fac, disp=disp, fglob=fglob,\ 
    2575                  callback=callback, 
    2576                  constraint=constraint)  
     2687                 p=p, fac=fac, fglob=fglob,constraint=constraint)  
    25772688    solution = solver.Solution() 
    25782689 
    2579     # code below here pushes output to scipy.optimize interface 
    25802690    x = solver.bestSolution 
    25812691    fval = solver.bestEnergy 
     
    26012711            retlist = (x, allvecs) 
    26022712 
     2713    if retuct: 
     2714        retlist += (solver.uncertainty(),) 
     2715 
    26032716    return retlist 
    26042717 
Note: See TracChangeset for help on using the changeset viewer.