dorun.go

// Aug 2021
package ackwork

import (
	"fmt"
	"io"
	"math"
	"math/rand"
	"os"
	"sync"

	"ackley_mc/ackwork/ackley"
)

var seedLocker sync.Mutex

// getSeed returns the random number seed, but we have to put a lock
// around it so multiple runs get consistent successive seeds
func getSeed() int64 {
	seedLocker.Lock()
	r := seed
	seed++
	seedLocker.Unlock()
	return r
}

type cprm struct { // parameters calculated from input
	rand     *rand.Rand
	coolMult float64 // multiplier for temperature
	coolme   bool    // Are we cooling or just doing constant temperature ?
	plotme   bool    // Are we making output for plotting
	nEvery   uint32  // Output every nEvery steps
	fOut     io.WriteCloser
	plotOut  io.WriteCloser
	plotx    []float64 // Why float 64 ? Because the plotting libraries
	ploty    []float64 // generally want this
}

// setupRun does things like get the cooling rate, seed the random numbers.
func setupRun(mcPrm *mcPrm, cprm *cprm) error {
	var err error
	if mcPrm.dummy == true {
		return nil
	}
	cprm.rand = rand.New(rand.NewSource(getSeed()))

	if mcPrm.iniTmp != mcPrm.fnlTmp { // cooling or constant temperature ?
		var coolrate float64
		cprm.coolme = true
		if mcPrm.fnlTmp > mcPrm.iniTmp {
			return fmt.Errorf("Final temp %g higher than initial temp %g", mcPrm.fnlTmp, mcPrm.iniTmp)
		}
		nStep := float64(mcPrm.nStep)
		fnlTmp, iniTmp := float64(mcPrm.fnlTmp), float64(mcPrm.iniTmp)
		if fnlTmp == 0 {
			fnlTmp = math.SmallestNonzeroFloat32
		}
		coolrate = -1 / nStep
		coolrate *= math.Log(fnlTmp / iniTmp)
		cprm.coolMult = math.Exp(-coolrate)
	}

	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 {
		return err
	}
	if mcPrm.pOutName != "" {
		mcPrm.pOutName = makepngName(mcPrm.pOutName)
		if cprm.plotOut, err = os.Create(mcPrm.pOutName); err != nil {
			return err
		}
		cprm.plotme = true
		n_alloc := mcPrm.nStep / 5 // About a fifth of the number of steps
		cprm.plotx = make([]float64, 0, n_alloc)
		cprm.ploty = make([]float64, 0, n_alloc)
	}
	return nil
}

// newx gets a candidate x slice. We have n dimensions, so a candidate
// move is a slice of length n.
func newx(xold []float32, xT []float32, rand *rand.Rand, xDlta float32) {
	for i, x := range xold {
		t := 2.*rand.Float32() - 1
		t *= xDlta
		xT[i] = x + t
	}
}
func breaker() {}

// plotPnt 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) {
	if !plotme {
		return
	}
	*plotx = append(*plotx, float64(n))
	*ploty = append(*ploty, float64(fTrial))
}

// 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 {
	var cprm cprm
	if err := setupRun(mcPrm, &cprm); err != nil {
		return err
	}
	if cprm.fOut != nil {
		defer cprm.fOut.Close()
	}
	if cprm.plotOut != nil {
		defer cprm.plotOut.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
	nout := uint32(1)                      // counter for printing output
	runAcc := accRateIdeal                 // Running acceptance rate, exponentially weighted
	copy(x, mcPrm.xIni)
	fOld := ackley.Ackley(x) // Initial function value
	tmprtr := float64(mcPrm.iniTmp)
	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.plotme)
	for n := uint32(0); n < mcPrm.nStep; n++ {
		var acc bool
		nout--
		if nout == 0 { // Do we want to print out results on this step ?
			fmt.Fprintf(cprm.fOut, "%d,%.4g,%.5g", n, tmprtr, fOld)
			for _, xx := range x {
				fmt.Fprintf(cprm.fOut, ",%.2g", xx)
			}
			fmt.Fprintln(cprm.fOut)
			nout = cprm.nEvery
		}
		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
		}
		newx(x, xT, cprm.rand, xDlta)
		fTrial := ackley.Ackley(xT)

		if fOld < fTrial {
			r := cprm.rand.Float64()
			delta := fTrial - fOld
			t := math.Exp(-delta / tmprtr)
			if r < t {
				acc = true
			}
		} else {
			acc = true
		}
		if acc == true {
			nAcc++
			runAcc = runMult*runAcc + (1.0 - runMult)
			copy(x, xT)
			fOld = fTrial
			plotPnt(n+1, fTrial, &cprm.plotx, &cprm.ploty, cprm.plotme)
		} else { // update the running estimate of acceptance
			runAcc = runMult * runAcc
		}
		if cprm.coolme {
			tmprtr *= cprm.coolMult
		}
	}
	err := plotWrite(cprm.plotx, cprm.ploty, cprm.plotOut, cprm.plotme)
	fmt.Println("n accepted:", nAcc, "of", mcPrm.nStep+1)
	return err
}