Skip to content
Snippets Groups Projects
Commit efbd0c80 authored by Andrew E. Torda's avatar Andrew E. Torda
Browse files

Trying to simplify main loop. Moved a few lines into their own function.

parent 9cff8eaa
Branches
No related tags found
No related merge requests found
......@@ -12,6 +12,12 @@ import (
"ackley_mc/ackwork/ackley"
)
const (
accRateIdeal = 0.05 // target of 5% acceptance rate
nRunAccAdj = 200 // every n steps, check acceptance rate
maxxDlta = 2. // max step size
)
var seedLocker sync.Mutex
// getSeed returns the random number seed, but we have to put a lock
......@@ -65,9 +71,11 @@ func setupRun(mcPrm *mcPrm, cprm *cprm) error {
}
cprm.nEvery = uint32(math.Ceil(float64(mcPrm.nStep) / float64(mcPrm.nOutput)))
// change the next line to handle the case of multiple runs
ntmp := mcPrm.fOutName + ".csv"
if cprm.fOut, err = os.Create(ntmp); err != nil {
if mcPrm.fOutName, err = setSuffix(mcPrm.fOutName, ".csv"); err != nil {
return err
}
if cprm.fOut, err = os.Create(mcPrm.fOutName); err != nil {
return err
}
if mcPrm.fPltName != "" {
......@@ -97,10 +105,10 @@ func newx(xold []float32, xT []float32, rand *rand.Rand, xDlta float32) {
}
func breaker() {}
// plotPnt adds a point to the data to be plotted. At the start, we check
// plotFval adds a point to the data to be plotted. At the start, we check
// if we have to plot at all. This is done here, to make the main loop
// below a bit shorter.
func plotPnt(n uint32, fTrial float64, plotx, ploty *[]float64, plotme bool) {
func plotFval(n uint32, fTrial float64, plotx, ploty *[]float64, plotme bool) {
if !plotme {
return
}
......@@ -108,6 +116,24 @@ func plotPnt(n uint32, fTrial float64, plotx, ploty *[]float64, plotme bool) {
*ploty = append(*ploty, float64(fTrial))
}
// nRunAdj will try to adjust the step size, given our history of accept/reject.
func nRunAdj(mcPrm *mcPrm, runAcc float64, xDlta float32) float32 {
const step = "step"
if runAcc < accRateIdeal-0.01 { // acceptance rate too low
xDlta *= 0.9
if mcPrm.verbose {
fmt.Println(step, xDlta, "decrease")
}
} else if (xDlta < maxxDlta) && (runAcc > accRateIdeal+0.01) {
xDlta *= 1.1
if mcPrm.verbose {
fmt.Println(step, xDlta, "increase")
}
}
return xDlta
}
// doRun does a Monte Carlo run. Although single precision is fine for the
// coordinates and function, we use double precision for the temperature.
func doRun(mcPrm *mcPrm) error {
......@@ -124,9 +150,6 @@ func doRun(mcPrm *mcPrm) error {
if cprm.xPlt != nil {
defer cprm.xPlt.Close()
}
const accRateIdeal = 0.05 // Arbitrary target of 5% acceptance
const nRunAccAdj = 200 // every n steps, check acceptance rate
const maxxDlta = 2. // max step size
var nAcc int // Counter of number of accepted moves
x := make([]float32, len(mcPrm.xIni)) // current position
xT := make([]float32, len(mcPrm.xIni)) // trial position
......@@ -138,7 +161,7 @@ func doRun(mcPrm *mcPrm) error {
nRunAcc := nRunAccAdj // Every nRunAccAdj, try adjusting the step size.
const runMult = 0.99
xDlta := mcPrm.xDlta // Adaptable step size
plotPnt(0, fOld, &cprm.plotx, &cprm.ploty, cprm.fplotme)
plotFval(0, fOld, &cprm.plotx, &cprm.ploty, cprm.fplotme)
for n := uint32(0); n < mcPrm.nStep; n++ {
var acc bool
nout--
......@@ -152,26 +175,15 @@ func doRun(mcPrm *mcPrm) error {
}
nRunAcc--
if nRunAcc == 0 { // Do we want to try adjusting the step size ?
const step = "step"
if runAcc < accRateIdeal-0.01 { // acceptance rate too low
xDlta *= 0.9
if mcPrm.verbose {
fmt.Println(step, xDlta, "decrease")
}
} else if (xDlta < maxxDlta) && (runAcc > accRateIdeal+0.01) {
xDlta *= 1.1
if mcPrm.verbose {
fmt.Println(step, xDlta, "increase")
}
}
nRunAcc = nRunAccAdj
xDlta = nRunAdj(mcPrm, runAcc, xDlta) // mucking about here
nRunAcc = nRunAccAdj // Reset the counter
}
newx(x, xT, cprm.rand, xDlta)
fTrial := ackley.Ackley(xT)
if fOld < fTrial {
r := cprm.rand.Float64()
delta := fTrial - fOld
r := cprm.rand.Float64() // Make the decision as to whether to
delta := fTrial - fOld // accept or reject the trial
t := math.Exp(-delta / tmprtr)
if r < t {
acc = true
......@@ -184,7 +196,7 @@ func doRun(mcPrm *mcPrm) error {
runAcc = runMult*runAcc + (1.0 - runMult)
copy(x, xT)
fOld = fTrial
plotPnt(n+1, fTrial, &cprm.plotx, &cprm.ploty, cprm.fplotme)
plotFval(n+1, fTrial, &cprm.plotx, &cprm.ploty, cprm.fplotme)
} else { // update the running estimate of acceptance
runAcc = runMult * runAcc
}
......
......@@ -82,9 +82,9 @@ func digest(prmMap map[string]string, mcPrm *mcPrm) error {
mcPrm.nRun = getu(prmMap["n_run"])
mcPrm.nOutput = getu(prmMap["n_output"])
mcPrm.xIni = getx(prmMap["x_ini"])
mcPrm.fOutName = prmMap["fOutName"]
mcPrm.fPltName = prmMap["fPltName"]
mcPrm.xPltName = prmMap["xPltName"]
mcPrm.fOutName = prmMap["foutname"]
mcPrm.fPltName = prmMap["fpltname"]
mcPrm.xPltName = prmMap["xpltname"]
seed = int64(getu(prmMap["seed"]))
if err != nil { // Only returns the first error encountered
return err
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment