Changeset 104


Ignore:
Timestamp:
02/20/09 18:36:22 (7 years ago)
Author:
mmckerns
Message:

moved fmin_powell to scipy_optimize;
removed scipy dependency;
cleanup of Solver documentation

Files:
1 added
1 deleted
6 edited

Legend:

Unmodified
Added
Removed
  • DEV_NOTES

    r103 r104  
    55*/*/*py: 
    66 - move "circle" to mystic.models, and move more of "mogi" to mystic.models 
    7 mystic/mystic/fmin_powell.py: 
     7mystic/mystic/scipy_optimize.py: 
    88 - derive an AbstractSolver class for all solvers to inherit from 
    9  - add scipy solver "Brent", to _remove_ the scipy dependency (*) 
    10  - migrate all contents into scipy_optimize.py 
    11  - clean up & edit inline documentation 
    12  - better way to deal with callback, direc, maxfun, maxiter (?) 
    13 mystic/mystic/scipy_optimize.py: 
    14  - add the "__all__" scipy trick 
     9 - better way to deal with direc, maxfun, maxiter (?) 
    1510mystic/mystic/differential_evolution.py: 
    1611 - add ability to choose initial point 
  • mystic/examples/test_rosenbrock3.py

    r103 r104  
    2929 
    3030   #from scipy.optimize import fmin_powell 
    31     from mystic.fmin_powell import fmin_powell, PowellDirectionalSolver 
     31    from mystic.scipy_optimize import fmin_powell, PowellDirectionalSolver 
    3232   #print fmin_powell(rosen,x0,retall=0,full_output=0) 
    3333    solver = PowellDirectionalSolver(len(x0)) 
  • mystic/mystic/Make.mm

    r103 r104  
    2929    differential_evolution.py \ 
    3030    scipy_optimize.py \ 
    31     fmin_powell.py \ 
     31    _scipy060optimize.py \ 
    3232    termination.py \ 
    3333    strategy.py \ 
  • mystic/mystic/__init__.py

    r103 r104  
    1616# solvers 
    1717import differential_evolution, scipy_optimize 
    18 import fmin_powell #FIXME: move to scipy_optimize after remove scipy dependency 
    1918 
    2019# strategies, termination conditions 
  • mystic/mystic/differential_evolution.py

    r102 r104  
    4242 
    4343""" 
     44__all__ = ['DifferentialEvolutionSolver','DifferentialEvolutionSolver2'] 
    4445 
    4546from mystic.tools import Null, wrap_function 
  • mystic/mystic/scipy_optimize.py

    r103 r104  
    44# (derives from optimize.py module by Travis E. Oliphant) 
    55# 
    6 # adapted from scipy.optimize.fmin, scipy version 0.4.8 
     6# adapted scipy.optimize.fmin (from scipy version 0.4.8) 
    77# by Patrick Hung, Caltech. 
    88# 
    9 # adapted to class (& added bounds) 
    10 # updated to scipy version 0.6.0 
     9# adapted from function to class (& added bounds) 
     10# adapted scipy.optimize.fmin_powell 
     11# updated solvers to scipy version 0.6.0 
    1112# by mmckerns@caltech.edu 
    1213 
    1314""" 
    14 Algorithm adapted from scipy.optimize,  
     15Algorithms adapted from scipy.optimize,  
    1516 
    1617fmin        ---      Nelder-Mead Simplex algorithm. 
     
    1920                        -- StepMonitor 
    2021                      
     22fmin_powell ---      Powell's Directional Search method. 
     23                     Takes two additional input args 
     24                        -- EvaluationMonitor 
     25                        -- StepMonitor 
     26                      
    2127""" 
     28__all__ = ['NelderMeadSimplexSolver','PowellDirectionalSolver', 
     29           'fmin','fmin_powell'] 
     30 
    2231 
    2332from mystic.tools import Null, wrap_function 
     
    2736from numpy import atleast_1d, eye, zeros, shape, \ 
    2837     asarray, absolute, sqrt, Inf, asfarray 
    29 from numpy import clip 
     38from numpy import clip, squeeze 
    3039 
    3140abs = absolute 
     41 
     42from _scipy060optimize import brent #XXX: local copy to avoid dependency! 
     43 
    3244 
    3345class NelderMeadSimplexSolver(object): 
     
    176188      ExtraArgs -- extra arguments for func. 
    177189 
    178           """ 
     190    Further Inputs: 
     191 
     192      callback -- an optional user-supplied function to call after each 
     193                  iteration.  It is called as callback(xk), where xk is the 
     194                  current parameter vector 
     195      disp -- non-zero to print convergence messages. 
     196      radius -- percentage change for initial simplex values 
     197 
     198""" 
    179199        # set arg names to scipy.optimize.fmin names; set fixed inputs 
    180200        x0 = self.population[0] 
     
    185205        callback=None  #user-supplied function, called after each step 
    186206        radius=0.05    #percentage change for initial simplex values 
     207        if kwds.has_key('callback'): callback = kwds['callback'] 
    187208        if kwds.has_key('disp'): disp = kwds['disp'] 
    188209        if kwds.has_key('radius'): radius = kwds['radius'] 
    189         if kwds.has_key('callback'): callback = kwds['callback'] 
    190210        #------------------------------------------------------------- 
    191211 
     
    431451    return retlist 
    432452 
     453############################################################################ 
     454 
     455def _linesearch_powell(func, p, xi, tol=1e-3): 
     456    # line-search algorithm using fminbound 
     457    #  find the minimium of the function 
     458    #  func(x0+ alpha*direc) 
     459    def myfunc(alpha): 
     460        return func(p + alpha * xi) 
     461    alpha_min, fret, iter, num = brent(myfunc, full_output=1, tol=tol) 
     462    xi = alpha_min*xi 
     463    return squeeze(fret), p+xi, xi 
     464 
     465 
     466class PowellDirectionalSolver(NelderMeadSimplexSolver): #FIXME: not a simplex solver 
     467    """ 
     468    Powell Direction Search optimization adapted from scipy.optimize.fmin_powell. 
     469    """ 
     470     
     471    def __init__(self, dim): 
     472        """ 
     473 Takes one initial input:  
     474   dim      -- dimensionality of the problem 
     475        """ 
     476        NelderMeadSimplexSolver.__init__(self,dim) 
     477        self._direc = None #FIXME: this is the easy way to return 'direc'... 
     478       #FIXME: NO SIMPLEX, so maybe best not to inherit __init__? 
     479 
     480 
     481    def Solve(self, func, termination, 
     482              maxiter=None, maxfun=None, sigint_callback=None, 
     483              EvaluationMonitor=Null, StepMonitor=Null, ExtraArgs=(), **kwds): 
     484        """Minimize a function using modified Powell's method. 
     485 
     486    Description: 
     487 
     488      Uses a modified Powell Directional Search algorithm to find 
     489      the minimum of function of one or more variables. 
     490 
     491    Inputs: 
     492 
     493      func -- the Python function or method to be minimized. 
     494      termination -- callable object providing termination conditions. 
     495 
     496    Additional Inputs: 
     497 
     498      maxiter -- the maximum number of iterations to perform. 
     499      maxfun -- the maximum number of function evaluations. 
     500      sigint_callback -- callback function for signal handler. 
     501      EvaluationMonitor -- a callable object that will be passed x, fval 
     502           whenever the cost function is evaluated. 
     503      StepMonitor -- a callable object that will be passed x, fval 
     504           after the end of a simplex iteration. 
     505      ExtraArgs -- extra arguments for func. 
     506 
     507    Further Inputs: 
     508 
     509      callback -- an optional user-supplied function to call after each 
     510                  iteration.  It is called as callback(xk), where xk is the 
     511                  current parameter vector 
     512      direc -- initial direction set 
     513      xtol -- line-search error tolerance. 
     514      disp -- non-zero to print convergence messages. 
     515 
     516""" 
     517 
     518        # set arg names to scipy.optimize.fmin_powell names; set fixed inputs 
     519        x0 = self.population[0] 
     520        args = ExtraArgs 
     521        full_output=1  #non-zero if fval and warnflag outputs are desired. 
     522        disp=0         #non-zero to print convergence messages. 
     523        retall=0       #non-zero to return all steps 
     524        direc=None 
     525        callback=None  #user-supplied function, called after each step 
     526        xtol=1e-4      #line-search error tolerance 
     527        if kwds.has_key('callback'): callback = kwds['callback'] 
     528        if kwds.has_key('direc'): direc = kwds['direc']  #XXX: best interface? 
     529        if kwds.has_key('xtol'): xtol = kwds['xtol'] 
     530        if kwds.has_key('disp'): disp = kwds['disp'] 
     531        #------------------------------------------------------------- 
     532 
     533        import signal 
     534        import mystic.termination as detools 
     535        detools.EARLYEXIT = False 
     536 
     537        fcalls, func = wrap_function(func, args, EvaluationMonitor) 
     538        if self._useStrictRange: 
     539            x0 = self._setGuessWithinRangeBoundary(x0) 
     540            func = wrap_bounds(func, self._strictMin, self._strictMax) 
     541 
     542        def handler(signum, frame): 
     543            import inspect 
     544            print inspect.getframeinfo(frame) 
     545            print inspect.trace() 
     546            while 1: 
     547                s = raw_input(\ 
     548""" 
     549  
     550 Enter sense switch. 
     551 
     552   sol: Write current best solution. 
     553   cont: Continue calculation. 
     554   call: Executes sigint_callback [%s]. 
     555   exit: Exits with current best solution. 
     556 
     557 >>> """ % sigint_callback) 
     558                if s.lower() == 'sol':  
     559                    print "sw1." 
     560                    print self.bestSolution 
     561                elif s.lower() == 'cont':  
     562                    return 
     563                elif s.lower() == 'call':  
     564                    # sigint call_back 
     565                    if sigint_callback is not None: 
     566                        sigint_callback(self.bestSolution) 
     567                elif s.lower() == 'exit':  
     568                    detools.EARLYEXIT = True 
     569                    return 
     570                else: 
     571                    print "unknown option : %s ", s 
     572 
     573        self.signal_handler = handler 
     574 
     575        if self._handle_sigint: signal.signal(signal.SIGINT, self.signal_handler) 
     576        #------------------------------------------------------------- 
     577 
     578        x = asfarray(x0).flatten() 
     579        if retall: 
     580            allvecs = [x] 
     581        N = len(x) #XXX: this should be equal to self.nDim 
     582        rank = len(x.shape) 
     583        if not -1 < rank < 2: 
     584            raise ValueError, "Initial guess must be a scalar or rank-1 sequence." 
     585        if maxiter is None: 
     586            maxiter = N * 1000 
     587        if maxfun is None: 
     588            maxfun = N * 1000 
     589        self._maxiter = maxiter #XXX: better to just copy the code? 
     590        self._maxfun = maxfun   #XXX: better to just copy the code? 
     591 
     592        if direc is None: 
     593            direc = eye(N, dtype=float) 
     594        else: 
     595            direc = asarray(direc, dtype=float) 
     596        fval = squeeze(func(x)) 
     597        x1 = x.copy() 
     598 
     599        self._direc = direc #XXX: instead, use a monitor? 
     600        self.bestSolution = x 
     601        self.bestEnergy = fval 
     602        self.population = [x]    #XXX: pointless, if simplex not used 
     603        self.popEnergy = [fval]  #XXX: pointless, if simplex not used 
     604        self.energy_history.append(self.bestEnergy) 
     605        StepMonitor(x,fval) # get initial values 
     606 
     607        iter = 0; 
     608        ilist = range(N) 
     609 
     610        CONTINUE = True 
     611        while CONTINUE: 
     612            fx = fval 
     613            bigind = 0 
     614            delta = 0.0 
     615            for i in ilist: 
     616                direc1 = direc[i] 
     617                fx2 = fval 
     618                fval, x, direc1 = _linesearch_powell(func, x, direc1, tol=xtol*100) 
     619                if (fx2 - fval) > delta: 
     620                    delta = fx2 - fval 
     621                    bigind = i 
     622 
     623            iter += 1 
     624            if callback is not None: 
     625                callback(x) 
     626            if retall: 
     627                allvecs.append(x) 
     628 
     629            self.energy_history.append(fval) #XXX: the 'best' for now... 
     630            if detools.EARLYEXIT or termination(self): CONTINUE = False #break 
     631            elif fcalls[0] >= maxfun: CONTINUE = False #break 
     632            elif iter >= maxiter: CONTINUE = False #break 
     633 
     634            else: # Construct the extrapolated point 
     635                direc1 = x - x1 
     636                x2 = 2*x - x1 
     637                x1 = x.copy() 
     638                fx2 = squeeze(func(x2)) 
     639     
     640                if (fx > fx2): 
     641                    t = 2.0*(fx+fx2-2.0*fval) 
     642                    temp = (fx-fval-delta) 
     643                    t *= temp*temp 
     644                    temp = fx-fx2 
     645                    t -= delta*temp*temp 
     646                    if t < 0.0: 
     647                        fval, x, direc1 = _linesearch_powell(func, x, direc1, tol=xtol*100) 
     648                        direc[bigind] = direc[-1] 
     649                        direc[-1] = direc1 
     650 
     651                self.energy_history[-1] = fval #...update to 'best' energy 
     652 
     653            self._direc = direc #XXX: instead, use a monitor? 
     654            self.bestSolution = x 
     655            self.bestEnergy = fval 
     656            self.population = [x]    #XXX: pointless, if simplex not used 
     657            self.popEnergy = [fval]  #XXX: pointless, if simplex not used 
     658            StepMonitor(x,fval) # get ith values; #XXX: should be [x],[fval] ? 
     659     
     660        self.generations = iter 
     661        signal.signal(signal.SIGINT,signal.default_int_handler) 
     662 
     663        # code below here is dead, unless disp!=0 
     664        warnflag = 0 
     665 
     666        if fcalls[0] >= maxfun: 
     667            warnflag = 1 
     668            if disp: 
     669                print "Warning: Maximum number of function evaluations has "\ 
     670                      "been exceeded." 
     671        elif iter >= maxiter: 
     672            warnflag = 2 
     673            if disp: 
     674                print "Warning: Maximum number of iterations has been exceeded" 
     675        else: 
     676            if disp: 
     677                print "Optimization terminated successfully." 
     678                print "         Current function value: %f" % fval 
     679                print "         Iterations: %d" % iter 
     680                print "         Function evaluations: %d" % fcalls[0] 
     681     
     682        x = squeeze(x) 
     683 
     684        if full_output: 
     685            retlist = x, fval, direc, iter, fcalls[0], warnflag 
     686            if retall: 
     687                retlist += (allvecs,) 
     688        else: 
     689            retlist = x 
     690            if retall: 
     691                retlist = (x, allvecs) 
     692 
     693        return #retlist 
     694 
     695 
     696def fmin_powell(func, x0, args=(), xtol=1e-4, ftol=1e-4, maxiter=None, 
     697                maxfun=None, full_output=0, disp=1, retall=0, callback=None, 
     698                direc=None): 
     699    """fmin_powell using the 'original' scipy.optimize.fmin_powell interface""" 
     700 
     701    #FIXME: need to resolve "direc" 
     702    #        - should just pass 'direc', and then hands-off ?  How return it ? 
     703 
     704    from mystic.tools import Sow 
     705    stepmon = Sow() 
     706    evalmon = Sow() 
     707    from mystic.termination import NormalizedChangeOverGeneration as NCOG 
     708 
     709    solver = PowellDirectionalSolver(len(x0)) 
     710    solver.SetInitialPoints(x0) 
     711   #solver.enable_signal_handler() 
     712    solver.Solve(func,termination=NCOG(ftol),\ 
     713                 maxiter=maxiter,maxfun=maxfun,\ 
     714                 EvaluationMonitor=evalmon,StepMonitor=stepmon,\ 
     715                 xtol=xtol, callback=callback, \ 
     716                 disp=disp, direc=direc)   #XXX: last two lines use **kwds 
     717    solution = solver.Solution() 
     718 
     719    # code below here pushes output to scipy.optimize.fmin_powell interface 
     720   #x = list(solver.bestSolution) 
     721    x = solver.bestSolution 
     722    fval = solver.bestEnergy 
     723    warnflag = 0 
     724    fcalls = len(evalmon.x) 
     725    iterations = len(stepmon.x) - 1 
     726    allvecs = [] 
     727    for i in range(len(stepmon.x)): 
     728       #allvecs.append(list(stepmon.x[i])) 
     729        allvecs.append(stepmon.x[i]) 
     730    direc = solver._direc #FIXME: better way to get direc from Solve() ? 
     731 
     732    if fcalls >= solver._maxfun: 
     733        warnflag = 1 
     734    elif iterations >= solver._maxiter: 
     735        warnflag = 2 
     736 
     737    x = squeeze(x) #FIXME: write squeezed x to stepmon instead? 
     738 
     739    if full_output: 
     740        retlist = x, fval, direc, iterations, fcalls, warnflag 
     741        if retall: 
     742            retlist += (allvecs,) 
     743    else: 
     744        retlist = x 
     745        if retall: 
     746            retlist = (x, allvecs) 
     747 
     748    return retlist 
     749 
    433750 
    434751if __name__=='__main__': 
Note: See TracChangeset for help on using the changeset viewer.