tMOPSO

The tuning multi-objective particle swarm optimization (tMOPSO) algorithm tunes an optimization algorithms to a single problem for multiple objective function evaluation (OFE) budgets. tMOPSO finds multiple control parameter value (CPV) tuples each of which is optimal for a different OFE budget.

class optTune.tMOPSO(optAlg, CPV_lb, CPV_ub, CPV_validity_checks, sampleSizes, gammaBudget, OFE_budgets=None, OFE_max=None, extra_termination_critea=, []N=10, w=0.2, c_g=2.0, c_p=2.0, c_beta=0.1, resampling_interruption_confidence=0.9, resampling_interruption_mode='reduce_max', OFE_assessment_overshoot_function=<linearFunction y= 1.500000 * x + 100.000000 >, OFE_assessment_undershoot_function=<linearFunction y= 0.000000 * x + 0.000000 >, constrain_to_initilization_bounds=False, saveTo=None, saveInterval=10, paretoArchive_local_use_history_info=True, printFunction=<function to_stdout at 0x2ac3c2ca30c8>, printLevel=2, addtoBatch=<function _passfunction at 0x2ac3bb995230>, processBatch=<function _passfunction at 0x2ac3bb995230>, post_iteration_function=<function _passfunction at 0x2ac3bb995230>, record_V_hist=True, record_X_hist=True)
Required Args
  • optAlg - function which calls optimization algorithm or numerical method and returns two lists. The first list optAlg should return the utility measure such as the solution error (i.e f_found_opt - f_min) and should be decreasing, and the second list the number of objective function evaluations (OFEs) used in order to determine each element in the first list and should be increasing. The input arguments passed from tMOPSO to optAlg are ( numpy.array([CPV_1, CPV_2, ..., ]), OFE_budgets, randomSeed ). If OFE_budgets is a numpy.array then a solution error for each value in the list is desired, else if OFE_budgets is a integer, the solutions errors for every iteration up until OFE_budgets is desired. The type of OFE_budgets depends upon the OFE_budgets == None, if so then integer, else values from OFE_budgets is passed into optAlg.
  • CPV_lb - initialization lower bound for CPV tuning, i.e. numpy.array([CPV_1_init_lb, CPV_2_init_lb, ...])
  • CPV_ub - numpy.array([CPV_1_init_ub, CPV_2_init_ub, ...])
  • CPV_validity_checks - function used to check validity of candidate CPVs. Usage CPV_validity_checks(CPV_array, OFE_budget) returns tuple (valid, msg) where msg is string of why invalid. Should be a cheap function, as candidate CPVs are regenerated if they do not satisfy it. Use to prevent negative population sizes, populations size larger then OFE_budget checks, etcetera.
  • sampleSizes - sample sizes used to generate and refine CPV utility values. For example if the sampleSizes are [5,15,30] then all candidate CPVs will be sampled 5 times, then the possibly not dominated CPVs are then sampled another 15 times, and if still promising the for another 30 times. CPV making it to the final iteration are therefore averaged over 50 independent runs.
  • gammaBudget - the number of application layer evaluations (evaluation of the function optAlg optimizes) allocated for the tuning. NB include repeats, i.e. assessing optAlg for on OFE budget of 100 at 5 repeats, counts as a gamma of 500.
  • OFE_budgets - numpy.array of OFE budgets for which the optAlg is to be tuned under. If None then algorithm is tuned under every OFE budget upto OFE_max.
  • OFE_max - maximum OFE budget of interest. Need not if specified if OFE_budgets specified
Optional Args
  • extra_termination_critea - termination criteria in addition to gammaBudget termination criteria.
  • N - tMOPSO population size
  • w - tMOPSO particle inertia factor
  • c_g - parameter controlling the attraction towards the global guide
  • c_p - parameter controlling the attraction towards the particles personal guide
  • c_beta - particle target OFE budget perturbation factor [0,1], influences each particle velocity in the OFE budget dimension, and the local and global guide selection points.
  • resampling_interruption_confidence - re-sampling interruption confidence level used by paretoArchive2D
  • resampling_interruption_mode - choices=[‘reduce_max’, ‘check_all’]
  • OFE__assessment_overshoot_function - when assessing a CPV tuple for a OFE budget of beta, this factor is used to control overshoot, beta_actual = OFE__assessment_undershot_function(beta)
  • OFE__assessment_undershoot_function - like OFE__assessment_overshoot_function except control minimum value
  • saveTo - save optimization to this file after the optimization has been complete, at the interval specified by save_interval, use None to disable saving
  • saveInterval - optimization is saved every save_interval iterations. Zero for no saving during optimization
  • boundConstrained - should CPV be constrained between initialization bounds CPV_lb and CPV_ub
  • paretoArchive_local_use_history_info - use history info from solution error calculations to update Particles local approximations of the Pareto front. This boolean if True, may speed up tuning by increasing the quality of the each particles approximation of the Pareto Front. However, do this increase the computational overhead of tMOPSO.
  • printLevel - 0 no printing, 1 only warnings, 2 overall info, 3 lots of info, 4 verbose (not implemented)
  • addtoBatch - optAlg inputs are passed to this function before optAlg is called,
  • processBatch - function is called after all addtoBatch calls have been made, and before any optAlg calls have been made. If used then optAlg, should be retrieve solutions from this functions results
  • post_iteration_function - at the end of each iteration this function is called with tMOPSO instance as the only arg.

Examples

pure Python example

"""
Tune the simulated annealing algorithm from the scipy package to the 5D Rosenbrock problem, for multiple objective function evaluation budgets simulatenously.
The contorl parameter values (CPVs)  tuned are those which control are those which control the rate of exploration versus the rate of exploitation, namely 
 * m - control how quickly the Temperature cools down, i.e. higher m = higher exploration and lower exploitation.
        c = m * exp(-n * quench)
        T_new = T0 * exp(-c * k**quench)
 * dwell - evaluation at each tempreture
"""
import numpy
from scipy import optimize
from optTune import tMOPSO, evaluation_history_recording_wrapper, get_F_vals_at_specified_OFE_budgets, linearFunction


def Ros_ND(x) :
    "gerneralised Rossenbroch function"
    return sum([100*(x[ii+1]-x[ii]**2)**2 + (1-x[ii])**2 for ii in range(len(x)-1)])

def solution_valid(x):
    return ( -2.048 <= x ).all() and ( x <= 2.048 ).all()

def run_simulated_annealing(CPVs, OFE_budgets, randomSeed):
    dwell = int(CPVs[0]) #equibavent to population size in evolutionary algorithms
    func = evaluation_history_recording_wrapper( Ros_ND, dwell, solution_valid )
    optimize.anneal(func, 
                    x0 = -0.5 * numpy.random.rand(5),
                    #x0 = -2.048 + 2*2.048*numpy.random.rand(10), #if used make sure tMOPSO sample size greater than 100
                    m = CPVs[1],
                    T0 = 500.0,
                    lower= -2.048,
                    upper=  2.048,
                    dwell=dwell, 
                    maxeval = max( OFE_budgets ), #termination criteria
                    feps = 0.0, 
                    Tf = 0.0)
    return get_F_vals_at_specified_OFE_budgets(F=func.f_hist, E=func.OFE_hist, E_desired=OFE_budgets)

def CPV_valid(CPVs, OFE_budget):
    if CPVs[0] < 5:
        return False,'dwell,CPVs[0] < 5'
    if CPVs[1] < 0.0001:
        return False,'CPVs[1] < 0.0001'
    return True,''

tuningOpt = tMOPSO( 
    optAlg = run_simulated_annealing, 
    CPV_lb = numpy.array([10, 0.0]), 
    CPV_ub = numpy.array([50, 5.0]),
    CPV_validity_checks = CPV_valid,
    OFE_budgets=numpy.logspace(1,3,30).astype(int),
    sampleSizes = [2,8,20], #resampling size of 30
    resampling_interruption_confidence = 0.6,
    gammaBudget = 30*1000*50, #increase to get a smoother result ...
    OFE_assessment_overshoot_function = linearFunction(2, 100 ),
    N = 10,
    )
print(tuningOpt)


Fmin_values = [ d.fv[1] for d in tuningOpt.PFA.designs ] 
OFE_budgets = [ d.fv[0] for d in tuningOpt.PFA.designs ]
dwell_values =  [ int(d.xv[1]) for d in tuningOpt.PFA.designs ]
m_values = [ d.xv[2] for d in tuningOpt.PFA.designs ]

print('OFE budget     Fmin      dwell       m     ')
for a,b,c,d in zip(OFE_budgets, Fmin_values, dwell_values, m_values):
    print('  %i       %6.4f      %i       %4.2f' % (a,b,c,d))

from matplotlib import pyplot
p1 = pyplot.semilogx(OFE_budgets, dwell_values, 'g.')[0]
pyplot.ylabel('dwell')
pyplot.ylim( min(dwell_values) - 1, max( dwell_values) + 1)
pyplot.twinx()
p2 = pyplot.semilogx(OFE_budgets, m_values, 'bx')[0]
pyplot.ylim( 0, max(m_values)*1.1)
pyplot.ylabel('m (rate of cool)')
pyplot.legend([p1,p2],['dwell','m'], loc='upper center')
pyplot.xlabel('OFE budget')
pyplot.xlim(min(OFE_budgets)-1,max(OFE_budgets)+60)
pyplot.title('Optimal CPVs for different OFE budgets')
#pyplot.figure()
#tuningOpt.plot_HV_hist()
#pyplot.title('Hyper Volume History of multi-objective tuning optimization')

pyplot.show()

f2py example

A faster way of the doing the above examples is to use f2py, and then call the resulting shared object from Python.

Fortran code

! written by Antoine Dymond, Jan 2012, based upon the Scipy.optimize.anneal code, the scipy license agreement is as follows:
!!$
!!$Copyright (c) 2001, 2002 Enthought, Inc.
!!$All rights reserved.
!!$
!!$Copyright (c) 2003-2009 SciPy Developers.
!!$All rights reserved.
!!$
!!$Redistribution and use in source and binary forms, with or without
!!$modification, are permitted provided that the following conditions are met:
!!$
!!$  a. Redistributions of source code must retain the above copyright notice,
!!$     this list of conditions and the following disclaimer.
!!$  b. Redistributions in binary form must reproduce the above copyright
!!$     notice, this list of conditions and the following disclaimer in the
!!$     documentation and/or other materials provided with the distribution.
!!$  c. Neither the name of the Enthought nor the names of its contributors
!!$     may be used to endorse or promote products derived from this software
!!$     without specific prior written permission.
!!$
!!$
!!$THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
!!$AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
!!$IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
!!$ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
!!$ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
!!$DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
!!$SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
!!$CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
!!$LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
!!$OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
!!$DAMAGE.
!
! compile using
! $ f2py -c -m anneal_fortran anneal.f90


module anneal_module
  ! to compile using f2py, leave real type as default, needs to be compadiable between numpy and fortran.
  implicit None  ! fortran is case insensitive.
  !f2py forces all variable to lower case in linking module
  integer :: problem_id, objFun_evals
  real(8), allocatable, dimension(:) :: fval_hist, x_min
  integer, allocatable, dimension(:) :: eval_hist

contains

  function objFun(x, D)
    integer :: D, i
    real(8) :: x(D), objFun
    objFun = 0
    objFun_evals = objFun_evals + 1
    select case(problem_id)
    case (1) !general Rossenbrock
       do i = 1,D-1
          objFun = objFun + 100.0*(x(i+1) - x(i)**2.0)**2.0 + (1.0 - x(i))**2.0
       end do
    case (2) !sphere
       objFun = sum(x**2)
    end select
  end function objFun

  SUBROUTINE set_random_seed(seed_offset)
    !taken from the gfortran site. 
    !http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html#RANDOM_005fSEED
    INTEGER :: i, n, seed_offset
    INTEGER, DIMENSION(:), ALLOCATABLE :: seed   
    CALL RANDOM_SEED(size = n)
    ALLOCATE(seed(n))
    seed = seed_offset + 37 * (/ (i - 1, i = 1, n) /)
    CALL RANDOM_SEED(PUT = seed)
    DEALLOCATE(seed)
  END SUBROUTINE set_random_seed


  function rand_uniform()
    real(8) :: rand_uniform
    call random_number(rand_uniform)
  end function rand_uniform

  subroutine fast_sa_run(prob_id, x0, D, T0, dwell, m, n, quench, boltzmann, maxEvals, lower, upper, random_seed)
    integer :: prob_id, D, dwell, maxEvals, random_seed
    real(8) :: T0, m, n, quench, boltzmann
    real(8), dimension(D) :: x0, u, lower, upper 
    real(8) :: T, c, dE, y_j
    integer :: k, kMax, acceptTest, i, j
    real(8) :: current_state_fv, last_state_fv, best_state_fv
    real(8), dimension(D) ::current_state_xv, last_state_xv, best_state_xv


    call set_random_seed( random_seed )
    problem_id  = prob_id
    kMax = maxEvals/ dwell
    if (allocated(fval_hist)) deallocate(fval_hist, eval_hist )
    allocate(fval_hist(kMax), eval_hist(kMax))

    fval_hist = 0.0
    objFun_evals = 0
    T = T0
    c = m * exp(n * quench)
    last_state_xv = x0
    last_state_fv = objFun(x0, D)
    best_state_xv = x0
    best_state_fv = last_state_fv
    do k = 1,kMax
       do i = 1, dwell
          ! schedule.update_guess
          call random_number(u)
          do j = 1,D
             y_j = T * ( ( 1+1.0/T) ** abs(2*u(i)-1) - 1.0 )
             if ( (u(j) - 0.5) < 0 ) y_j = - y_j
             current_state_xv(j) = last_state_xv(j) + y_j*(upper(j) - lower(j))
             if (current_state_xv(j) < lower(j)) current_state_xv(j) = lower(j)
             if (current_state_xv(j) > upper(j)) current_state_xv(j) = upper(j)
          end do
          current_state_fv =  objFun(current_state_xv, D)
          dE = current_state_fv - last_state_fv
          !schedule.accept_test(dE)
          acceptTest = 0
          if (dE < 0) then
             acceptTest = 1
          else
             if  ( exp(-dE* 1.0 / boltzmann /T ) > rand_uniform() ) acceptTest = 1
          end if
          if ( acceptTest == 1) then
             last_state_xv = current_state_xv
             last_state_fv = current_state_fv
             if (last_state_fv < best_state_fv) then
                best_state_xv = last_state_xv
                best_state_fv = last_state_fv
             end if
          end if
       end do
       !tempreture update
       T = T0 * exp(-c*k*quench)
       !book keeping
       fval_hist(k) = best_state_fv
       eval_hist(k) = objFun_evals
    end do
    if (allocated(x_min)) deallocate(x_min)
    allocate(x_min(D))
    x_min = best_state_xv
  end subroutine fast_sa_run

end module anneal_module

Python code

"""
Tune the simulated annealing algorithm from the scipy package to the generalized Rosenbrock problem, for multiple objective function evaluation (OFE) budgets simulatenously.
Same as the other example, except a fortran version of fast sa is used.
"""
import numpy, os
from optTune import tMOPSO, get_F_vals_at_specified_OFE_budgets, linearFunction
print('Please note this example only works on Linux, and requires gfortran')
if not os.path.exists('anneal_fortran.so'):
    os.system('f2py -c -m anneal_fortran anneal.f90')
from anneal_fortran import anneal_module
D = 5 #number of dimensions for Rosenbrock problem

def anneal(CPVs, OFE_budgets, randomSeed):
    #fast_sa_run - Function signature:
    #  fast_sa_run(prob_id,x0,t0,dwell,m,n,quench,boltzmann,maxevals,lower,upper,random_seed,[d])
    anneal_module.fast_sa_run(prob_id = 1 ,
                              x0 = -2.048 + 2*2.048*numpy.random.rand(D),
                              t0 = 500.0,
                              dwell = int(CPVs[0]),
                              m = CPVs[1],
                              n = 1.0,
                              quench = 1.0,
                              boltzmann = 1.0,
                              maxevals = max(OFE_budgets),
                              lower = -2.048*numpy.ones(D),
                              upper =  2.048*numpy.ones(D),
                              random_seed = randomSeed)
    return get_F_vals_at_specified_OFE_budgets(F=anneal_module.fval_hist.copy(), E=anneal_module.eval_hist.copy(), E_desired=OFE_budgets)

def CPV_valid(CPVs, OFE_budget):
    if CPVs[0] < 5:
        return False,'dwell,CPVs[0] < 5'
    if CPVs[1] < 0.0001:
        return False,'CPVs[1] < 0.0001'
    return True,''

tuningOpt = tMOPSO( 
    optAlg = anneal, 
    CPV_lb = numpy.array([10, 0.0]), 
    CPV_ub = numpy.array([50, 5.0]),
    CPV_validity_checks = CPV_valid,
    OFE_budgets=numpy.logspace(1,3,30).astype(int),
    sampleSizes = [2,8,20], #resampling size of 30
    resampling_interruption_confidence = 0.6,
    gammaBudget = 30*1000*50, #increase to get a smoother result ...
    OFE_assessment_overshoot_function = linearFunction(2, 100 ),
    N = 10,
    printLevel=1,
    )
print(tuningOpt)

Fmin_values = [ d.fv[1] for d in tuningOpt.PFA.designs ] 
OFE_budgets = [ d.fv[0] for d in tuningOpt.PFA.designs ]
dwell_values =  [ int(d.xv[1]) for d in tuningOpt.PFA.designs ]
m_values = [ d.xv[2] for d in tuningOpt.PFA.designs ]

print('OFE budget     Fmin      dwell       m     ')
for a,b,c,d in zip(OFE_budgets, Fmin_values, dwell_values, m_values):
    print('  %i       %6.4f      %i       %4.2f' % (a,b,c,d))

from matplotlib import pyplot
p1 = pyplot.semilogx(OFE_budgets, dwell_values, 'g.')[0]
pyplot.ylabel('dwell')
pyplot.ylim( min(dwell_values) - 1, max( dwell_values) + 1)
pyplot.twinx()
p2 = pyplot.semilogx(OFE_budgets, m_values, 'bx')[0]
pyplot.ylim( 0, max(m_values)*1.1)
pyplot.ylabel('m (rate of cool)')
pyplot.legend([p1,p2],['dwell','m'], loc='best')
pyplot.xlabel('OFE budget')
pyplot.xlim(min(OFE_budgets)-1,max(OFE_budgets)+60)
pyplot.title('Optimal CPVs for different OFE budgets')

pyplot.show()

Table Of Contents

Previous topic

Examples

Next topic

tPSO

This Page