Changeset 175


Ignore:
Timestamp:
08/09/09 13:04:00 (7 years ago)
Author:
altafang
Message:

Editing termination conditions and making various edits to annealing and snobfit

Location:
branches/alta/mystic-0.1a2/mystic
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/alta/mystic-0.1a2/mystic/anneal_solver.py

    r160 r175  
    2020    - StepMonitor = Sow() 
    2121    - enable_signal_handler() 
    22     - termination = AnnealTermination() 
     22    - termination = AveragePopulationImprovement(tolerance) 
    2323 
    2424Usage 
    … …  
    189189        return 
    190190 
    191 # Replace uses of this with Sow()? 
    192191class _state(object): 
    193192    def __init__(self): 
    … …  
    238237    T0           -- Initial Temperature (estimated as 1.2 times the largest 
    239238                    cost-function deviation over random points in the range) 
    240                     [default=None] 
    241     Tf           -- Final goal temperature 
    242                     [default=1e-12] 
    243     maxaccept    -- Maximum changes to accept 
    244239                    [default=None] 
    245240    learn_rate   -- scale constant for adjusting guesses 
    … …  
    252247    dwell        -- The number of times to search the space at each temperature. 
    253248                    [default=50] 
    254     disp         -- Non-zero to print convergence messages. 
    255                     [default=0] 
     249 
     250    Optional Termination Conditions: 
     251    Tf           -- Final goal temperature 
     252                    [default=1e-12] 
     253    maxaccept    -- Maximum changes to accept 
     254                    [default=None] 
    256255        """ 
    257256        #allow for inputs that don't conform to AbstractSolver interface 
    … …  
    259258        x0 = self.population[0] 
    260259 
    261         callback = None 
    262260        schedule = "fast" 
    263261        T0 = None 
    264         Tf = 1e-12 
    265262        boltzmann = 1.0 
    266263        learn_rate = 0.5 
    … …  
    269266        m = 1.0 
    270267        n = 1.0 
    271         disp = 0 
    272         self._maxaccept = None    
    273  
    274         if kwds.has_key('callback'): callback = kwds['callback'] 
     268 
     269        Tf = 1e-12 # or None? 
     270        self._maxaccept = None 
     271 
     272        self.disp = 0 
     273        self.callback = None 
     274 
    275275        if kwds.has_key('schedule'): schedule = kwds['schedule'] 
    276276        if kwds.has_key('T0'): T0 = kwds['T0'] 
    277         if kwds.has_key('Tf'): Tf = kwds['Tf'] 
    278277        if kwds.has_key('boltzmann'): boltzmann = kwds['boltzmann'] 
    279278        if kwds.has_key('learn_rate'): learn_rate = kwds['learn_rate'] 
    … …  
    282281        if kwds.has_key('m'): m = kwds['m'] 
    283282        if kwds.has_key('n'): n = kwds['n']     
    284         if kwds.has_key('disp'): disp = kwds['disp'] 
     283 
     284        if kwds.has_key('Tf'): Tf = kwds['Tf'] 
    285285        if kwds.has_key('maxaccept'): self._maxaccept = kwds['maxaccept'] 
     286 
     287        if kwds.has_key('disp'): self.disp = kwds['disp'] 
     288        if kwds.has_key('callback'): self.callback = kwds['callback'] 
     289 
    286290        #------------------------------------------------------------- 
    287291 
    … …  
    293297        if self._useStrictRange: 
    294298            x0 = self._clipGuessWithinRangeBoundary(x0) 
     299            # Note: wrap_bounds changes the results slightly from the original 
    295300            func = wrap_bounds(func, self._strictMin, self._strictMax) 
    296301 
    … …  
    314319            self.bestSolution = x0 
    315320        else: 
    316             self.bestSolution = None 
     321            #self.bestSolution = None 
     322            self.bestSolution = x0 
    317323            self.bestEnergy = 300e8 
    318324 
    … …  
    327333        schedule.T = schedule.T0 
    328334        fqueue = [100, 300, 500, 700] 
     335        self.population = asarray(fqueue)*1.0 
    329336        iters = 0 
    330337        while 1: 
     338            StepMonitor(self.bestSolution, self.bestEnergy) 
    331339            for n in range(dwell): 
    332340                current_state.x = schedule.update_guess(last_state.x) 
    … …  
    348356            af = asarray(fqueue)*1.0 
    349357 
     358            # Update monitors/variables 
    350359            self.population = af 
    351360            self.energy_history.append(self.bestEnergy) 
    352             StepMonitor(self.bestSolution[:], self.bestEnergy) 
    353  
    354             if callback is not None: 
    355                 callback(self.bestSolution) 
     361 
     362            if self.callback is not None: 
     363                self.callback(self.bestSolution) 
    356364 
    357365            # Stopping conditions 
    … …  
    373381            if (Tf is not None) and (schedule.T < Tf): 
    374382                break 
    375             if (self._maxfun is not None) and (schedule.feval > self._maxfun): 
     383           # if (self._maxfun is not None) and (schedule.feval > self._maxfun): 
     384           #     retval = 1 
     385           #     break 
     386            if (self._maxfun is not None) and (fcalls[0] > self._maxfun): 
    376387                retval = 1 
    377388                break 
    … …  
    382393                break 
    383394 
     395        signal.signal(signal.SIGINT,signal.default_int_handler) 
     396 
    384397        # Store some information. Is there a better way to do this? 
    385398        self.generations = iters        # Number of iterations 
    … …  
    387400        self.accept = schedule.accepted # Number of tests accepted 
    388401 
    389         signal.signal(signal.SIGINT,signal.default_int_handler) 
    390  
    391402        # code below here pushes output to scipy.optimize interface 
    392403 
    393         if disp: 
     404        if self.disp: 
    394405            if retval == 1:  
    395406                print "Warning: Maximum number of function evaluations has "\ 
    … …  
    452463    retval  -- Flag indicating stopping condition: 
    453464                1 : Maximum function evaluations 
    454                 2 : Maximum cooling iterations reached 
     465                2 : Maximum iterations reached 
    455466                3 : Maximum accepted query locations reached 
    456467    allvecs -- a list of solutions at each iteration. 
    … …  
    458469 
    459470    from mystic.tools import Sow 
     471    from mystic.termination import AveragePopulationImprovement 
    460472    stepmon = Sow() 
    461473    evalmon = Sow() 
    462     from mystic.termination import AnnealTermination 
    463  
     474     
    464475    x0 = asarray(x0) 
    465  
    466476    solver = AnnealSolver(len(x0)) 
    467477    solver.SetInitialPoints(x0) 
    … …  
    470480     
    471481    # Default bounds cause errors when function is not 1-dimensional 
    472     lower = asarray(lower) 
    473     upper = asarray(upper) 
    474     solver.SetStrictRanges(lower,upper)  
    475     solver.Solve(func,termination=AnnealTermination(feps),\ 
     482    solver.SetStrictRanges(lower,upper) 
     483 
     484    solver.Solve(func,termination=AveragePopulationImprovement(tolerance=feps),\ 
    476485                 EvaluationMonitor=evalmon,StepMonitor=stepmon,\ 
    477486                 ExtraArgs=args, callback=callback,\ 
    … …  
    483492 
    484493    # code below here pushes output to scipy.optimize interface 
    485     x = list(solver.bestSolution) 
    486     #x = solver.bestSolution 
     494    #x = list(solver.bestSolution) 
     495    x = solver.bestSolution 
    487496    fval = solver.bestEnergy 
    488497    fcalls = len(evalmon.x) 
    … …  
    493502        allvecs.append(stepmon.x[i]) 
    494503 
    495     # to be fixed? 
     504    # to be fixed 
    496505    T = solver.T 
    497506    accept = solver.accept 
    … …  
    516525    return retlist 
    517526 
     527if __name__=='__main__': 
     528    help(__name__) 
     529 
     530# end of file 

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

    r174 r175  
    2727    - StepMonitor = Sow() 
    2828    - enable_signal_handler() 
    29     - termination = SnobfitTermination(rtol, ftol, generations) 
     29    - termination = FglobTermination(fglob, tolerance, generations) for fglob !=0. 
     30                  = VTRChangeOverGenerations(costtolerance,gentolerance,generations)  
     31                    for fglob == 0 
     32 
     33                where fglob is the user-specified global function value (default=0). 
    3034 
    3135Usage 
    … …  
    123127    p -- probability of generating a evaluation point of Class 4. 
    124128    fac     -- Factor for multiplicative perturbation of the data. 
    125     fglob -- the user specified global function value. 
    126129 
    127130        """ 
    … …  
    130133        self.p = 0.5 
    131134        self.fac = 0.0 
    132         self.fglob = 0 
    133         self.disp = False 
    134135        self.callback = None 
     136        self.disp = 0 
    135137 
    136138        # Initialize some variables depending on input 
    137139        if kwds.has_key('p'): self.p = kwds['p'] 
    138140        if kwds.has_key('fac'): self.fac = kwds['fac'] 
    139         if kwds.has_key('fglob'): self.fglob = kwds['fglob'] 
     141        if kwds.has_key('callback'): self.callback = kwds['callback'] 
    140142        if kwds.has_key('disp'): self.disp = kwds['disp'] 
    141         if kwds.has_key('callback'): self.callback = kwds['callback'] 
     143 
     144        # For convenience (and uncertainty calculations): 
     145        self.func = func 
    142146 
    143147        #------------------------------------------------------------- 
    … …  
    145149        self._EARLYEXIT = False 
    146150 
    147         fcalls, func = wrap_function(func, ExtraArgs, EvaluationMonitor) 
     151        fcalls, self.func = wrap_function(self.func, ExtraArgs, EvaluationMonitor) 
    148152        if self._useStrictRange: 
    149             #for i in range(self.nPop): 
    150             #    self.population[i] = self._clipGuessWithinRangeBoundary(self.population[i]) 
    151             #func = wrap_bounds(func, self._strictMin, self._strictMax) 
    152  
    153             # Fix for when bounds were set after initial points.  
    154             # Snobfit's initialization uses the bounds not just for clipping, so  
    155             # the above code wouldn't really be quite right. 
    156             # Problem with this: what if SetRandomInitialPoints was called and  
    157             # the self.population[0] was something weird? 
     153            for i in range(self.nPop): 
     154                self.population[i] = self._clipGuessWithinRangeBoundary(self.population[i]) 
     155            self.func = wrap_bounds(self.func, self._strictMin, self._strictMax) 
     156 
     157            # Fix for when bounds were set after initial points: reinitialize 
     158            # This may give different results from if SetStrictRanges was called first 
    158159            x0 = self.population[0] 
    159160            self.SetInitialPoints(x0) 
    … …  
    164165        #------------------------------------------------------------- 
    165166 
    166         # to be fixed 
    167         self.func = func 
    168  
    169167        x0 = self.population[0] 
    170168        self.x = self.population 
    171169 
     170        self.fglob = 0 
     171 
    172172        # The number of safeguarded nearest neighbors 
    173173        self.snn = self.nPop - 1 
    … …  
    178178        # Other initializations 
    179179        self.constraint = None 
    180         (self.f, self.df) = self._setInitial_fdf(func) 
     180        (self.f, self.df) = self._setInitial_fdf() 
    181181        self.nx = self.x.shape[0] 
    182182        self.setInit()    
    … …  
    210210             if fcalls[0] == initfcalls: 
    211211                 # initial call 
    212                  self.fit(func, initialCall=True) 
     212                 self.fit(initialCall=True) 
    213213 
    214214             else: 
    215215                 # continuation call 
    216                  self.fit(func, x=x, f=f) 
     216                 self.fit(x=x, f=f) 
    217217 
    218218             [xopt, fopt] = (self.xbest, self.fbest) 
    … …  
    220220             for j in xrange( self.nPop ): 
    221221                 x[j,:] = self.request[j,0:self.nDim] 
    222                  f[j]   = func(x[j,:]) 
     222                 f[j]   = self.func(x[j,:]) 
    223223 
    224224             # update function call counter 
    225              #ncall += self.nPop    # <-- original snobfit function call counter 
     225             #ncall += self.nPop    
    226226 
    227227             # best of the new function values 
    … …  
    232232                 fopt = fbestn 
    233233                 xopt = x[jbest,:] 
    234  
    235              if self.callback is not None: 
    236                  self.callback(x[jbest]) 
    237234 
    238235             # Monitor progress 
    … …  
    242239             self.popEnergy = self.f 
    243240             self.energy_history.append(self.bestEnergy) 
     241 
     242             if self.callback is not None: 
     243                 #self.callback(x[jbest]) 
     244                 self.callback(self.bestSolution) 
    244245 
    245246             i += 1 
    … …  
    270271 
    271272    # Overwrite abstract_solver's SetRandomInitialPoints() 
    272  
     273     
    273274    # If SetStrictRange was called first, use those bounds. 
    274275    # If SetInitialPoints was called first, temporarily use defaults, 
    … …  
    282283            self._useStrictRange = False 
    283284        else: 
    284             # Defaults copied from abstract_solver. They do not make a difference 
    285             # because Solve() will re-initialize with new min and max 
     285            # Defaults copied from abstract_solver.  
     286            # Solve() will re-initialize with new min and max 
    286287            if min is None: min = [-1e3]*self.nDim   
    287288            if max is None: max = [1e3]*self.nDim 
    288  
    289         self.population = self._setInitRecommendedEvalPoints(min,max) 
     289        self.population = self._setInitRecommendedEvalPoints(min,max)  
     290        print self.population 
     291   # def SetInitialPoints(self,x0): 
     292   #     self.population = self._setInitRecommendedEvalPoints(self._strictMin,self._strictMax) 
     293   #     self.population[0] = x0 
    290294         
    291 ################################################################################# 
    292  
     295############################################################################ 
    293296# Helper methods for SnobfitSolver 
    294297 
    295     def _setInitRecommendedEvalPoints(self, min, max): 
     298    def _setInitRecommendedEvalPoints(self,min,max): 
    296299        """ 
    297300        Snobfit's initialization of the population 
    … …  
    301304                       ) + \ 
    302305            numpy.array( [min] ).repeat(self.nPop,axis=0) 
    303  
    304         #x = numpy.append( [x0], x, axis=0 ) 
     306       # x = numpy.append( [x0], x, axis=0 ) 
    305307        return x 
    306308 
    307     def _setInitial_fdf(self, func): 
     309    def _setInitial_fdf(self): 
    308310        """ 
    309311        The computation of the function values of the first set of the 
    … …  
    319321 
    320322        for i in xrange(self.nPop): 
    321             f[i]  = func( self.x[i,:] ) 
     323            f[i]  = self.func( self.x[i,:] ) 
    322324            df[i] = max(3*self.fac, numpy.sqrt(eps)) 
    323325 
    … …  
    683685 
    684686 
    685     def _calc_SplitPoint(self, j, x1, x2, func): 
     687    def _calc_SplitPoint(self, j, x1, x2): 
    686688        """ 
    687689        Compute the split point 
    688690        """ 
    689         _lambda = self._calc_Lambda(func(self.x[j,:]), 
    690                                     func(self.x[j+1,:]) 
     691        _lambda = self._calc_Lambda(self.func(self.x[j,:]), 
     692                                    self.func(self.x[j+1,:]) 
    691693                                    ) 
    692694        return _lambda*x1 + (1.0 - _lambda)*x2 
    693695 
    694696 
    695     def _split(self, x, f, df, xl0, xu0, nspl, func): 
     697    def _split(self, x, f, df, xl0, xu0, nspl): 
    696698        """ Split a box 
    697699         
    … …  
    764766        j = numpy.argmax( y[1:len(y)]-y[0:len(y)-1] ) 
    765767        if self.constraint is None: 
    766            ymid = self._calc_SplitPoint(j, y[j], y[j+1], func) 
     768           ymid = self._calc_SplitPoint(j, y[j], y[j+1]) 
    767769        else: 
    768770           ymid = 0.5*(y[j]+y[j+1]) 
    … …  
    804806            k = numpy.argmax( d ) 
    805807            if self.constraint is None: 
    806                ymid = self._calc_SplitPoint(k, y[k], y[k+1], func) 
     808               ymid = self._calc_SplitPoint(k, y[k], y[k+1]) 
    807809            else: 
    808810               ymid = 0.5*(y[k] + y[k+1])  
    … …  
    984986                             fnew, # new function values and 
    985987                             dfnew, # their variations 
    986                              func 
    987988                             ): 
    988989        """ 
    … …  
    11301131                                       xl[j,:], 
    11311132                                       xu[j,:], 
    1132                                        self.nsplit[j,:], 
    1133                                        func 
     1133                                       self.nsplit[j,:] 
    11341134                                       ) 
    11351135            k = self._sumFind( x0 == \ 
    … …  
    13411341    # Step 3 
    13421342    #------------------------------------ 
    1343     def branch(self, func): 
     1343    def branch(self): 
    13441344        """ 
    13451345        All current boxes containing more than one point are split 
    … …  
    13571357                                  self.u, 
    13581358                                  self.v, 
    1359                                   numpy.zeros( self.nDim), 
    1360                                   func 
     1359                                  numpy.zeros( self.nDim) 
    13611360                                  ) 
    13621361        else: 
    … …  
    15381537 
    15391538 
    1540     def updateWorkspace(self, func): 
     1539    def updateWorkspace(self): 
    15411540        """ update the current work space """ 
    15421541        oldxbest = self.xbest 
    … …  
    15471546                                           self.xnew, 
    15481547                                           self.fnew, 
    1549                                            self.dfnew, 
    1550                                            func 
     1548                                           self.dfnew 
    15511549                                           ) 
    15521550 
    … …  
    21102108 
    21112109 
    2112     def _q(self, f, f0, Delta): 
    2113         return ( f - f0 ) / ( Delta + abs(f - f0) ) 
    2114  
    2115      
    2116  
    21172110    # ----------------------------------------------------------------- 
    21182111 
    2119     def fit(self, func, x=None, f=None, df=None, initialCall=False): 
     2112    def fit(self, x=None, f=None, df=None, initialCall=False): 
    21202113        """ 
    21212114        One call of snobfit 
    … …  
    21412134            self.update_uv()     # Step 1 
    21422135            self.update_input()  # Step 2 
    2143             self.branch(func)        # step 3 
     2136            self.branch()        # step 3 
    21442137 
    21452138            (fmax, fmin) = self._getExtrmeFunc() 
    … …  
    21622155               self.dfnew = numpy.ones(len(f))*max(3*self.fac,numpy.sqrt(eps)) 
    21632156 
    2164             self.updateWorkspace(func) 
     2157            self.updateWorkspace() 
    21652158 
    21662159            if  isEmpty(self.near): # This is equal go to Step 11 
    … …  
    21872180# matrix, uncertainty: 
    21882181# Note: uncertainty() sometimes returns NaN or causes errors.... 
     2182 
     2183    def _q(self, f, f0, Delta): 
     2184        return ( f - f0 ) / ( Delta + abs(f - f0) ) 
    21892185 
    21902186    def _calcCovarAtSolution(self, x, f, df=None): 
    … …  
    22022198            self.dfnew = numpy.ones(len(f))*max(3*self.fac,numpy.sqrt(eps)) 
    22032199 
    2204         self.updateWorkspace(self.func) 
     2200        self.updateWorkspace() 
    22052201 
    22062202        # Current best solution 
    … …  
    23042300        return uct 
    23052301 
    2306  
    2307  
    23082302############################################################################### 
    23092303# If using an instance of SoftConstraints (constraint is not None), 
    … …  
    23432337        self.p = 0.5 
    23442338        self.fac = 0.0 
     2339        self.disp = 0 
     2340        self.callback = None 
     2341 
     2342        # Fix me? SnobfitSoftSolver requires fglob to be entered in both  
     2343        # termination and Solve().... 
    23452344        self.fglob = 0 
    2346         self.disp = False 
    2347         self.callback = None 
    23482345 
    23492346        # Initialize some variables depending on input 
    … …  
    23512348        if kwds.has_key('fac'): self.fac = kwds['fac'] 
    23522349        if kwds.has_key('fglob'): self.fglob = kwds['fglob'] 
     2350        if kwds.has_key('callback'): self.callback = kwds['callback'] 
    23532351        if kwds.has_key('disp'): self.disp = kwds['disp'] 
    2354         if kwds.has_key('callback'): self.callback = kwds['callback'] 
    23552352 
    23562353        # SoftConstraints requires a constraint passed to Solve().... 
    23572354        if kwds.has_key('constraint'): self.constraint = kwds['constraint'] 
     2355 
     2356        # Make func an instance variable 
     2357        self.func = func 
    23582358 
    23592359        #------------------------------------------------------------- 
    … …  
    23612361        self._EARLYEXIT = False 
    23622362 
    2363         fcalls, func = wrap_function(func, ExtraArgs, EvaluationMonitor) 
     2363        fcalls, self.func = wrap_function(self.func, ExtraArgs, EvaluationMonitor) 
    23642364        if self._useStrictRange: 
    2365             #for i in range(self.nPop): 
    2366             #    self.population[i] = self._clipGuessWithinRangeBoundary(self.population[i]) 
    2367             #func = wrap_bounds(func, self._strictMin, self._strictMax) 
     2365            for i in range(self.nPop): 
     2366                self.population[i] = self._clipGuessWithinRangeBoundary(self.population[i]) 
     2367            self.func = wrap_bounds(self.func, self._strictMin, self._strictMax) 
    23682368 
    23692369            # Fix for when bounds were set after initial points.  
    2370             # Snobfit's initialization uses the bounds not just for clipping, so  
    2371             # the above code wouldn't really be quite right. 
    2372             # Problem with this: what if SetRandomInitialPoints was called and  
    2373             # the self.population[0] was something weird? 
    23742370            x0 = self.population[0] 
    23752371            self.SetInitialPoints(x0) 
    … …  
    23902386        # Other initialization 
    23912387        self.x = self.population 
    2392         (self.f, self.df) = self._setInitial_fdf(func) 
     2388        (self.f, self.df) = self._setInitial_fdf() 
    23932389        self.nx = self.x.shape[0] 
    23942390        self.setInit()    
    … …  
    24022398        # Function call counter 
    24032399        #ncall = self.nPop 
    2404         initfcalls = fcalls[0] 
    24052400 
    24062401        xopt  = inf 
    … …  
    24142409        _f = vector( self.nPop ) 
    24152410        for j in xrange( self.nPop ): 
    2416           _f[j] = func( _x[j] ) 
     2411          _f[j] = self.func( _x[j] ) 
    24172412 
    24182413 
    … …  
    24362431        self.bestEnergy = self._init_popEnergy 
    24372432 
     2433        initfcalls = fcalls[0] 
     2434 
    24382435        # Repeat until the limit on function calls is reached 
    24392436        #while ncall < self._maxiter: 
    … …  
    24452442             if fcalls[0] == initfcalls: 
    24462443                 # initial call 
    2447                  self.fit(func, initialCall=True) 
     2444                 self.fit(initialCall=True) 
    24482445 
    24492446             else: 
    24502447                 # continuation call 
    2451                  self.fit(func, x=x, f=fm) 
     2448                 self.fit(x=x, f=fm) 
    24522449 
    24532450             [xopt, fopt] = (self.xbest, self.fbest) 
    … …  
    24592456             for j in xrange( self.nPop ): 
    24602457                 x[j,:] = self.request[j,0:self.nDim] 
    2461                  f[j]   = func(x[j,:]) 
     2458                 f[j]   = self.func(x[j,:]) 
    24622459                 FF    = self.constraint.F(x[j,:]) 
    24632460                 F[j]  = FF 
    … …  
    24962493                 fopt = fbestn 
    24972494                 xopt = x[jbest,:] 
    2498  
    2499              if  self.callback is not None : 
    2500                  self.callback(x[jbest]) 
    25012495 
    25022496             if  fopt < 0 and change == 0: 
    … …  
    25172511                 x = _x.copy() 
    25182512 
    2519              _fopt = func(xopt) 
     2513             _fopt = self.func(xopt) 
    25202514 
    25212515             # Monitor progress 
    … …  
    25262520             self.energy_history.append(self.bestEnergy) 
    25272521 
     2522             if  self.callback is not None : 
     2523                 #self.callback(x[jbest]) 
     2524                 self.callback(self.bestSolution) 
     2525 
    25282526             i+=1 
    25292527 
    … …  
    25312529                 break 
    25322530 
    2533         import signal 
    25342531        self.generations = i 
    25352532 
    … …  
    25382535        # code below here pushes output to scipy.optimize interface 
    25392536        fval = self.bestEnergy 
    2540         warnflag = 0 
    25412537 
    25422538        if fcalls[0] >= self._maxfun: 
    … …  
    26882684    uncertainty -- list - Returned if retuct is True. The uncertainty of  
    26892685                   the fitting parameters. 
    2690  
    26912686    """ 
    26922687 
    26932688    from mystic.tools import Sow 
     2689    from mystic.termination import FglobTermination, VTRChangeOverGenerations 
     2690 
    26942691    stepmon = Sow() 
    26952692    evalmon = Sow() 
    2696     from mystic.termination import SnobfitTermination 
    26972693 
    26982694    if constraint is None: 
    … …  
    27072703            solver = SnobfitSoftSolver(len(x0)) 
    27082704 
     2705    if fglob != 0: 
     2706        termination = FglobTermination(fglob=fglob, tolerance=rtol, generations=nstop) 
     2707    else: 
     2708        termination = VTRChangeOverGenerations(costtolerance=ftol, gentolerance=0,\ 
     2709                                               generations=nstop) 
     2710 
    27092711    # Bounds are required 
    27102712    solver.SetStrictRanges(bounds[0],bounds[1]) 
    … …  
    27122714    solver.enable_signal_handler() 
    27132715    solver.SetEvaluationLimits(maxiter,maxfun) 
    2714     solver.Solve(func,termination=SnobfitTermination(rtol, ftol, nstop),\ 
     2716    solver.Solve(func,termination,\ 
    27152717                 EvaluationMonitor=evalmon,StepMonitor=stepmon,\ 
    27162718                 ExtraArgs=args,disp=disp, callback=callback,\ 
    … …  
    27212723    fval = solver.bestEnergy 
    27222724    warnflag = 0 
    2723     fcalls = len(evalmon.x)  # Different from ncalls of original snobfit 
     2725    fcalls = len(evalmon.x) 
    27242726    iterations = len(stepmon.x) 
    27252727    allvecs = [] 

  • branches/alta/mystic-0.1a2/mystic/termination.py

    r169 r175  
    6060    return _ 
    6161 
    62 def SnobfitTermination(rtol=1e-6, ftol=1e-6, generations=250): 
     62def FglobTermination(fglob=0, tolerance=1e-6, generations = None): 
    6363    """ 
    64 Parameters: 
    65 rtol -- a relative error for checking stopping criterion. 
    66 ftol -- acceptable relative error in func(xopt) for convergence. 
    67 generations  -- number of times no improvement is tolerated. 
     64Input Parameters: 
     65fglob     -- User-specified global function value 
     66tolerance -- a relative error for checking stopping criterion. 
     67generations -- number of times no improvement is tolerated. 
    6868 
    69 Terminate if no improvement after iterations = generations or  if  
    70 abs((fopt-fglob)/fglob) < rtol for fglob != 0, abs(fopt) < ftol  
    71 for fglob == 0""" 
     69Terminate if abs((bestEnergy-fglob)/fglob < tolerance)""" 
    7270 
    7371    def _(inst): 
    74          answer = False 
    75          hist = inst.energy_history 
    76          if len(hist) > generations and (hist[-generations]-hist[-1]) < 0: 
    77              answer = True 
    78          else: 
    79              fopt = inst.bestEnergy 
    80              if inst.fglob: 
    81                  if  abs((fopt-inst.fglob)/inst.fglob) < rtol: 
    82                      answer = True  
    83              else: # safeguard if functions with fglob=0 are added by user 
    84                  if  abs(fopt) < ftol: 
    85                      answer = True  
    86          return answer 
     72         if generations: 
     73             hist = inst.energy_history 
     74             lg = len(hist) 
     75             return lg > generations and (hist[-generations]-hist[-1]) < 0 
     76         return fglob != 0 and \ 
     77               abs((inst.bestEnergy-fglob)/fglob) < tolerance 
    8778    return _ 
    8879 
    89 def AnnealTermination(feps=1e-6): 
    90     """Terminate if abs((population-x0)/x0) < feps""" 
     80def VTRChangeOverGenerations(costtolerance = 0.005, gentolerance = 1e-6, \ 
     81                             generations = 250): 
     82    """Terminate if change in cost is < gentolerance over a number of  
     83generations, or cost of last iteration is < costtolerance: 
     84 
     85cost[-g] - cost[-1] < gentolerance, where g=generations or  
     86cost[-1] < costtolerance.""" 
     87    def _(inst): 
     88         hist = inst.energy_history 
     89         lg = len(hist) 
     90         return (lg > generations and (hist[-generations]-hist[-1]) < gentolerance)\ 
     91                or ( hist[-1] < costtolerance ) 
     92    return _ 
     93 
     94def AveragePopulationImprovement(tolerance=1e-6): 
     95    """Terminate if abs((population-population[0])/population[0]) < tolerance""" 
    9196    def _(inst): 
    9297         sim = inst.population 
    93          return all(abs((sim-sim[0])/sim[0]) < feps) 
     98         return all(abs((sim-sim[0])/sim[0]) < tolerance) 
    9499    return _ 
    95100 
    96 def GradientTermination(gtol=1e-5, norm=Inf):  
    97     """Terminate if abs(gfk)**norm**1/norm <= gtol 
    98 where gfk is the gradient value""" 
     101def GradientTermination(tolerance=1e-5, norm=Inf):  
     102    """Terminate based on the current gradient value gfk.""" 
    99103    def _(inst): 
    100104        gfk = inst.gfk 
    … …  
    105109        else: 
    106110            gnorm= numpy.sum(abs(gfk)**norm,axis=0)**(1.0/norm) 
    107         return gnorm <= gtol 
     111        return gnorm <= tolerance 
    108112          
    109113    return _ 
    110114 
    111 def XTermination(tolerance=1e-5):  # Rename me? 
     115def SolutionImprovement(tolerance=1e-5):   
    112116    """Terminate if the difference between the last two  
    113 best solutions abs(bestSolution - trialSolution) is less than tolerance.""" 
     117best solutions, abs(bestSolution - trialSolution), is less than tolerance.""" 
    114118    def _(inst): 
    115119        update = inst.bestSolution - inst.trialSolution 
Note: See TracChangeset for help on using the changeset viewer.