diff --git a/errors.go b/errors.go index 702430b..07bf5df 100644 --- a/errors.go +++ b/errors.go @@ -65,3 +65,8 @@ func (err ErrGrad) Error() string { panic("optimize: bad ErrGrad") } } + +// List of shared panic strings +var ( + badProblem = "optimize: objective function is undefined" +) diff --git a/global.go b/global.go new file mode 100644 index 0000000..3ca8f5e --- /dev/null +++ b/global.go @@ -0,0 +1,335 @@ +// 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" + "sync" + "time" + + "github.com/gonum/matrix/mat64" +) + +// GlobalMethod is a global optimizer. Typically will require more function +// evaluations and no sense of local convergence +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 + // Global method may assume that the same task id always has the same pointer with it. + IterateGlobal(task int, loc *Location) (Operation, error) + Needser + // Done communicates to the optimization method that the optimization has + // concluded to allow for shutdown. + Done() +} + +// Global uses a global optimizer to search for the gloabl minimum of a +// function. A maximization problem can be transformed into a +// minimization problem by multiplying the function by -1. +// +// The first argument represents the problem to be minimized. Its fields are +// routines that evaluate the objective function, gradient, and other +// quantities related to the problem. The objective function, p.Func, must not +// be nil. The optimization method used may require other fields to be non-nil +// as specified by method.Needs. Global will panic if these are not met. The +// method can be determined automatically from the supplied problem which is +// described below. +// +// If p.Status is not nil, it is called before every evaluation. If the +// returned Status is not NotTerminated or the error is not nil, the +// optimization run is terminated. +// +// The third argument contains the settings for the minimization. The +// DefaultGlobalSettings function can be called for a Settings struct with the +// default values initialized. If settings == nil, the default settings are used. +// Global optimization methods typically do not make assumptions about the number +// and location of local minima. Thus, the only convergence metric used is the +// function values found at major iterations of the optimization. Bounds on the +// length of optimization are obeyed, such as the number of allowed function +// evaluations. +// +// The final argument is the optimization method to use. If method == nil, then +// an appropriate default is chosen based on the properties of the other arguments +// (dimension, gradient-free or gradient-based, etc.). +// +// If method implements Statuser, method.Status is called before every call +// to method.Iterate. If the returned Status is not NotTerminated or the +// error is non-nil, the optimization run is terminated. +// +// Global returns a Result struct and any error that occurred. See the +// documentation of Result for more information. +// +// Be aware that the default behavior of Global is to find the minimum. +// For certain functions and optimization methods, this process can take many +// function evaluations. If you would like to put limits on this, for example +// maximum runtime or maximum function evaluations, modify the Settings +// input struct. +// +// 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{} + optLoc := newLocation(dim, method) + optLoc.F = math.Inf(1) + + if settings.FunctionConverge != nil { + settings.FunctionConverge.Init(optLoc.F) + } + + 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 + if settings.Recorder != nil { + err = settings.Recorder.Record(optLoc, InitIteration, stats) + // TODO(btracey): Handle this error? Fix when merge with Local. + } + + 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{ + Location: *optLoc, + Stats: *stats, + Status: status, + }, err +} + +func minimizeGlobal(p *Problem, method GlobalMethod, settings *Settings, stats *Stats, optLoc *Location, startTime time.Time) (status Status, err error) { + dim := len(optLoc.X) + gs := &globalStatus{ + mux: &sync.RWMutex{}, + stats: stats, + status: NotTerminated, + p: p, + startTime: startTime, + optLoc: optLoc, + settings: settings, + } + + nTasks := settings.Concurrent + nTasks = method.InitGlobal(nTasks) + + // Launch optimization workers + var wg sync.WaitGroup + for task := 0; task < nTasks; task++ { + wg.Add(1) + go func(task int) { + defer wg.Done() + loc := newLocation(dim, method) + x := make([]float64, dim) + globalWorker(task, method, gs, loc, x) + }(task) + } + wg.Wait() + method.Done() + return gs.status, gs.err +} + +type globalStatus struct { + mux *sync.RWMutex + stats *Stats + status Status + p *Problem + startTime time.Time + optLoc *Location + settings *Settings + method GlobalMethod + err error +} + +func globalWorker(task int, m GlobalMethod, g *globalStatus, loc *Location, x []float64) { + for { + // Find Evaluation location + op, err := m.IterateGlobal(task, loc) + if err != nil { + // TODO(btracey): Figure out how to handle errors properly. Shut + // everything down? Pass to globalStatus so it can shut everything down? + g.mux.Lock() + g.err = err + g.status = Failure + g.mux.Unlock() + break + } + + // Evaluate location and/or update stats. + status := g.globalOperation(op, loc, x) + if status != NotTerminated { + break + } + } +} + +// globalOperation updates handles the status received by an individual worker. +// 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. + g.mux.RLock() + s := g.status + g.mux.RUnlock() + if s != NotTerminated { + return s + } + switch op { + case NoOperation: + case InitIteration: + panic("optimize: GlobalMethod return 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) + } + g.mux.Unlock() + g.mux.RLock() + status := checkConvergenceGlobal(g.optLoc, g.settings) + 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() + } + } + + // 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) + if g.status == NotTerminated { + g.status = s + } + 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 +} + +func DefaultSettingsGlobal() *Settings { + return &Settings{ + FunctionThreshold: math.Inf(-1), + FunctionConverge: &FunctionConverge{ + Absolute: 1e-10, + Iterations: 100, + }, + } +} diff --git a/guessandcheck.go b/guessandcheck.go new file mode 100644 index 0000000..40d615c --- /dev/null +++ b/guessandcheck.go @@ -0,0 +1,38 @@ +// 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 "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. +type GuessAndCheck struct { + Rander distmv.Rander + + eval []bool +} + +func (g *GuessAndCheck) Needs() struct{ Gradient, Hessian bool } { + return struct{ Gradient, Hessian bool }{false, false} +} + +func (g *GuessAndCheck) Done() { + // No cleanup needed +} + +func (g *GuessAndCheck) InitGlobal(tasks int) int { + g.eval = make([]bool, tasks) + return tasks +} + +func (g *GuessAndCheck) IterateGlobal(task int, loc *Location) (Operation, error) { + if g.eval[task] { + g.eval[task] = false + return MajorIteration, nil + } + g.eval[task] = true + g.Rander.Rand(loc.X) + return FuncEvaluation, nil +} diff --git a/guessandcheck_test.go b/guessandcheck_test.go new file mode 100644 index 0000000..e751581 --- /dev/null +++ b/guessandcheck_test.go @@ -0,0 +1,34 @@ +// 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 ( + "testing" + + "github.com/gonum/matrix/mat64" + "github.com/gonum/optimize/functions" + "github.com/gonum/stat/distmv" +) + +func TestGuessAndCheck(t *testing.T) { + dim := 3000 + problem := Problem{ + Func: functions.ExtendedRosenbrock{}.Func, + } + mu := make([]float64, dim) + sigma := mat64.NewSymDense(dim, nil) + for i := 0; i < dim; i++ { + sigma.SetSym(i, i, 1) + } + d, ok := distmv.NewNormal(mu, sigma, nil) + if !ok { + panic("bad test") + } + Global(problem, dim, nil, &GuessAndCheck{Rander: d}) + settings := DefaultSettingsGlobal() + settings.Concurrent = 5 + settings.MajorIterations = 15 + Global(problem, dim, settings, &GuessAndCheck{Rander: d}) +} diff --git a/interfaces.go b/interfaces.go index 36dc888..e28ed3a 100644 --- a/interfaces.go +++ b/interfaces.go @@ -44,6 +44,10 @@ type Method interface { // updates loc and returns the next operation. Iterate(loc *Location) (Operation, error) + Needser +} + +type Needser interface { // Needs specifies information about the objective function needed by the // optimizer beyond just the function value. The information is used // internally for initialization and must match evaluation types returned diff --git a/types.go b/types.go index f02c7c0..bf224f5 100644 --- a/types.go +++ b/types.go @@ -142,7 +142,7 @@ type Problem struct { // TODO(btracey): Think about making this an exported function when the // constraint interface is designed. -func (p Problem) satisfies(method Method) error { +func (p Problem) satisfies(method Needser) error { if method.Needs().Gradient && p.Grad == nil { return errors.New("optimize: problem does not provide needed Grad function") } @@ -228,6 +228,9 @@ type Settings struct { HessEvaluations int Recorder Recorder + + // Concurrent represents how many concurrent evaluations are possible. + Concurrent int } // DefaultSettings returns a new Settings struct containing the default settings.