diff --git a/errors.go b/errors.go index 702430b..a77f3e4 100644 --- a/errors.go +++ b/errors.go @@ -28,6 +28,8 @@ var ( // due to floating-point arithmetic. ErrNoProgress = errors.New("linesearch: no change in location after Linesearcher step") + ErrFunctionConvergeNil = errors.New("global: function convergence must not be nil") + // ErrLinesearcherBound signifies that a Linesearcher reached a step that // lies out of allowed bounds. ErrLinesearcherBound = errors.New("linesearch: step out of bounds") @@ -65,3 +67,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..ff8d60c --- /dev/null +++ b/global.go @@ -0,0 +1,248 @@ +package optimize + +import ( + "math" + "time" + + "github.com/gonum/matrix/mat64" +) + +// Global uses a global seeks to find the gloabl minimum of a minimization problem. +// A maximization problem can be transformed into a minimization problem by +// multiplying the function by -1. +// +// Global is different than local because it tries to find the global minimum, +// which can be anywhere. The only reasonable convergence metric is the updates +// on major iteration. GradientThreshold is not checked during convergence, and +// initial values are not used. +// +// TODO(btracey): Split into LocalSettings and GlobalSettings with a common Settings +// struct that is embedded? In initializeOptimization +// TODO(btracey): LocalMethod / GlobalMethod or Method and GlobalMethod? +func Global(p Problem, dim int, settings *Settings, method GlobalMethod) (*Result, error) { + if method == nil { + // use getDefaultMethod with wrapper + panic("need method") + } + if settings == nil { + settings = DefaultSettings() + } + startTime, err := initializeOptimization(&p, settings, method) + if err != nil { + return nil, err + } + if settings.FunctionConverge == nil { + return nil, ErrFunctionConvergeNil + } + stats := &Stats{} + optLoc := &Location{ + F: math.Inf(1), + X: make([]float64, dim), + } + if method.Needs().Gradient { + optLoc.Gradient = make([]float64, dim) + } + if method.Needs().Hessian { + optLoc.Hessian = mat64.NewSymDense(dim, nil) + } + + status, err := minimizeGlobal(p, method, settings, stats, optLoc, startTime) + + if settings.Recorder != nil && err == nil { + // Send the optimal location to Recorder. + err = settings.Recorder.Record(optLoc, PostIteration, stats) + } + 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) { + concurrent := settings.Concurrent + if concurrent == 0 { + concurrent = 1 + } + + jobs := make(chan LocationSender) + quit := make(chan struct{}) + updater := make(chan LocationSender, settings.Concurrent) + closer := make(chan struct{}, settings.Concurrent) + + go func() { + method.Schedule(len(optLoc.X), concurrent, jobs, quit) + }() + + // High level: Method continually sends locations and operations on the + // jobs channel. These are received by the workers, who evaluate the function, + // and send the result to the updater. The updater updates optLoc if necessary + // and tests for convergence. + // On convergence, globalUpdater closes the quit channel. It the the methods + // job to stop sending values and close the jobs channel. Worker sees the jobs + // channel has closed, and signaled it has completed. When all of the workers + // have completed, the updater channel is closed so the updater knows to return. + + // Local has the property that the Method keeps feeding locations until done. + // This is an important property. + + // Launch worker goroutines + for i := 0; i < concurrent; i++ { + go globalWorker(p, jobs, updater, closer) + } + + // This closes the updater channel once all of the workers have quit. This + // then causes updater to quit. + go func() { + for i := 0; i < concurrent; i++ { + <-closer + } + close(updater) + }() + + // Keep updating until quit + status = globalUpdater(stats, settings, updater, quit, optLoc, startTime) + return status, nil +} + +// Worker receives a location and an operation. Worker sends the result of that +// operation to globalUpdater (as stats updates and checking for convergence +// must happen in serial). Once the needed information from Location is used, +// it sends the location back to the GlobalMethod so it can use the result. + +// either we have to copy location or we have to deal with the fact this could +// be a bottle neck. +func globalUpdater(stats *Stats, settings *Settings, updater <-chan LocationSender, quit chan struct{}, optLoc *Location, startTime time.Time) Status { + first := true + var status Status + var closed bool + for update := range updater { + // Get all of the needed information on the global update so can let + // the worker continue as fast as possible + eval := update.Operation + var f float64 + if eval == MajorIteration { + f = update.Location.F + if f < optLoc.F { + copyLocation(optLoc, update.Location) + } + } + update.Evaluated <- struct{}{} + + // TODO(btracey): Check that it's a valid operation + if eval == MajorIteration { + if first { + settings.FunctionConverge.Init(f) + first = false + } else { + status = settings.FunctionConverge.FunctionConverged(f) + } + } else { + if eval&FuncEvaluation != 0 { + stats.FuncEvaluations++ + } + if eval&GradEvaluation != 0 { + stats.GradEvaluations++ + } + if eval&HessEvaluation != 0 { + stats.HessEvaluations++ + } + } + + // TODO(btracey): Add statuser check for both Problem and Optimize. + settings.Runtime = time.Since(startTime) + if status == NotTerminated { + if optLoc.F < settings.FunctionThreshold { + status = FunctionThreshold + } else { + status = checkLimits(optLoc, stats, settings) + } + } + if !closed && status != NotTerminated { + closed = true + close(quit) + } + } + return status +} + +// globalWorker receives a location and operation to evaluate from the global +// queue, evaluates the function, and sends the result to the global updater. +// OptLoc is ONLY to keep the sizing code in one place. Location is temporary +// storage. +// Method's job to close queue once quit happens +func globalWorker(p Problem, jobs <-chan LocationSender, updater chan<- LocationSender, closer chan<- struct{}) { + var x []float64 + for job := range jobs { + eval := job.Operation + if x == nil { + x = make([]float64, len(job.Location.X)) + } + copy(x, job.Location.X) + if eval&FuncEvaluation != 0 { + job.Location.F = p.Func(x) + } + if eval&GradEvaluation != 0 { + p.Grad(job.Location.Gradient, x) + } + if eval&HessEvaluation != 0 { + p.Hess(job.Location.Hessian, x) + } + updater <- job + } + // Communicate that has finished + closer <- struct{}{} +} + +func initializeOptimization(p *Problem, settings *Settings, method GlobalMethod) (start time.Time, err error) { + if p.Func == nil { + panic(badProblem) + } + startTime := time.Now() + if err := p.satisfies(method); err != nil { + return startTime, err + } + if p.Status != nil { + _, err := p.Status() + if err != nil { + return startTime, err + } + } + 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 startTime, err + } + } + return startTime, nil +} + +func DefaultSettingsGlobal() *Settings { + return &Settings{ + FunctionThreshold: math.Inf(-1), + FunctionConverge: &FunctionConverge{ + Absolute: 1e-10, + Iterations: 20, + }, + } +} + +// LocationSender sends a location to be evaluated and its iteration. +// The location with values populated will be returned on the receive channel. +type LocationSender struct { + Location *Location + Operation Operation + Evaluated chan<- struct{} +} + +// Init initializes the GlobalMethod. It continues to send new locations on send +// until a message from quit is sent. +// Method must close jobs when a receive from quit happens +type GlobalMethod interface { + Schedule(dim, concurrent int, jobs chan<- LocationSender, quit <-chan struct{}) + + Needser +} diff --git a/guessandcheck.go b/guessandcheck.go new file mode 100644 index 0000000..ab030a4 --- /dev/null +++ b/guessandcheck.go @@ -0,0 +1,59 @@ +package optimize + +import "sync" + +// TODO(btracey): Add Rander interface to distmv. +type Rander interface { + Rand([]float64) []float64 +} + +// 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 Rander +} + +func (g GuessAndCheck) Needs() struct{ Gradient, Hessian bool } { + return struct{ Gradient, Hessian bool }{false, false} +} + +func (g GuessAndCheck) Schedule(dim, concurrent int, jobs chan<- LocationSender, quit <-chan struct{}) { + var wg sync.WaitGroup + wg.Add(concurrent) + for i := 0; i < concurrent; i++ { + go func() { + defer wg.Done() + done := make(chan struct{}) + location := &Location{ + X: make([]float64, dim), + } + for { + g.Rander.Rand(location.X) + + // See if it's time to quit + select { + case <-quit: + return + default: + } + + // Evaluate the function and then send a major iteration update. + jobs <- LocationSender{ + Operation: FuncEvaluation, + Location: location, + Evaluated: done, + } + <-done + + jobs <- LocationSender{ + Operation: MajorIteration, + Location: location, + Evaluated: done, + } + <-done + } + }() + } + wg.Wait() + close(jobs) +} diff --git a/guessandcheck_test.go b/guessandcheck_test.go new file mode 100644 index 0000000..1d2fcc3 --- /dev/null +++ b/guessandcheck_test.go @@ -0,0 +1,29 @@ +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 = 10 + 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/local.go b/local.go index db58cb7..c32f07c 100644 --- a/local.go +++ b/local.go @@ -64,7 +64,7 @@ import ( // input struct. func Local(p Problem, initX []float64, settings *Settings, method Method) (*Result, error) { if p.Func == nil { - panic("optimize: objective function is undefined") + panic(badProblem) } if len(initX) == 0 { panic("optimize: initial X has zero length") diff --git a/randomlocal.go b/randomlocal.go new file mode 100644 index 0000000..55a1617 --- /dev/null +++ b/randomlocal.go @@ -0,0 +1,59 @@ +package optimize + +/* +// RestartLocal performs global optimization through random-restart local +// optimization. +// TODO(btracey): Some way of making Method non-nil? Otherwise can't do +// &RestartLocal{} +type RestartLocal struct { + Rander Rander + MethodFactory func() Method // Method for the local optimization + SettingFactory func() *Settings // Settings for the local optimization +} + +func (r *RestartLocal) Needs() { + return r.Method.Needs() +} + +func (r *RestartLocal) Schedule(dim, concurrent int, jobs chan<- LocationSender, quit chan<- struct{}) { + // Launch concurrent local optimizers + for i := 0; i < concurrent; i++ { + go func() { + initialize := true + method := r.MethodFactory() + + var location Location + location.X = make([]float64, dim) + if method.Needs().Gradient { + location.Gradient = make([]float64, dim) + } + if method.Needs().Hessian { + location.Hessian = mat64.NewSymDense(dim, nil) + } + evaluated := make(chan struct{}) + + locationSend := LocationSender{ + Location: location, + Evaluated: evaluated, + } + + for { + if initialize { + // Starting interation of an optimization + // Use rander to find an initialX + r.Rander.Rand(initX) + if method.Needs().Gradient { + locationSend.Operation = GradEvaluation + jobs <- locationSend + + } + + stats := &Stats{} + settings := r.SettingFactory() + method.Init(loc) + } + } + }() + } +} +*/ diff --git a/types.go b/types.go index f02c7c0..f28980b 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 + + // For global, how many concurrent. + Concurrent int } // DefaultSettings returns a new Settings struct containing the default settings.