diff --git a/global.go b/global.go index 3ca8f5e..f80701f 100644 --- a/global.go +++ b/global.go @@ -5,12 +5,9 @@ package optimize import ( - "fmt" "math" "sync" "time" - - "github.com/gonum/matrix/mat64" ) // GlobalMethod is a global optimizer. Typically will require more function @@ -18,7 +15,7 @@ import ( type GlobalMethod interface { // Global tells method the max number of tasks, method returns how many it wants. // This is needed to sync the Global goroutines and inside goroutines. - InitGlobal(tasks int) int + InitGlobal(dim, tasks int) int // Global method may assume that the same task id always has the same pointer with it. IterateGlobal(task int, loc *Location) (Operation, error) Needser @@ -72,40 +69,19 @@ type GlobalMethod interface { // Something about Global cannot guarantee strict bounds on function evaluations, // iterations, etc. in the precense of concurrency. func Global(p Problem, dim int, settings *Settings, method GlobalMethod) (*Result, error) { - if p.Func == nil { - panic("optimize: objective function is undefined") - } - if dim <= 0 { - panic("optimize: impossible problem dimension") - } startTime := time.Now() if method == nil { method = &GuessAndCheck{} } - if err := p.satisfies(method); err != nil { - return nil, err - } - if p.Status != nil { - _, err := p.Status() - if err != nil { - return nil, err - } - } - if settings == nil { settings = DefaultSettingsGlobal() } - - if settings.Recorder != nil { - // Initialize Recorder first. If it fails, we avoid the (possibly - // time-consuming) evaluation at the starting location. - err := settings.Recorder.Init() - if err != nil { - return nil, err - } + stats := &Stats{} + err := checkOptimization(p, dim, method, settings.Recorder) + if err != nil { + return nil, err } - stats := &Stats{} optLoc := newLocation(dim, method) optLoc.F = math.Inf(1) @@ -115,22 +91,21 @@ func Global(p Problem, dim int, settings *Settings, method GlobalMethod) (*Resul stats.Runtime = time.Since(startTime) - // Don't need to check convergence because it can't possibly have converged. - // (No function evaluations and no starting location). - var err error + // Send initial location to Recorder if settings.Recorder != nil { err = settings.Recorder.Record(optLoc, InitIteration, stats) - // TODO(btracey): Handle this error? Fix when merge with Local. + if err != nil { + return nil, err + } } + // Run optimization var status Status status, err = minimizeGlobal(&p, method, settings, stats, optLoc, startTime) // Cleanup and collect results if settings.Recorder != nil && err == nil { - // Send the optimal location to Recorder. err = settings.Recorder.Record(optLoc, PostIteration, stats) - // TODO(btracey): Handle this error? Fix when merge with Local. } stats.Runtime = time.Since(startTime) return &Result{ @@ -142,6 +117,7 @@ func Global(p Problem, dim int, settings *Settings, method GlobalMethod) (*Resul func minimizeGlobal(p *Problem, method GlobalMethod, settings *Settings, stats *Stats, optLoc *Location, startTime time.Time) (status Status, err error) { dim := len(optLoc.X) + statuser, _ := method.(Statuser) gs := &globalStatus{ mux: &sync.RWMutex{}, stats: stats, @@ -150,10 +126,11 @@ func minimizeGlobal(p *Problem, method GlobalMethod, settings *Settings, stats * startTime: startTime, optLoc: optLoc, settings: settings, + statuser: statuser, } nTasks := settings.Concurrent - nTasks = method.InitGlobal(nTasks) + nTasks = method.InitGlobal(dim, nTasks) // Launch optimization workers var wg sync.WaitGroup @@ -180,6 +157,7 @@ type globalStatus struct { optLoc *Location settings *Settings method GlobalMethod + statuser Statuser err error } @@ -209,119 +187,46 @@ func globalWorker(task int, m GlobalMethod, g *globalStatus, loc *Location, x [] // It uses a mutex to protect updates where necessary. func (g *globalStatus) globalOperation(op Operation, loc *Location, x []float64) Status { // Do a quick check to see if one of the other workers converged in the meantime. + var status Status + var err error g.mux.RLock() - s := g.status + status = g.status g.mux.RUnlock() - if s != NotTerminated { - return s + if status != NotTerminated { + return status } switch op { case NoOperation: case InitIteration: - panic("optimize: GlobalMethod return InitIteration") + panic("optimize: Method returned InitIteration") case PostIteration: panic("optimize: Method returned PostIteration") case MajorIteration: g.mux.Lock() g.stats.MajorIterations++ - if loc.F < g.optLoc.F { - copyLocation(g.optLoc, loc) - } + copyLocation(g.optLoc, loc) g.mux.Unlock() + g.mux.RLock() - status := checkConvergenceGlobal(g.optLoc, g.settings) + status = checkConvergence(g.optLoc, g.settings, false) g.mux.RUnlock() - if status != NotTerminated { - // Update g.status, preserving the first termination status. - g.mux.Lock() - if g.status == NotTerminated { - g.status = status - } - status = g.status - g.mux.Unlock() - return status - } - default: - if !op.isEvaluation() { - panic(fmt.Sprintf("optimize: invalid evaluation %v", op)) - } - copy(x, loc.X) - if op&FuncEvaluation != 0 { - loc.F = g.p.Func(x) - g.mux.Lock() - g.stats.FuncEvaluations++ - g.mux.Unlock() - } - if op&GradEvaluation != 0 { - g.p.Grad(loc.Gradient, x) - g.mux.Lock() - g.stats.GradEvaluations++ - g.mux.Unlock() - } - if op&HessEvaluation != 0 { - g.p.Hess(loc.Hessian, x) - g.mux.Lock() - g.stats.HessEvaluations++ - g.mux.Unlock() - } + default: // Any of the Evaluation operations. + status, err = evaluate(g.p, loc, op, x) + g.mux.Lock() + updateStats(g.stats, op) + g.mux.Unlock() } - // TODO(btracey): Need to fix all these things to avoid deadlock. - // When re-do, need to make sure aren't overwritting a converged status. g.mux.Lock() - g.stats.Runtime = time.Since(g.startTime) - if g.settings.Recorder != nil { - err := g.settings.Recorder.Record(loc, op, g.stats) - if err != nil { - if g.status == NotTerminated && g.err != nil { - g.status = Failure - g.err = err - } - } - } - s = checkLimits(loc, g.stats, g.settings) + status, err = iterCleanup(status, err, g.stats, g.settings, g.statuser, g.startTime, loc, op) + // Update the termination status if it hasn't already terminated. if g.status == NotTerminated { - g.status = s + g.status = status + g.err = err } - methodStatus, methodIsStatuser := g.method.(Statuser) - if methodIsStatuser { - s, err := methodStatus.Status() - if err != nil && g.status == NotTerminated { - g.status = s - g.err = err - } - } - s = g.status g.mux.Unlock() - return s -} - -func newLocation(dim int, method Needser) *Location { - // TODO(btracey): combine this with Local. - loc := &Location{ - X: make([]float64, dim), - } - loc.F = math.Inf(1) - if method.Needs().Gradient { - loc.Gradient = make([]float64, dim) - } - if method.Needs().Hessian { - loc.Hessian = mat64.NewSymDense(dim, nil) - } - return loc -} -func checkConvergenceGlobal(loc *Location, settings *Settings) Status { - if loc.F < settings.FunctionThreshold { - return FunctionThreshold - } - if settings.FunctionConverge != nil { - status := settings.FunctionConverge.FunctionConverged(loc.F) - if status != NotTerminated { - return NotTerminated - } - } - return NotTerminated + return status } func DefaultSettingsGlobal() *Settings { diff --git a/guessandcheck.go b/guessandcheck.go index 40d615c..1ec0789 100644 --- a/guessandcheck.go +++ b/guessandcheck.go @@ -4,7 +4,12 @@ package optimize -import "github.com/gonum/stat/distmv" +import ( + "math" + "sync" + + "github.com/gonum/stat/distmv" +) // GuessAndCheck is a global optimizer that evaluates the function at random // locations. Not a good optimizer, but useful for comparison and debugging. @@ -12,6 +17,10 @@ type GuessAndCheck struct { Rander distmv.Rander eval []bool + + mux *sync.Mutex + bestF float64 + bestX []float64 } func (g *GuessAndCheck) Needs() struct{ Gradient, Hessian bool } { @@ -22,14 +31,27 @@ func (g *GuessAndCheck) Done() { // No cleanup needed } -func (g *GuessAndCheck) InitGlobal(tasks int) int { +func (g *GuessAndCheck) InitGlobal(dim, tasks int) int { g.eval = make([]bool, tasks) + g.bestF = math.Inf(1) + g.bestX = resize(g.bestX, dim) + g.mux = &sync.Mutex{} return tasks } func (g *GuessAndCheck) IterateGlobal(task int, loc *Location) (Operation, error) { + // Task is true if it contains a new function evaluation. if g.eval[task] { g.eval[task] = false + g.mux.Lock() + if loc.F < g.bestF { + g.bestF = loc.F + copy(g.bestX, loc.X) + } else { + loc.F = g.bestF + copy(loc.X, g.bestX) + } + g.mux.Unlock() return MajorIteration, nil } g.eval[task] = true diff --git a/local.go b/local.go index db58cb7..706229d 100644 --- a/local.go +++ b/local.go @@ -5,12 +5,8 @@ package optimize import ( - "fmt" "math" "time" - - "github.com/gonum/floats" - "github.com/gonum/matrix/mat64" ) // Local finds a local minimum of a minimization problem using a sequential @@ -63,43 +59,22 @@ import ( // maximum runtime or maximum function evaluations, modify the Settings // input struct. func Local(p Problem, initX []float64, settings *Settings, method Method) (*Result, error) { - if p.Func == nil { - panic("optimize: objective function is undefined") - } - if len(initX) == 0 { - panic("optimize: initial X has zero length") - } - startTime := time.Now() - + dim := len(initX) if method == nil { method = getDefaultMethod(&p) } - if err := p.satisfies(method); err != nil { - return nil, err - } - - if p.Status != nil { - _, err := p.Status() - if err != nil { - return nil, err - } - } - if settings == nil { settings = DefaultSettings() } - if settings.Recorder != nil { - // Initialize Recorder first. If it fails, we avoid the (possibly - // time-consuming) evaluation at the starting location. - err := settings.Recorder.Init() - if err != nil { - return nil, err - } + stats := &Stats{} + + err := checkOptimization(p, dim, method, settings.Recorder) + if err != nil { + return nil, err } - stats := &Stats{} optLoc, err := getStartingLocation(&p, method, initX, stats, settings) if err != nil { return nil, err @@ -109,15 +84,20 @@ func Local(p Problem, initX []float64, settings *Settings, method Method) (*Resu settings.FunctionConverge.Init(optLoc.F) } - // Runtime is the only Stats field that needs to be updated here. stats.Runtime = time.Since(startTime) - // Send optLoc to Recorder before checking it for convergence. + + // Send initial location to Recorder if settings.Recorder != nil { err = settings.Recorder.Record(optLoc, InitIteration, stats) + if err != nil { + return nil, err + } } // Check if the starting location satisfies the convergence criteria. - status := checkConvergence(optLoc, settings) + status := checkConvergence(optLoc, settings, true) + + // Run optimization if status == NotTerminated && err == nil { // The starting location is not good enough, we need to perform a // minimization. The optimal location will be stored in-place in @@ -125,6 +105,7 @@ func Local(p Problem, initX []float64, settings *Settings, method Method) (*Resu status, err = minimize(&p, method, settings, stats, optLoc, startTime) } + // Cleanup and collect results if settings.Recorder != nil && err == nil { // Send the optimal location to Recorder. err = settings.Recorder.Record(optLoc, PostIteration, stats) @@ -142,7 +123,7 @@ func minimize(p *Problem, method Method, settings *Settings, stats *Stats, optLo copyLocation(loc, optLoc) x := make([]float64, len(loc.X)) - methodStatus, methodIsStatuser := method.(Statuser) + statuser, _ := method.(Statuser) var op Operation op, err = method.Init(loc) @@ -157,60 +138,24 @@ func minimize(p *Problem, method Method, settings *Settings, stats *Stats, optLo switch op { case NoOperation: - case InitIteration: panic("optimize: Method returned InitIteration") - case PostIteration: panic("optimize: Method returned PostIteration") - case MajorIteration: - stats.MajorIterations++ copyLocation(optLoc, loc) - status = checkConvergence(optLoc, settings) - + stats.MajorIterations++ + status = checkConvergence(optLoc, settings, true) default: // Any of the Evaluation operations. - if !op.isEvaluation() { - panic(fmt.Sprintf("optimize: invalid evaluation %v", op)) - } - - if p.Status != nil { - status, err = p.Status() - if err != nil || status != NotTerminated { - return - } - } - evaluate(p, loc, op, stats, x) - } - - if settings.Recorder != nil { - stats.Runtime = time.Since(startTime) - err = settings.Recorder.Record(loc, op, stats) - if err != nil { - if status == NotTerminated { - status = Failure - } - return - } - } - - if status != NotTerminated { - return + status, err = evaluate(p, loc, op, x) + updateStats(stats, op) } - stats.Runtime = time.Since(startTime) - status = checkLimits(loc, stats, settings) - if status != NotTerminated { + status, err = iterCleanup(status, err, stats, settings, statuser, startTime, loc, op) + if status != NotTerminated || err != nil { return } - if methodIsStatuser { - status, err = methodStatus.Status() - if err != nil || status != NotTerminated { - return - } - } - op, err = method.Iterate(loc) if err != nil { status = Failure @@ -220,23 +165,6 @@ func minimize(p *Problem, method Method, settings *Settings, stats *Stats, optLo panic("optimize: unreachable") } -func copyLocation(dst, src *Location) { - dst.X = resize(dst.X, len(src.X)) - copy(dst.X, src.X) - - dst.F = src.F - - dst.Gradient = resize(dst.Gradient, len(src.Gradient)) - copy(dst.Gradient, src.Gradient) - - if src.Hessian != nil { - if dst.Hessian == nil || dst.Hessian.Symmetric() != len(src.X) { - dst.Hessian = mat64.NewSymDense(len(src.X), nil) - } - dst.Hessian.CopySym(src.Hessian) - } -} - func getDefaultMethod(p *Problem) Method { if p.Grad != nil { return &BFGS{} @@ -247,16 +175,8 @@ func getDefaultMethod(p *Problem) Method { // getStartingLocation allocates and initializes the starting location for the minimization. func getStartingLocation(p *Problem, method Method, initX []float64, stats *Stats, settings *Settings) (*Location, error) { dim := len(initX) - loc := &Location{ - X: make([]float64, dim), - } + loc := newLocation(dim, method) copy(loc.X, initX) - if method.Needs().Gradient { - loc.Gradient = make([]float64, dim) - } - if method.Needs().Hessian { - loc.Hessian = mat64.NewSymDense(dim, nil) - } if settings.UseInitialData { loc.F = settings.InitialValue @@ -289,7 +209,8 @@ func getStartingLocation(p *Problem, method Method, initX []float64, stats *Stat eval |= HessEvaluation } x := make([]float64, len(loc.X)) - evaluate(p, loc, eval, stats, x) + evaluate(p, loc, eval, x) + updateStats(stats, eval) } if math.IsInf(loc.F, 1) || math.IsNaN(loc.F) { @@ -303,79 +224,3 @@ func getStartingLocation(p *Problem, method Method, initX []float64, stats *Stat return loc, nil } - -// checkConvergence returns NotTerminated if the Location does not satisfy the -// convergence criteria given by settings. Otherwise a corresponding status is -// returned. -// Unlike checkLimits, checkConvergence is called by Local only at MajorIterations. -func checkConvergence(loc *Location, settings *Settings) Status { - if loc.Gradient != nil { - norm := floats.Norm(loc.Gradient, math.Inf(1)) - if norm < settings.GradientThreshold { - return GradientThreshold - } - } - - if loc.F < settings.FunctionThreshold { - return FunctionThreshold - } - - if settings.FunctionConverge != nil { - return settings.FunctionConverge.FunctionConverged(loc.F) - } - - return NotTerminated -} - -// checkLimits returns NotTerminated status if the various limits given by -// settings has not been reached. Otherwise it returns a corresponding status. -// Unlike checkConvergence, checkLimits is called by Local at _every_ iteration. -func checkLimits(loc *Location, stats *Stats, settings *Settings) Status { - // Check the objective function value for negative infinity because it - // could break the linesearches and -inf is the best we can do anyway. - if math.IsInf(loc.F, -1) { - return FunctionNegativeInfinity - } - - if settings.MajorIterations > 0 && stats.MajorIterations >= settings.MajorIterations { - return IterationLimit - } - - if settings.FuncEvaluations > 0 && stats.FuncEvaluations >= settings.FuncEvaluations { - return FunctionEvaluationLimit - } - - if settings.GradEvaluations > 0 && stats.GradEvaluations >= settings.GradEvaluations { - return GradientEvaluationLimit - } - - if settings.HessEvaluations > 0 && stats.HessEvaluations >= settings.HessEvaluations { - return HessianEvaluationLimit - } - - // TODO(vladimir-ch): It would be nice to update Runtime here. - if settings.Runtime > 0 && stats.Runtime >= settings.Runtime { - return RuntimeLimit - } - - return NotTerminated -} - -// evaluate evaluates the routines specified by the Operation at loc.X, stores -// the answer into loc and updates stats. loc.X is copied into x before -// evaluating in order to prevent the routines from modifying it. -func evaluate(p *Problem, loc *Location, eval Operation, stats *Stats, x []float64) { - copy(x, loc.X) - if eval&FuncEvaluation != 0 { - loc.F = p.Func(x) - stats.FuncEvaluations++ - } - if eval&GradEvaluation != 0 { - p.Grad(loc.Gradient, x) - stats.GradEvaluations++ - } - if eval&HessEvaluation != 0 { - p.Hess(loc.Hessian, x) - stats.HessEvaluations++ - } -} diff --git a/minimize.go b/minimize.go new file mode 100644 index 0000000..035a9bd --- /dev/null +++ b/minimize.go @@ -0,0 +1,202 @@ +// Copyright ©2016 The gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package optimize + +import ( + "fmt" + "math" + "time" + + "github.com/gonum/floats" + "github.com/gonum/matrix/mat64" +) + +// newLocation allocates a new locatian structure of the appropriate size. It +// allocates memory based on the dimension and the values in Needs. The initial +// function value is set to math.Inf(1). +func newLocation(dim int, method Needser) *Location { + // TODO(btracey): combine this with Local. + loc := &Location{ + X: make([]float64, dim), + } + loc.F = math.Inf(1) + if method.Needs().Gradient { + loc.Gradient = make([]float64, dim) + } + if method.Needs().Hessian { + loc.Hessian = mat64.NewSymDense(dim, nil) + } + return loc +} + +func copyLocation(dst, src *Location) { + dst.X = resize(dst.X, len(src.X)) + copy(dst.X, src.X) + + dst.F = src.F + + dst.Gradient = resize(dst.Gradient, len(src.Gradient)) + copy(dst.Gradient, src.Gradient) + + if src.Hessian != nil { + if dst.Hessian == nil || dst.Hessian.Symmetric() != len(src.X) { + dst.Hessian = mat64.NewSymDense(len(src.X), nil) + } + dst.Hessian.CopySym(src.Hessian) + } +} + +func checkOptimization(p Problem, dim int, method Needser, recorder Recorder) error { + if p.Func == nil { + panic("optimize: objective function is undefined") + } + if dim <= 0 { + panic("optimize: impossible problem dimension") + } + if err := p.satisfies(method); err != nil { + return err + } + if p.Status != nil { + _, err := p.Status() + if err != nil { + return err + } + } + if recorder != nil { + err := recorder.Init() + if err != nil { + return err + } + } + return nil +} + +// evaluate evaluates the routines specified by the Operation at loc.X, and stores +// the answer into loc. loc.X is copied into x before +// evaluating in order to prevent the routines from modifying it. +func evaluate(p *Problem, loc *Location, op Operation, x []float64) (Status, error) { + if !op.isEvaluation() { + panic(fmt.Sprintf("optimize: invalid evaluation %v", op)) + } + if p.Status != nil { + status, err := p.Status() + if err != nil || status != NotTerminated { + return status, err + } + } + copy(x, loc.X) + if op&FuncEvaluation != 0 { + loc.F = p.Func(x) + } + if op&GradEvaluation != 0 { + p.Grad(loc.Gradient, x) + } + if op&HessEvaluation != 0 { + p.Hess(loc.Hessian, x) + } + return NotTerminated, nil +} + +// checkConvergence returns NotTerminated if the Location does not satisfy the +// convergence criteria given by settings. Otherwise a corresponding status is +// returned. +// Unlike checkLimits, checkConvergence is called only at MajorIterations. +// +// If local is true, gradient convergence is also checked. +func checkConvergence(loc *Location, settings *Settings, local bool) Status { + if local && loc.Gradient != nil { + norm := floats.Norm(loc.Gradient, math.Inf(1)) + if norm < settings.GradientThreshold { + return GradientThreshold + } + } + if loc.F < settings.FunctionThreshold { + return FunctionThreshold + } + if settings.FunctionConverge != nil { + return settings.FunctionConverge.FunctionConverged(loc.F) + } + return NotTerminated +} + +// updateStats updates the statistics based on the operation. +func updateStats(stats *Stats, op Operation) { + if op&FuncEvaluation != 0 { + stats.FuncEvaluations++ + } + if op&GradEvaluation != 0 { + stats.GradEvaluations++ + } + if op&HessEvaluation != 0 { + stats.HessEvaluations++ + } +} + +// checkLimits returns NotTerminated status if the various limits given by +// settings have not been reached. Otherwise it returns a corresponding status. +// Unlike checkConvergence, checkLimits is called by Local and Global at _every_ +// iteration. +func checkLimits(loc *Location, stats *Stats, settings *Settings) Status { + // Check the objective function value for negative infinity because it + // could break the linesearches and -inf is the best we can do anyway. + if math.IsInf(loc.F, -1) { + return FunctionNegativeInfinity + } + + if settings.MajorIterations > 0 && stats.MajorIterations >= settings.MajorIterations { + return IterationLimit + } + + if settings.FuncEvaluations > 0 && stats.FuncEvaluations >= settings.FuncEvaluations { + return FunctionEvaluationLimit + } + + if settings.GradEvaluations > 0 && stats.GradEvaluations >= settings.GradEvaluations { + return GradientEvaluationLimit + } + + if settings.HessEvaluations > 0 && stats.HessEvaluations >= settings.HessEvaluations { + return HessianEvaluationLimit + } + + // TODO(vladimir-ch): It would be nice to update Runtime here. + if settings.Runtime > 0 && stats.Runtime >= settings.Runtime { + return RuntimeLimit + } + + return NotTerminated +} + +// TODO(btracey): better name +func iterCleanup(status Status, err error, stats *Stats, settings *Settings, statuser Statuser, startTime time.Time, loc *Location, op Operation) (Status, error) { + if status != NotTerminated || err != nil { + return status, err + } + + if settings.Recorder != nil { + stats.Runtime = time.Since(startTime) + err = settings.Recorder.Record(loc, op, stats) + if err != nil { + if status == NotTerminated { + status = Failure + } + return status, err + } + } + + stats.Runtime = time.Since(startTime) + status = checkLimits(loc, stats, settings) + if status != NotTerminated { + return status, nil + } + + if statuser != nil { + status, err = statuser.Status() + if err != nil || status != NotTerminated { + return status, err + } + } + return status, nil +}