/*********************************************************
Implements the simplex-simulated annealing scheme for
continuous, bounded optimization described in:
Cardoso, M.F., Salcedo, R.L., Feyo de Azevedo, S.
The simplex-simulated annealing approach to continuous
non-linear optimization. Computers chem. Engng
20(9): 1065-1080, (1996).
Unconstrained Simulated Annealing code was written by
Andrew Davison (andrew.davison@unic.cnrs-gif.fr).
This code was modified to implement the Cardoso algorithm
by Christina Weaver (Dec. 2004), christina.weaver@mssm.edu.
This implementation is described in
Weaver, C.M. and Wearne, S.L. The role of action potential shape
and parameter constraints in optimization of compartment models.
Neurocomputing 69: 1053-1057, (2006).
*********************************************************/
mulfit_optimizers_append("Constrained Simulated Annealing", "MulfitSA_Rctr", 0)
begintemplate MulfitSA_Rctr
public prun, sa_efun, minerr, nefun, time
public saveflag, showopt, verbose
public estimate_inittemp
public alpha, beta, gamma, lambda, ftol, start_temptr, gtol
public restart_factor, cool_rate, iter_per_step, quit_temptr
public savepath_name, bd_pnlty, rest_inname, use_rst, get_rstname
// testing wraparound BCs
public wrapfile
objref wrapfile
objref r, simplex, old_simplex
objref simplex_errors, new_point, best_point, mean_point, wrap_point
objref savepath, tl, pf, verbfile, annealpath, rest_file
objref simperrors_cpy, simplex_noise, centrd, tmp_vec
strdef savefit_name, cmnd, savepath_name, savetmp_name, ann_str, verbname
strdef annealfit_name
objref code_low, code_high
// There are two restart filenames: one to read for initialization (when directed),
// and one to write. Filenames may be different, but not necessarily! Be careful
// of this.
strdef rest_name, rest_inname
//
// N.B. cmw nov 04:
//
// "eb" corresponds to the error of the "best point"
// "elo" is the lowest error of the simplex vertices
// "minerr" is the value returned to the MulRunFitter wrapper
//
// in the noisy_simplex proc, comments were added as in
// Numerical Recipes in C.
proc init() {
printf("Simulated Annealing + Recenter, last updated 12 Jun 06\n")
pf = $o1
r = new Random()
r.uniform(1e-20,1)
savepath = new File()
annealpath = new File()
verbfile = new File()
rest_file = new File()
new_point = new Vector()
nefun = 0
st = 0
lambda = 2
ftol = 1e-5
quit_temptr = 1e-7
gtol = ftol
alpha = 1.0
beta = 2.0
gamma = 0.5
iter_per_step = 100
start_temptr = 1000
cool_rate = 0.9
restart_factor = -1
eb = 1e9
enow = 1e9
saveflag = 1
verbose = 0
sprint(savepath_name,"annealBD")
sprint(rest_inname,"annealBD.rst")
bd_pnlty = 1e5
set_annsched(0)
// for 'complex' cooling schedule from A. Davison, use set_annsched(1)
use_NE = 0
use_rst = 0
Kbound = Kbound_m2 = 1
anneal_count = 1
}
/**************************************************
Evaluate the error between model & target
**************************************************/
func sa_efun() { local e, ttime, i
// $3 = 1 if this is centroid/out-of-bounds, 0 otherwise
nefun += 1
if( verbose ) {
verbfile.printf("SA_EFUN : centroid/OOB = %d\n",$3)
new_point.printf(verbfile)
verbfile.printf("vs. sa_efun sent ")
for i = 0,$1-1 {
verbfile.printf("%-12.8g ",$&2[i])
}
verbfile.printf("\n")
verbfile.printf("vs. pf.pmlst point ")
for i = 0,pf.parmlist.count-1 {
verbfile.printf("%-12.8g ",pf.parmlist.object(i).val)
}
verbfile.printf("\n")
}
e = pf.efun($1, &$&2)
if (!stoprun) {
if (minerr == -1 || (($3 != 1) && (e < minerr))) {
minerr = e
if (saveflag) {
ttime = startsw() - st
savepath.aopen(savefit_name)
lines += 1
if( verbose ) {
if( $3 == 1 ) { savepath.printf("C") }
savepath.printf("EF \t")
}
savepath.printf("%g %d %d %-12.8g ", temptr, nefun, ttime, minerr)
saveval(savepath)
savepath.close
if( verbose ) {
verbfile.printf("WROTE new_point ")
new_point.printf(verbfile)
verbfile.printf("vs. sa_efun sent ")
for i = 0,$1-1 {
verbfile.printf("%-12.8g ",$&2[i])
}
verbfile.printf("\n")
verbfile.printf("vs. pf.pmlst point ")
for i = 0,pf.parmlist.count-1 {
verbfile.printf("%-12.8g ",pf.parmlist.object(i).val)
}
verbfile.printf("\n")
}
}
}
}
doNotify()
if (verbose) {
verbfile.printf("%d %14.12g ",nefun, e)
for i = 0,$1-1 {
verbfile.printf("%14.12g ",$&2[i])
}
verbfile.printf("\n")
}
enow = e
return e
}
/****************************************************************
Use the data stored in savepath.rst file to initialize
the simplex, and set the starting temperature.
****************************************************************/
func init_from_restart_file() { local i, j, rst_temp, n_par, nrow, tmpval, idx
rest_file.ropen(rest_inname)
rst_temp = rest_file.scanvar()
n_par = rest_file.scanvar()
nrow = rest_file.scanvar()
tmp_vec = new Vector(n_par)
tl = pf.parmlist
/*** first, read which variables were optimized in the last run. ***/
for i = 0, n_par-1 {
tmp_vec.x[i] = rest_file.scanvar()
if( tl.object(i).doarg != tmp_vec.x[i] ) {
printf("Parameters optimized in restart file %s do not \
match those specified here. Abort init_restart\n", \
rest_inname)
rest_file.close()
if( verbose ) {
verbfile.printf("Parameters optimized in restart file %s do not \
match those specified here. Abort init_restart\n", \
rest_inname)
}
return 0
}
}
/*** now, read the actual values of each of the simplex points. ***/
for j = 0, nrow-1 {
idx = 0
for i = 0, n_par-1 {
tmpval = rest_file.scanvar()
if( ( tmp_vec.x[i] == 0 ) && ((tl.object(i).val - tmpval) > ftol) ) {
printf("Value of %s = %g in restart file, != current value %g\n",\
tl.object(i).name,tmpval,tl.object(i).val)
if( verbose ) {
verbfile.printf("Value of %s = %g in restart file, != current value %g\n",\
tl.object(i).name,tmpval,tl.object(i).val)
}
rest_file.close()
return 0
} else {
if( tmp_vec.x[i] ) {
simplex.x[j][idx] = code_value(tmpval,idx)
idx += 1
}
}
}
}
rest_file.close()
start_temptr = rst_temp
printf("Restart file %s read; Begin at temperature %g\n",rest_inname,rst_temp)
printf("Simplex initialized to: \n")
simplex.printf
if( verbose ) {
verbfile.printf("Restart file %s read; Begin at temperature %g\n",rest_inname,rst_temp)
simplex.fprintf(verbfile)
}
return 1
}
/**************************************************************
Initialize the simplex, define best_point and mean_point,
and evaluate the error function at simplex points.
**************************************************************/
proc initialize_simplex() { local rval, inbound, abval, i,j,bdval, did_rst
abval = 0
for i = 0, new_point.size-1 { abval += abs(true_value(new_point.x[i],i)) }
// initialise the simplex
nparam = new_point.size()
npoints = nparam + 1
simplex = new Matrix(npoints,nparam)
did_rst = 0
if( use_rst ) {
// Use the data stored in savepath.rst file to initialize
// the simplex, and set the starting temperature.
did_rst = init_from_restart_file()
}
// Either no restart was specified, or there was an error reading the
// restart file; initialize the simplex randomly.
//
// NOTE: cmw, jun 06: This initialization does NOT work well when the
// points are converted to log space. Initialize the points in linear
// space first, then convert them (if appropriate) to log space before
// they are added to the simplex.
//
if( !use_rst || !did_rst ) {
for j = 0, nparam-1 {
if( (pf.parmlist_use.object(j).low == -1e6) || \
(pf.parmlist_use.object(j).high == 1e6) ) {
bdval = 2*abval
} else {
bdval = Kbound*(pf.parmlist_use.object(j).high-pf.parmlist_use.object(j).low)
}
for i = 0,npoints-1 {
inbound = 0
while( !inbound ) {
rval = 0.5-r.repick()
simplex.x[i][j] = code_value(true_value(new_point.x[j],j) + rval*bdval,j)
if( (simplex.x[i][j] >= code_low.x[j] ) && \
(simplex.x[i][j] <= code_high.x[j]) ) {
inbound = 1
if( verbose ) {
verbfile.printf("Set row %d, %s to %g\n",j, \
pf.parmlist_use.object(j).name,\
true_value(simplex.x[i][j],j))
}
} else {
if( verbose ) {
verbfile.printf("Initialized row %d, %s to %g, out of range [%g, %g]\n",\
j,pf.parmlist_use.object(j).name,\
true_value(simplex.x[i][j],j),\
pf.parmlist_use.object(j).low,pf.parmlist_use.object(j).high)
}
}
}
}
}
}
best_point = new Vector(nparam)
mean_point = new Vector(nparam)
centrd = new Vector(nparam)
// initialise the error vector
simplex_errors = new Vector(npoints)
for i = 0,npoints-1 {
if (stoprun) return
new_point = simplex.getrow(i)
if( verbose ) {
verbfile.printf("INIT (OOB %d) new_point is ",out_of_bounds(new_point))
new_point.printf(verbfile,"%-12.8 ")
}
simplex_errors.x[i] = sa_efun(new_point.size(),&new_point.x[0],out_of_bounds(new_point))
}
//initialize the vector which stores the applied noise.
simplex_noise = new Vector(npoints,0)
}
func prun() { local lines
minerr = -1
restart_temptr = restart_factor * start_temptr
if (saveflag) {
sprint(savefit_name,"%s.tmp",savepath_name)
sprint(annealfit_name,"%s.ann",savepath_name)
lines = 0
}
if( verbose ) {
printf("Opening %s.vrb\n",savepath_name)
sprint(verbname,"%s.vrb",savepath_name)
verbfile.wopen(verbname)
}
// flag: whether the simplex has 'collapsed', and the trial number at which this occurs
simplx_stop = -1
simplx_stoptrial = 0
pf.doarg_get(new_point)
set_coded_bounds()
nefun = 0
eb = 1e9
st = startsw()
time = st
if (new_point.size == 0) {
minerr = pf.efun(0, &time) // time is dummy here
} else {
if (stoprun) {return minerr}
minerr = fit_simanneal_seq()
new_point = best_point.c
if( verbose ) {
verbfile.printf("PRUN before sa_efun, new_point = ")
new_point.printf(verbfile)
}
minerr = sa_efun(new_point.size(),&new_point.x[0],out_of_bounds(new_point))
}
time = startsw() - st
if (saveflag) {
if( savepath.isopen ) { savepath.close }
// write the final temp & params to the .tmp file, then copy data over
// to the .fit file
ttime = startsw() - st
savepath.aopen(savefit_name)
savepath.printf("%g %d %d %-12.8g ", temptr, nefun, ttime, minerr)
save_best_params(savepath)
savepath.close
sprint(savefit_name,"%s.fit",savepath_name)
savepath.aopen(savefit_name)
savepath.printf("start\n")
// First line of .fit file records the number of columns in each line of the data.
savepath.printf("%d %d\n", lines, pf.parmlist.count + 4)
savepath.close
sprint(cmnd,"cat %s.tmp >> %s.fit",savepath_name,savepath_name)
system(cmnd)
if( annealpath.isopen ) { annealpath.close }
sprint(annealfit_name,"%s.annfit",savepath_name)
annealpath.aopen(annealfit_name)
annealpath.printf("start\n")
annealpath.printf("%d %d\n", lines, pf.parmlist.count + 4)
annealpath.close
sprint(cmnd,"cat %s.ann >> %s.annfit",savepath_name,savepath_name)
system(cmnd)
print_info_file()
}
printf("Final temperature: %9.5g\tnefun: %d\tFinal error: %9.5g\tat point ",temptr,nefun,eb)
best_point.printf("%9.4g")
if( verbose ) {
verbfile.printf("Final temperature: %9.5g\tnefun: %d\tFinal error: %9.5g\tat point ",temptr,nefun,eb)
best_point.printf(verbfile)
verbfile.close
}
return minerr
}
proc showopt() {
variable_domain(&start_temptr,0,1e12)
variable_domain(&cool_rate,0,1)
variable_domain(&iter_per_step,1,1e12)
xpanel("")
xbutton("Estimate initial temp","estimate_inittemp()")
xcheckbox("Verbose output",&verbose)
xcheckbox("Restart from .rst file",&use_rst)
xvarlabel(rest_inname)
xbutton("Change restart filename","get_rstname(rest_inname)")
xlabel("Simplex parameters")
xpvalue("alpha",&alpha,1)
xpvalue("beta",&beta,1)
xpvalue("gamma",&gamma,1)
xpvalue("lambda",&lambda,1)
xlabel("Simulated annealing parameters")
xmenu("Select cooling schedule")
xbutton("Proportional decrease","set_annsched(0)")
xmenu()
xvarlabel(ann_str)
xpvalue("Cooling rate",&cool_rate,1)
xpvalue("Start temperature",&start_temptr,1)
xpvalue("Iterations per cooling step",&iter_per_step,1)
xpvalue("ftol",&ftol,1)
xpvalue("Norm grad tol",>ol,1)
xpvalue("Restart factor",&restart_factor,1)
xpanel()
xpanel("")
xpvalue("Current temperature",&temptr)
xpvalue("Current lowest error",&elo)
xpvalue("Current tolerance",&rtol)
xpanel()
}
proc set_annsched() {
ann_sched = $1
if( ann_sched == 0 ) {
sprint(ann_str,"Proportional decrease")
}
}
func noisy_simplex() { local i, e, iter, esave, gstop
iter = $1
Ni = pf.parmlist_use.count * 5
Ci = Ci_avg = CiM1_avg = 0
m1 = m2 = delf = 0
nchk = Ni
ndone = 1
gstop = 0
while (1) {
if( verbose ) {
verbfile.printf("\n\nTop of loop, ")
print_simplex_info(verbfile)
}
if (stoprun) break
// find the highest (worst), next highest and lowest (best) points
//
// whenever we "look" at a vertex, it gets a random thermal fluctuation.
// simplex_errors stores the error of the simplex points,
// plus these fluctuations.
//
simperrors_cpy = simplex_errors.c.sub(simplex_noise)
simplex_errors = simperrors_cpy.c
for i = 0, npoints-1 {
// By storing the noise applied at each step, in addition to the
// noise-amplified error, we are able to retain the actual function
// values observed, which is crucial to the Press algorithm. Otherwise
// we keep adding noise on top of itself, which is a real problem when the
// noise is large!!
//
if( verbose ) verbfile.printf("\tError[%d] E0 %g\t",i,simplex_errors.x[i])
simplex_noise.x[i] = -1*temptr*log(r.repick())
simplex_errors.x[i] += simplex_noise.x[i]
if( verbose ) verbfile.printf("N %g\tEN %g\n",simplex_noise.x[i],simplex_errors.x[i])
}
ilo = simplex_errors.min_ind() // index of smallest error (incl fluct'n)
ihi = simplex_errors.max_ind()
ehi = simplex_errors.x[ihi]
simplex_errors.x[ihi] = -1
inhi = simplex_errors.max_ind()
simplex_errors.x[ihi] = ehi
elo = simplex_errors.x[ilo]
enhi = simplex_errors.x[inhi]
if( verbose ) {
verbfile.printf("Worst point ")
for i = 0, simplex.ncol-1 {
verbfile.printf("%g ",simplex.getrow(ihi).x[i])
}
verbfile.printf("E0 %g\tN %g EN %g\n",simplex_errors.x[ihi]-simplex_noise.x[ihi], \
simplex_noise.x[ihi],simplex_errors.x[ihi])
}
//
// track 'uphill' steps.
//
// error values without random fluct'n
ehi_det = simperrors_cpy.x[ihi]
elo_det = simperrors_cpy.x[ilo]
enhi_det = simperrors_cpy.x[inhi]
// test for completion
rtol = 2 * abs(ehi - elo)/( abs(ehi) + abs(elo) )
if (rtol < ftol || iter < 0) {
break
}
// Compute vector average of all points except highest
mean_point.fill(0)
for i = 0, npoints-1 {
if (i != ihi) mean_point.add(simplex.getrow(i))
}
mean_point.div(npoints-1)
iter -= 2
old_simplex = simplex.c()
if (stoprun) break
esave = ehi
// Begin a new iteration. First extrapolate by alpha through the face of
// the simplex across from the high point, i.e. reflect the simplex from the
// high point
e = new_step(alpha)
Ci += eval_centroid()
ndone += 1
if( ndone >= nchk ) {
gstop = estimate_normgrad()
if( gstop ) {
printf("Normalized gradient shows little change. Stop annealing.\n")
break
}
}
if (e < elo) { // if lower than the lowest... N.B Press et al use <= but this seems to lead to loops
// Gives a better result than the best point, so try an additional
// extrapolation by a factor of beta.
if (stoprun) break
if (verbose) verbfile.printf("New best point. Try a further extrapolation.\n")
e = new_step(beta) // ...Extrapolate in the same direction
// add Ci calculation here.
if( verbose ) print_simplex_info(verbfile)
Ci += eval_centroid()
ndone += 1
if( ndone >= nchk ) {
gstop = estimate_normgrad()
if( gstop ) {
printf("Normalized gradient shows little change. Stop annealing.\n")
break
}
}
e_det = simplex_errors.x[ihi]-simplex_noise.x[ihi]
if( e_det < ehi_det ) {
m1 += 1 // increment count of successful moves
if( verbose ) verbfile.printf("increment m1 %d\n",m1)
} else {
m2 += 1 // increment count of unsuccessful moves
if( verbose ) {
verbfile.printf("increment m2 %d curr_delf %g - %g = %g tot_delf ",\
m2,e_det,ehi_det,e_det-ehi_det)
}
delf += (e_det - ehi_det)
if( verbose ) { verbfile.printf("%g\n",delf) }
}
} else if (e >= enhi) { // else if higher than the second-highest...
// The reflected point is worse than the second-highest, so look for an
// intermediate lower point, i.e. do a 1-D contraction
//
// Contraction is either positive or negative: cf.
// Wright MH. Direct search methods: once scorned, now respectable.
// In Griffiths DF and Watson GA (eds.), Numerical Analysis 1995
// (Proceedings of the 1995 Dundee Biennial Conference in Numerical
// Analysis), 191-208, Addison Wesley Longman, Harlow, U.K., 1996.
if (stoprun) break
if (ehi < esave) { // new point is better than the worst point;
// do an outside (+) contraction
if (verbose) verbfile.printf("Do a one dimensional contraction (+)\n")
esave = ehi
e = new_step(gamma) // ...Do a 1D positive contraction\
} else { // new point is worse than the worst point in
// simplex; do an inside (-) contraction
if (verbose) verbfile.printf("Do a one dimensional contraction (-)\n")
e = new_step(-1*gamma) // ...Do a 1D negative contraction\
}
if (e >= esave) { // if the highest point is still there...
if (verbose) verbfile.printf("Contract around the best point\n")
for i = 0, npoints-1 { // ...contract about the lowest point
if (stoprun) break
if (i != ilo) {
new_point = simplex.getrow(i).add(simplex.getrow(ilo)).mul(gamma)
simplex.setrow(i,new_point)
if( verbose ) {
verbfile.printf("E>=ESAVE, before sa_efun new_point = ")
new_point.printf(verbfile)
}
simplex_errors.x[i] = sa_efun(new_point.size(),&new_point.x[0],out_of_bounds(new_point))
simplex_noise.x[i] = 0
if( simplex_errors.x[i] < simperrors_cpy.x[i] ) {
m1 += 1
if( verbose ) verbfile.printf("increment m1 %d\n",m1)
} else {
m2 += 1
if( verbose ) {
verbfile.printf("increment m2 %d curr_delf %g - %g = %g tot_delf ",\
m2,simplex_errors.x[i],simperrors_cpy.x[i],\
simplex_errors.x[i] - simperrors_cpy.x[i])
}
delf += (simplex_errors.x[i] - simperrors_cpy.x[i])
if( verbose ) { verbfile.printf("%g\n",delf) }
}
}
}
//Calculate Ci here.
Ci += eval_centroid()
ndone += 1
if( ndone >= nchk ) {
gstop = estimate_normgrad()
if( gstop ) {
printf("Normalized gradient shows little change. Stop annealing.\n")
break
}
}
iter -= npoints-1
} else {
//Calculate Ci here.
Ci += eval_centroid()
ndone += 1
if( ndone >= nchk ) {
gstop = estimate_normgrad()
if( gstop ) break
}
if( e_det < ehi_det ) {
m1 += 1 // increment count of successful moves
if( verbose ) verbfile.printf("increment m1 %d\n",m1)
} else {
m2 += 1 // increment count of unsuccessful moves
if( verbose ) {
verbfile.printf("increment m2 %d curr_delf %g - %g = %g tot_delf ",\
m2,e_det,ehi_det,e_det - ehi_det)
}
delf += (e_det - ehi_det)
if( verbose ) { verbfile.printf("%g\n",delf) }
}
}
} else {
//Calculate Ci here.
Ci += eval_centroid()
ndone += 1
if( ndone >= nchk ) {
gstop = estimate_normgrad()
if( gstop ) break
}
if( verbose ) verbfile.printf("Improved over the worst points. Accepted. \n")
if( e_det < ehi_det ) {
m1 += 1 // increment count of successful moves
if( verbose ) verbfile.printf("increment m1 %d\n",m1)
} else {
m2 += 1 // increment count of unsuccessful moves
if( verbose ) {
verbfile.printf("increment m2 %d curr_delf %g - %g = %g tot_delf ",\
m2,e_det,ehi_det,e_det - ehi_det)
}
delf += (e_det - ehi_det)
if( verbose ) { verbfile.printf("%g\n",delf) }
}
iter += 1
}
}
return iter
}
/*****************************************************************
Same as noisy_simplex above, except it calculates m1 and m2
for initial temperature estimation.
*****************************************************************/
func noisy_simplex_m1m2() { local i, e, iter, esave
iter = $1
m1 = m2 = delf = 0
while (1) {
printf("m1 %d + m2 %d = %d < %d\n",m1,m2,m1+m2,iter_per_step)
verbfile.printf("m1 %d + m2 %d = %d < %d\n",m1,m2,m1+m2,iter_per_step)
if( m1+m2 > iter_per_step ) break
// find the highest (worst), next highest and lowest (best) points
//
// whenever we "look" at a vertex, it gets a random thermal fluctuation.
// simplex_errors stores the error of the simplex points,
// plus these fluctuations.
//
simperrors_cpy = simplex_errors.c.sub(simplex_noise)
simplex_errors = simperrors_cpy.c
for i = 0, npoints-1 {
// By storing the noise applied at each step, in addition to the
// noise-amplified error, we are able to retain the actual function
// values observed, which is crucial to the Press algorithm. Otherwise
// we keep adding noise on top of itself, which is a real problem when the
// noise is large!!
//
if( verbose ) {
verbfile.printf("\tError[%d] E0 %g\t",i,simplex_errors.x[i])
}
simplex_noise.x[i] = -1*temptr*log(r.repick())
simplex_errors.x[i] += simplex_noise.x[i]
if( verbose ) {
verbfile.printf("N %g\tEN %g\n",simplex_noise.x[i],simplex_errors.x[i])
}
}
ilo = simplex_errors.min_ind() // index of smallest error (incl fluct'n)
ihi = simplex_errors.max_ind()
ehi = simplex_errors.x[ihi]
simplex_errors.x[ihi] = -1
inhi = simplex_errors.max_ind()
simplex_errors.x[ihi] = ehi
elo = simplex_errors.x[ilo]
enhi = simplex_errors.x[inhi]
if( verbose ) {
verbfile.printf("Worst point ")
for i = 0, simplex.ncol-1 {
verbfile.printf("%g ",simplex.getrow(ihi).x[i])
}
verbfile.printf("E0 %g\tN %g EN %g\n",simplex_errors.x[ihi]-simplex_noise.x[ihi], \
simplex_noise.x[ihi],simplex_errors.x[ihi])
}
//
// track 'uphill' steps.
//
// error values without random fluct'n
ehi_det = simperrors_cpy.x[ihi]
elo_det = simperrors_cpy.x[ilo]
enhi_det = simperrors_cpy.x[inhi]
// test for completion
rtol = 2 * abs(ehi - elo)/( abs(ehi) + abs(elo) )
// Compute vector average of all points except highest
mean_point.fill(0)
for i = 0, npoints-1 {
if (i != ihi) mean_point.add(simplex.getrow(i))
}
mean_point.div(npoints-1)
iter -= 2
old_simplex = simplex.c()
if( m1+m2 > iter_per_step ) break
esave = ehi
// Begin a new iteration. First extrapolate by alpha through the face of
// the simplex across from the high point, i.e. reflect the simplex from the
// high point
e = new_step(alpha)
if (e < elo) { // if lower than the lowest... N.B Press et al use <= but this seems to lead to loops
// Gives a better result than the best point, so try an additional
// extrapolation by a factor of beta.
if( m1+m2 > iter_per_step ) break
if (verbose) verbfile.printf("New best point. Try a further extrapolation.\n")
e = new_step(beta) // ...Extrapolate in the same direction
if( verbose ) print_simplex_info(verbfile)
e_det = simplex_errors.x[ihi]-simplex_noise.x[ihi]
if( e_det < ehi_det ) {
m1 += 1 // increment count of successful moves
if( verbose ) verbfile.printf("increment m1 %d\n",m1)
} else {
m2 += 1 // increment count of unsuccessful moves
if( verbose ) {
verbfile.printf("increment m2 %d curr_delf %g - %g = %g tot_delf ",\
m2,e_det,ehi_det,e_det-ehi_det)
}
delf += (e_det - ehi_det)
if( verbose ) { verbfile.printf("%g\n",delf) }
}
} else if (e >= enhi) { // else if higher than the second-highest...
// The reflected point is worse than the second-highest, so look for an
// intermediate lower point, i.e. do a 1-D contraction
//
// Contraction is either positive or negative: cf.
// Wright MH. Direct search methods: once scorned, now respectable.
// In Griffiths DF and Watson GA (eds.), Numerical Analysis 1995
// (Proceedings of the 1995 Dundee Biennial Conference in Numerical
// Analysis), 191-208, Addison Wesley Longman, Harlow, U.K., 1996.
if( m1+m2 > iter_per_step ) break
if (ehi < esave) { // new point is better than the worst point;
// do an outside (+) contraction
if (verbose) verbfile.printf("Do a one dimensional contraction (+)\n")
esave = ehi
e = new_step(gamma) // ...Do a 1D positive contraction\
} else { // new point is worse than the worst point in
// simplex; do an inside (-) contraction
if (verbose) verbfile.printf("Do a one dimensional contraction (-)\n")
e = new_step(-1*gamma) // ...Do a 1D negative contraction\
}
if (e >= esave) { // if the highest point is still there...
if (verbose) verbfile.printf("Contract around the best point\n")
for i = 0, npoints-1 { // ...contract about the lowest point
if( m1+m2 > iter_per_step ) break
if (i != ilo) {
new_point = simplex.getrow(i).add(simplex.getrow(ilo)).mul(gamma)
simplex.setrow(i,new_point)
if( verbose ) {
verbfile.printf("M1M2: E>=ESAVE, before sa_efun new_point = ")
new_point.printf(verbfile)
}
simplex_errors.x[i] = sa_efun(new_point.size(),&new_point.x[0],out_of_bounds(new_point))
simplex_noise.x[i] = 0
if( simplex_errors.x[i] < simperrors_cpy.x[i] ) {
m1 += 1
if( verbose ) verbfile.printf("increment m1 %d\n",m1)
} else {
m2 += 1
if( verbose ) {
verbfile.printf("increment m2 %d curr_delf %g - %g = %g tot_delf ",\
m2,simplex_errors.x[i],simperrors_cpy.x[i],\
simplex_errors.x[i] - simperrors_cpy.x[i])
}
delf += (simplex_errors.x[i] - simperrors_cpy.x[i])
if( verbose ) { verbfile.printf("%g\n",delf) }
}
}
}
iter -= npoints-1
} else {
if( e_det < ehi_det ) {
m1 += 1 // increment count of successful moves
if( verbose ) verbfile.printf("increment m1 %d\n",m1)
} else {
m2 += 1 // increment count of unsuccessful moves
if( verbose ) {
verbfile.printf("increment m2 %d curr_delf %g - %g = %g tot_delf ",\
m2,e_det,ehi_det,e_det - ehi_det)
}
delf += (e_det - ehi_det)
if( verbose ) { verbfile.printf("%g\n",delf) }
}
}
} else {
if( verbose ) verbfile.printf("Improved over the worst points. Accepted. \n")
if( e_det < ehi_det ) {
m1 += 1 // increment count of successful moves
if( verbose ) verbfile.printf("increment m1 %d\n",m1)
} else {
m2 += 1 // increment count of unsuccessful moves
if( verbose ) {
verbfile.printf("increment m2 %d curr_delf %g - %g = %g tot_delf ",\
m2,e_det,ehi_det,e_det - ehi_det)
}
delf += (e_det - ehi_det)
if( verbose ) { verbfile.printf("%g\n",delf) }
}
iter += 1
}
}
if( m2 > 0 ) delf /= m2
return iter
}
func new_step() { local e, ttime, fluc, i
new_point = mean_point.c.mul(1+$1).sub(old_simplex.getrow(ihi).mul($1))
if( out_of_bounds(new_point) ) {
e = bd_pnlty
if( wrapfile.isopen ) {
wrapfile.printf("Point was ")
new_point.printf(wrapfile)
}
// if the first attempt is out-of-bounds, recenter the
// point around the best point of the simplex so far.
if( verbose ) {
verbfile.printf("Before recenter new_point = ")
new_point.printf(verbfile)
}
recenter_point()
if( verbose ) {
verbfile.printf("After recenter new_point = ")
new_point.printf(verbfile)
}
}
if( verbose ) {
verbfile.printf("NEW_STEP, before sa_efun new_point = ")
new_point.printf(verbfile)
}
e = sa_efun(new_point.size(),&new_point.x[0],out_of_bounds(new_point))
if( out_of_bounds(new_point) ) {
// if point is still out of bounds, even after recentering,
// penalize it.
e += bd_pnlty
}
e_det = e
enow = e
if (e <= eb) { // Save the best-ever point
if( verbose ) {
verbfile.printf("Best point was = ")
best_point.printf(verbfile)
}
best_point = new_point
eb = e
if( verbose ) {
verbfile.printf("Best point NOW = ")
best_point.printf(verbfile)
}
}
fluc = -1*temptr*log(r.repick())
e -= fluc
if( verbose ) {
verbfile.printf("\tNew point proposed ")
for i = 0, simplex.ncol-1 {
verbfile.printf("%g ",new_point.x[i])
}
verbfile.printf("E0 %g\tN %g EN %g\n",e+fluc,fluc,e)
}
if (e < ehi) {
if (verbose) {
verbfile.printf("Replaced highest point\t")
if( e_det < ehi_det ) {
verbfile.printf("true downhill point %g (%g) replaces %g (%g)\n",e,e_det,ehi,ehi_det)
} else { verbfile.printf("an uphill point added to simplex %g (%g) replaces %g (%g)\n",e,e_det,ehi,ehi_det) }
}
ehi = e
if( verbose ) {
verbfile.printf("\tPoint was ")
simplex.getrow(ihi).printf(verbfile)
verbfile.printf("\tNow is ")
new_point.printf(verbfile)
}
simplex.setrow(ihi,new_point)
simplex_errors.x[ihi] = e
simplex_noise.x[ihi] = -1*fluc
if( verbose ) {
verbfile.printf("After replacement, ")
print_simplex_info(verbfile)
}
}
return e
}
proc recenter_point() { local j, abval, bdval, rval, tmplo
// Center about the best point, WITHOUT the noise.
tmplo = simplex_errors.c.sub(simplex_noise).min_ind()
if( verbose ) {
verbfile.printf("Recentering out-of-bounds point ")
for j = 0, new_point.size-1 {
verbfile.printf(" %g",new_point.x[j])
}
verbfile.printf("\nabout best point %d with error %g (%g) (K=%g)\n",tmplo,simplex_errors.x[tmplo]-simplex_noise.x[tmplo],simplex_errors.x[tmplo],Kbound)
for j = 0, new_point.size-1 {
verbfile.printf(" %g",simplex.getrow(tmplo).x[j])
}
verbfile.printf("\n")
}
for j = 0, nparam-1 {
if( (pf.parmlist_use.object(j).low <= -1e6) || \
(pf.parmlist_use.object(j).high >= 1e6) ) {
bdval = 2*simplex.getrow(tmplo).abs.x[j]
} else {
bdval = Kbound*(code_high.x[j]-code_low.x[j])
}
rval = 0.5-r.repick()
new_point.x[j] = simplex.x[tmplo][j] + rval*bdval
}
if( verbose ) {
verbfile.printf("\nNow using ")
for j = 0, new_point.size-1 {
verbfile.printf(" %g",new_point.x[j])
}
verbfile.printf("\n")
}
}
proc estimate_inittemp() { local lines
minerr = -1
restart_temptr = restart_factor * start_temptr
sprint(verbname,"%s.est",savepath_name)
if (saveflag) {
sprint(savefit_name,"%s.estfit",savepath_name)
lines = 0
}
verbfile.wopen(verbname)
pf.doarg_get(new_point)
set_coded_bounds()
nefun = 0
eb = 1e9
st = startsw()
time = st
if (new_point.size == 0) {
minerr = pf.efun(0, &time) // time is dummy here
} else {
if (stoprun) {return minerr}
// Estimate the initial starting temp as in Cardoso et al. 1996
pf.doarg_get(new_point)
minerr = sa_efun(new_point.size(),&new_point.x[0],out_of_bounds(new_point))
start_temptr = 1e5*minerr
printf("Initial error %g; start temp at %g\n",minerr,start_temptr)
verbfile.printf("Initial error %g; start temp at %g\n",minerr,start_temptr)
verbose = 0
start_temptr = temp_est()
new_point = best_point.c
minerr = sa_efun(new_point.size(),&new_point.x[0],out_of_bounds(new_point))
}
time = startsw() - st
if( savepath.isopen ) { savepath.close }
best_point.printf("%9.4g")
verbfile.printf("Best point is ")
best_point.printf(verbfile)
verbfile.close
}
proc temp_est() { local i, id, rval
wrapfile = new File()
printf("Estimating initial temperature\n")
if( verbfile.isopen ) verbfile.printf("Estimating initial temperature\n")
initialize_simplex()
//Start the optimisation
iter = -1
temptr = start_temptr
//iter_per_step = pf.parmlist_use.count * 100
printf("Try starting temp of %g; Do %d iterations\n",start_temptr,iter_per_step)
if( verbfile.isopen ) {
verbfile.printf("Try starting temp of %g; Do %d iterations\n",start_temptr,iter_per_step)
}
// Only do one iteration of the outer loop, i.e. iter_per_step function iterations
iter = iter_per_step
iter = noisy_simplex_m1m2(iter)
printf("\nSimplex complete.\n")
printf("Estimating initial temperature...\nm1 = %g\tm2 = %g\tdel F = %g\n",m1,m2,delf)
if( verbfile.isopen ) {
verbfile.printf("\nSimplex complete.\n")
verbfile.printf("Estimating initial temperature...\nm1 = %g\tm2 = %g\tdel F = %g\n",m1,m2,delf)
}
chitmp = 0.95
if( m2 > 0 ) { tmpval = (m2*chitmp - m1*(1-chitmp))/m2 } else { tmpval = 1e6 }
tempest = -delf / log(tmpval)
printf("Initial temperature estimate is %g\n",tempest)
printf("Ended noisy_simplex_m1m2\n")
if( verbfile.isopen ) {
verbfile.printf("Initial temperature estimate is %g\n",tempest)
verbfile.printf("Ended noisy_simplex_m1m2\n")
}
}
proc fit_simanneal_seq() { local i, id, rval
wrapfile = new File()
initialize_simplex()
//Start the optimisation
iter = -1
temptr = start_temptr
if( saveflag ) {
sprint(rest_name, "%s.rst",savepath_name)
write_restart_file()
}
while (iter < 0 && !stoprun && temptr > quit_temptr) {
iter = iter_per_step
iter = noisy_simplex(iter)
if( verbose ) {
verbfile.printf("Temp = %g, m1 %d m2 %d delf %g\n",temptr,m1,m2,delf)
}
if (stoprun) return
if( iter < 0 && simplx_stop < 0 ) {
simplx_stop = 1
simplx_stoptrial = nefun
simplx_temptr = temptr
}
anneal_step(ann_sched)
}
if( wrapfile.isopen ) { wrapfile.close }
}
/*****************************************************
Apply the annealing schedule
*****************************************************/
proc anneal_step() { local newTemptr, i, in_simplex, oldTemptr
time = startsw() - st
if (saveflag) {
write_restart_file()
// this is the beginning of the annealing step. Write the current
// best to the savepath.fit file.
savepath.aopen(savefit_name)
lines += 1
if( verbose ) { savepath.printf("AS\t") }
savepath.printf("%g %d %d %-12.8g ", temptr, nefun, ttime, minerr)
save_best_params(savepath)
savepath.close
// best_point doesn't give the actual *values* of the parameters being
// fitted; use saveval
//
// This writes the CURRENT error and parameter set, at the start of
// this annealing step.
annealpath.aopen(annealfit_name)
annealpath.printf("%g %d %d %-12.8g ", temptr, nefun, time, enow)
saveval(annealpath)
annealpath.close
}
printf("temptr=%g nefun=%d time=%d eb=%-12.8g\n", temptr, nefun, time, eb)
oldTemptr = temptr
if( $1 == 1 ) { // more complex schedule
newTemptr = 1.0*(elo-eb)
if (newTemptr < temptr) {
if (newTemptr > 0.5*temptr) {
temptr = newTemptr
} else {
temptr = 0.5*temptr
}
}
} else { // simplest annealing schedule
temptr *= cool_rate
}
// update Kbound, the compression parameter for bound constraints
anneal_count += 1
if( anneal_count/2 == int(anneal_count/2) ) { // an even annealing step
Kbound = Kbound_m2 * temptr / oldTemptr
if( verbose ) {
printf("Updating Kbound compression, anneal step %d\n",anneal_count)
printf("Kbound was 1, now Kbound = %g * %g / %g = %g\n", Kbound_m2,temptr,oldTemptr,Kbound)
verbfile.printf("Updating Kbound compression, anneal step %d\n",anneal_count)
verbfile.printf("Kbound was 1, now Kbound = %g * %g / %g = %g\n", Kbound_m2,temptr,oldTemptr,Kbound)
}
} else {
Kbound_m2 = Kbound
Kbound = 1
if( verbose ) {
printf("Updating Kbound compression, anneal step %d\n",anneal_count)
printf("Kbound was %g, now Kbound = %g\n",Kbound_m2,Kbound)
verbfile.printf("Updating Kbound compression, anneal step %d\n",anneal_count)
verbfile.printf("Kbound was %g, now Kbound = %g\n",Kbound_m2,Kbound)
}
}
if (temptr < restart_temptr) {
restart_temptr *= restart_factor
in_simplex = 0
for i = 0, npoints-1 {
in_simplex += simplex.getrow(i).eq(best_point)
}
if (in_simplex > 0) {
print "Best point is in simplex. No restart."
} else {
simplex.setrow(ihi,best_point)
printf("Restarting at temperature %g\n",temptr)
}
}
}
/***
The vector "best_point" only stores the parameters used in the search. To
get all parameters, you must retrieve data from pf.parmlist
***/
proc save_best_params() { local i, use_idx
// $o1 is the file to write
use_idx = 0
tl = pf.parmlist
printf("Writing best_point = ")
for i=0, tl.count-1 {
if( tl.object(i).doarg ) {
$o1.printf("%-12.8g ",true_value(best_point.x[use_idx],use_idx))
printf("%-12.8g ",true_value(best_point.x[use_idx],use_idx))
use_idx += 1
} else {
$o1.printf("%-12.8g ", tl.object(i).val)
printf("%-12.8g ", tl.object(i).val)
//$o1.printf("%-12.8g ", tl.object(i).val.x[0])
//printf("%-12.8g ", tl.object(i).val.x[0])
}
}
$o1.printf("\n")
printf("\n")
}
proc saveval() { local i
// $o1 is the file to write
if (numarg() == 2) $o1.printf("%s", $s2)
tl = pf.parmlist
for i=0, tl.count-1 {
$o1.printf("%-12.8g ", tl.object(i).val)
//$o1.printf("%-12.8g ", tl.object(i).val.x[0])
}
$o1.printf("\n")
}
func out_of_bounds() { local i
for i = 0, $o1.size-1 {
rval = $o1.x[i]
if( rval < code_low.x[i] || rval > code_high.x[i] ) {
rval = true_value($o1.x[i],i)
printf("Param %s = %g (%g) out of range [%g, %g]; adjust point\n",\
pf.parmlist_use.object(i).name,$o1.x[i],rval,\
pf.parmlist_use.object(i).low,pf.parmlist_use.object(i).high)
if( verbose ) verbfile.printf("Param %s (%g) out of range [%g, %g]; adjust point\n",\
pf.parmlist_use.object(i).name,rval,\
pf.parmlist_use.object(i).low,pf.parmlist_use.object(i).high)
return 1
}
}
return 0
}
func true_value() {
if (pf.parmlist_use.object($2).uselog) {
if (pf.parmlist_use.object($2).low > 0) {
return exp($1)
}else if (pf.parmlist_use.object($2).high < 0) {
return exp(-$1)
}
return $1
}else{
return $1
}
}
func code_value() {
if (pf.parmlist_use.object($2).uselog) {
if (pf.parmlist_use.object($2).low > 0) {
return log($1)
}else if (pf.parmlist_use.object($2).high < 0) {
return log(-$1)
}
return $1
}else{
return $1
}
}
/**********************************************************
write_point2file()
writes the specified point (vector) $o1 to the file $o2.
if the point is logscale it is converted back to its
true value.
**********************************************************/
proc write_point2file() { local i
// $o2 is the file to write
// $o1 is the point to write
for i=0, $o1.size-1 {
$o2.printf("%-12.8g ", true_value($o1.x[i],i))
}
$o2.printf("\n")
}
proc print_simplex_info() { local i, j
$o1.printf("Simplex info:\n")
for i = 0, simplex.nrow-1 {
$o1.printf("[%d]\t",i)
for j = 0, simplex.ncol-1 {
$o1.printf("%g ",simplex.getrow(i).x[j])
}
$o1.printf("E0 %g\tN %g\tEN %g\n",simplex_errors.x[i]-simplex_noise.x[i],\
simplex_noise.x[i],simplex_errors.x[i])
}
}
func eval_centroid() { local i, cent_err
centrd.fill(0)
for i = 0, npoints-1 {
centrd.add(simplex.getrow(i))
}
centrd.div(npoints)
if( verbose ) {
verbfile.printf("CNTRD, before sa_efun new_point = ")
centrd.printf(verbfile)
}
cent_err = sa_efun(centrd.size(),¢rd.x[0],1)
if( verbose ) {
verbfile.printf("At this step, centroid is (")
for i = 0, centrd.size-1 {
verbfile.printf(" %g",centrd.x[i])
}
verbfile.printf(") error %g\n",cent_err)
}
return cent_err
}
proc print_info_file() { local i, col
if( savepath.isopen ) savepath.close
sprint(savefit_name,"%s.info",savepath_name)
savepath.wopen(savefit_name)
savepath.printf("Sim Annealing, improved Bounded (Cardoso algorithm)\n")
savepath.printf("Initial temp was %g\n",start_temptr)
savepath.printf("Cooling rate was %g\n",cool_rate)
savepath.printf("Iterations per cooling step %d\n",iter_per_step)
savepath.printf("ftol = %g\n",ftol)
savepath.printf("gtol = %g\n",gtol)
savepath.printf("Simplex collapsed after %d function evaluations, temp %g\n",simplx_stoptrial,simplx_temptr)
col = 1
savepath.printf("Columns of file %s.fit are as follows: \n",savepath_name)
savepath.printf("[%d]\tCurrent temperature\n",col)
col += 1
savepath.printf("[%d]\tFunction evaluation count\n",col)
col += 1
savepath.printf("[%d]\tRun time\n",col)
col += 1
savepath.printf("[%d]\tMin Error\n",col)
col += 1
savepath.printf("\nThe remaining columns give the values of the following parameters:\n")
for i = 0, pf.parmlist.count-1 {
savepath.printf("[%d]\t%s\n",col,pf.parmlist.object(i).name)
col += 1
}
savepath.printf("\nThe following %d parameters were optimized:\n",pf.parmlist_use.count)
for i = 0, pf.parmlist_use.count-1 {
savepath.printf("[%d]\t%s\n",i,pf.parmlist_use.object(i).name)
}
savepath.close
}
func estimate_normgrad() { local gval, gsmall
//we've just passed a multiple of 5XN; recompute Ci & Ci-1
CiM1_avg = Ci_avg
Ci_avg = Ci / Ni
gval = (Ci_avg-CiM1_avg)/(Ci_avg*Ni)
if( gval < gtol ) { gsmall = 1 } else { gsmall = 0 }
printf("\tNormalized gradient %g < %g\n",gval,gtol)
if( verbose ) {
verbfile.printf("...Testing normalized gradient: ")
verbfile.printf("Iteration %d, multiple of %d\n",ndone,Ni)
verbfile.printf("Ci = %g, CiM1 = %g, grad ~ %g...",\
Ci_avg,CiM1_avg,(Ci_avg-CiM1_avg)/(Ci_avg*Ni))
if( !gsmall ) verbfile.printf("Gradient test failed (%g)",ftol)
verbfile.printf("\n")
}
nchk += Ni
return 0 // at the moment, we want to stop only when temperature is small enough.
//return gsmall
}
proc write_restart_file() { local i, j, use_idx
// write current simplex to restart file
printf("Writing output file %s\n",rest_name)
tl = pf.parmlist
rest_file.wopen(rest_name)
rest_file.printf("%g %d %d\n",temptr,tl.count,simplex.nrow)
for i = 0, tl.count-1 {
if( tl.object(i).doarg ) {
rest_file.printf("1 ")
} else { rest_file.printf("0 ") }
}
rest_file.printf("\n")
for j = 0, simplex.nrow-1 {
use_idx = 0
for i=0, tl.count-1 {
if( tl.object(i).doarg ) {
rest_file.printf("%-12.8g ",true_value(simplex.getrow(j).x[use_idx],use_idx))
use_idx += 1
} else {
rest_file.printf("%-12.8g ", tl.object(i).val)
}
}
rest_file.printf("\n")
}
rest_file.close
}
proc get_rstname() {
string_dialog("Enter restart filename",$s1)
}
/************************************************************
set_coded_bounds()
Often throughout the search we need to check how a point
compares to its bounds. The user may have specified that
a log transform of the data is desired, complicating this
comparison. Rather than (possibly) applying a log transform
to the bounds each time such a check is done, set up the
transformed bounds up front.
************************************************************/
proc set_coded_bounds() { local i
code_low = new Vector()
code_high = new Vector()
for i = 0, pf.parmlist_use.count-1 {
code_low.append( code_value(pf.parmlist_use.object(i).low, i))
code_high.append(code_value(pf.parmlist_use.object(i).high,i))
}
if( verbose ) {
verbfile.printf("Lower bounds were : ")
for i = 0, pf.parmlist_use.count-1 {
verbfile.printf("%g ",pf.parmlist_use.object(i).low)
}
verbfile.printf("\nTransformed to : ")
code_low.printf(verbfile)
verbfile.printf("Upper bounds were : ")
for i = 0, pf.parmlist_use.count-1 {
verbfile.printf("%g ",pf.parmlist_use.object(i).high)
}
verbfile.printf("\nTransformed to : ")
code_high.printf(verbfile)
}
}
endtemplate MulfitSA_Rctr