source: branches/decorate/restarts.py @ 855

#!/usr/bin/env python
#
# Author: Mike McKerns (mmckerns @caltech and @uqfoundation)
# Copyright (c) 2012-2016 California Institute of Technology.
# License: 3-clause BSD.  The full license text is available at:
#  - http://mmckerns.github.io/project/mystic/browser/mystic/LICENSE

def sumt(cost, solver, term, eps = 1e-4, max_SUMT_iters = 10, **kwds):

    from mystic.tools import Null
    from numpy import inf
    # standard solver config
    StepMonitor=kwds.get('StepMonitor', Null)
    EvaluationMonitor=kwds.get('EvaluationMonitor', Null)
    sigint_callback=kwds.get('sigint_callback', None)
    disp = kwds.get('disp',False)

    # other config
    randomstate = kwds.get('randomstate',None)

    # SUMT penalty restart
    F_values = []
    E_values = []
    x_values = []
    total_iterations = 0.
    total_energy_history = []
    exitflag = 0
    # The outer SUMT for loop. In each iteration, an optimization will be 
    # performed.
    for q in range(max_SUMT_iters):
        if q == 0:
            x0 = solver.population[0]
        else:
            x0 = x_values[-1]

        if disp:
            print 'SUMT iteration %d: x = ' % q, x0

        if randomstate:
            import random
            random.setstate(randomstate)

        # 'Clear' some of the attributes of the solver
        solver.population = [[0.0 for i in range(solver.nDim)] for j in range(solver.nPop)]
       #solver.generations = 0
        solver.bestEnergy = None
        solver.bestSolution = None
        solver.trialSolution = [0.0] * solver.nDim
        solver.popEnergy = [solver._init_popEnergy] * solver.nPop
       #solver.energy_history = []

        # individual solver __init__
        if solver.__class__.__name__ == 'NelderMeadSimplexSolver':
            simplex = solver.nDim + 1
            solver.popEnergy = [solver._init_popEnergy] * simplex
            solver.population = [[0.0 for i in range(solver.nDim)] for j in range(simplex)]
        elif solver.__class__.__name__ == 'PowellDirectionalSolver':
            solver._direc = None
        elif solver.__class__.__name__ in ['DifferentialEvolutionSolver','DifferentialEvolutionSolver2']:
            solver.genealogy = [ [] for j in range(solver.nPop)]
        else: pass

        # Perform the 'inner' optimization
        solver.SetInitialPoints(x0)
        solver.SetEvaluationMonitor(EvaluationMonitor)
        solver.SetGenerationMonitor(StepMonitor)
        solver.Solve(cost, term, sigint_callback=sigint_callback, **kwds)
        x = solver.bestSolution
        y = solver.bestEnergy

        # When a user sets maxiter and maxfun, that is actually the maxiter 
        # and maxfun *each* SUMT iteration. However, it is difficult to check 
        # total maxfun and maxiter each SUMT iteration because then what 
        # would the inner maxiter and maxfun be? Better to just leave maxiter 
        # and maxfun as limits for each inner optimization...
        total_iterations += solver.generations
        total_energy_history += solver.energy_history

        F_values.append(cost(x)) #XXX: cost = no penalty; y = penalty
        E_values.append(y)
        x_values.append(x)

        # Update penalty parameters
        solver._penalty.store(x) #XXX: may want to combine store with iter

        # break due to the signal handler
        if solver._EARLYEXIT:
            exitflag = 3
            break

        # breaking if constraints satisfied may cause premature exit. so don't
        # if False:
        #     exitflag = 1
        #     break

        # break if cost(x) not improving after SUMT loop
        if q != 0:
            if abs(F_values[q] - F_values[q-1]) <= eps: #XXX: F or E ?
                exitflag = 2
                break
            # break for no improvement in x's after SUMT loop
            #if linalg.norm(asarray(x_values[q])-asarray(x_values[q-1])) <= eps:
            #    exitflag = 4
            #    break

        # Update penalty parameters
        solver._penalty.iter()   #XXX: may need 'assurance' to have iter()
        
    if disp:
        if exitflag == 0:
            print "SUMT exited because maximum number of SUMT iterations reached."
        elif exitflag == 1:
            print "SUMT exited because constraints were satisfied."
        elif exitflag == 2:
            print "SUMT exited because F values reached tolerance for improving."
        elif exitflag == 3:
            print "SUMT exited because signal handler 'exit' was selected."
        #elif exitflag == 4:
        #    print "SUMT exited because x values reached tolerance for improving."    
        else:
            print "SUMT exited for other reason."
        print "%d SUMT iterations were performed." % (q + 1)

    # Return the last solver, which contains all of the important information
    # First modify the solver to include total data.
   #solver.generations = total_iterations
   #solver.energy_history = total_energy_history

    #XXX solver.bestEnergy is not exactly cost(solver.bestSolution),
    # since it is cost(bestSolution) with a penalty. Compensate for that
    # here? Then, it is inconsistent with solver.energy_history, which 
    # also records cost(bestSolution) with a penalty....
    #solver.bestEnergy = cost(solver.bestSolution) 
    return solver