optimize

go get github.com/go-numerics/optimize@v0.1.0

// scipy.optimize.minimize(rosen, x0, jac=rosen_der, method='BFGS')
r, err := optimize.Minimize(optimize.Problem{Func: rosen, Grad: rosenDer},
	[]float64{-1.2, 1}, optimize.BFGS, nil)
// r.X = [1 1], r.Nit = 32 (SciPy: 32), "Optimization terminated successfully."

// scipy.optimize.minimize(f, x0, method='L-BFGS-B', bounds=[(-2, 0.5), (None, None)])
r, _ = optimize.Minimize(optimize.Problem{Func: f}, []float64{0, 0}, optimize.LBFGSB,
	&optimize.Options{Bounds: []optimize.Bound{{Lo: -2, Hi: 0.5}, optimize.Unbounded()}})

// popt, pcov = scipy.optimize.curve_fit(lambda x, a, b, c: a*np.exp(-b*x) + c, x, y)
fit, err := optimize.CurveFit(func(x float64, p []float64) float64 {
	return p[0]*math.Exp(-p[1]*x) + p[2]
}, xdata, ydata, []float64{1, 1, 1}, nil)
// fit.Popt, fit.PcovArray(), fit.Perr()  ->  a = 2.500 ± 0.004, b = 1.300 ± 0.005, c = 0.500 ± 0.002

// scipy.optimize.brentq(lambda x: x**3 - 2*x - 5, 2, 3)
root, _ := optimize.Brentq(func(x float64) float64 { return x*x*x - 2*x - 5 }, 2, 3, nil)
// root.Root = 2.094551481542…, 7 iterations (as SciPy)

// scipy.optimize.root(fun, [0, 0])        (fun writes f(x) into its second argument)
sol, _ := optimize.Root(optimize.SystemProblem{Func: fun}, []float64{0, 0}, nil)

// scipy.optimize.least_squares(fun, x0, bounds=([0, 0], [2, 0.5]), loss='soft_l1', f_scale=0.1)
ls, _ := optimize.LeastSquares(optimize.LSProblem{Func: resid, M: m}, []float64{1, 0.1}, optimize.TRF,
	&optimize.LSOptions{Bounds: []optimize.Bound{{0, 2}, {0, 0.5}}, Loss: optimize.LossSoftL1, FScale: 0.1})

// scipy.optimize.linprog(c, A_ub=A, b_ub=b, bounds=[(None, None), (-3, None)])
lp, _ := optimize.Linprog([]float64{-1, 4}, optimize.LPConstraints{AUb: A, BUb: b,
	Bounds: []optimize.Bound{optimize.Unbounded(), {Lo: -3, Hi: math.Inf(1)}}}, nil)
// lp.X = [10 -3], lp.Fun = -22, lp.IneqLin.Marginals (HiGHS sign convention)

package main

import (
	"fmt"
	"math"

	"github.com/go-ndarray/ndarray"
	"github.com/go-numerics/optimize"
	"github.com/go-plot/plot"
)

func init() { register("curvefit", curvefit) }

// The README's CurveFit example: 50 points of 2.5 exp(-1.3x) + 0.5 with a
// deterministic ripple, fitted from p0 = (1, 1, 1). The legend prints the
// parameters and standard errors CurveFit returned.
func curvefit() *plot.Figure {
	x, _ := ndarray.Linspace(0, 4, 50)
	y := x.Map(func(v float64) float64 { return 2.5*math.Exp(-1.3*v) + 0.5 + 0.01*math.Sin(37*v) })
	model := func(x float64, p []float64) float64 { return p[0]*math.Exp(-p[1]*x) + p[2] }
	r, err := optimize.CurveFit(model, x, y, []float64{1, 1, 1}, nil)
	if err != nil {
		panic(err)
	}
	perr := r.Perr()
	xs := linspace(0, 4, 200)
	fit := make([]float64, len(xs))
	for i, v := range xs {
		fit[i] = model(v, r.Popt)
	}
	fig, axs := plot.Subplots(1, 1, plot.FigSize(6.4, 4))
	ax := axs[0][0]
	ax.Plot(x, y, plot.Marker("o"), plot.LineStyle("None"), plot.MarkerSize(4), plot.Label("data"))
	ax.Plot(xs, fit, plot.Label(fmt.Sprintf(`$a e^{-bx} + c$: a = %.3f ± %.3f, b = %.3f ± %.3f, c = %.3f ± %.3f`,
		r.Popt[0], perr[0], r.Popt[1], perr[1], r.Popt[2], perr[2])))
	ax.SetTitle(fmt.Sprintf("CurveFit (Levenberg–Marquardt): %d evaluations", r.Nfev))
	ax.SetXLabel(`$x$`)
	ax.Legend(plot.FontSize(8))
	return fig
}

package main

import (
	"math"

	"github.com/go-numerics/optimize"
	"github.com/go-plot/plot"
)

func init() { register("robust", robust) }

// The README's LeastSquares example, fitted twice: with the ordinary squared
// loss, which the outlier at t = 5 drags upwards, and with soft_l1, which
// keeps it from dominating. Same bounds, same start.
func robust() *plot.Figure {
	t := []float64{0, 1, 2, 3, 4, 5}
	y := []float64{1.0, 0.62, 0.36, 0.23, 0.13, 3.0} // the last point is an outlier
	p := optimize.LSProblem{M: len(t), Func: func(x, f []float64) {
		for i := range t {
			f[i] = x[0]*math.Exp(-x[1]*t[i]) - y[i]
		}
	}}
	bounds := []optimize.Bound{{Lo: 0, Hi: 2}, {Lo: 0, Hi: 0.5}}
	linear, err := optimize.LeastSquares(p, []float64{1, 0.1}, optimize.TRF,
		&optimize.LSOptions{Bounds: bounds})
	if err != nil {
		panic(err)
	}
	soft, err := optimize.LeastSquares(p, []float64{1, 0.1}, optimize.TRF,
		&optimize.LSOptions{Bounds: bounds, Loss: optimize.LossSoftL1, FScale: 0.1})
	if err != nil {
		panic(err)
	}
	xs := linspace(0, 5, 200)
	curve := func(x []float64) []float64 {
		v := make([]float64, len(xs))
		for i, s := range xs {
			v[i] = x[0] * math.Exp(-x[1]*s)
		}
		return v
	}
	fig, axs := plot.Subplots(1, 1, plot.FigSize(6.4, 4))
	ax := axs[0][0]
	ax.Plot(t, y, plot.Marker("o"), plot.LineStyle("None"), plot.Label("data"))
	ax.Plot(xs, curve(linear.X), plot.LineStyle("--"), plot.Label("squared loss"))
	ax.Plot(xs, curve(soft.X), plot.Label("soft_l1 loss, f_scale = 0.1"))
	ax.SetTitle("LeastSquares (TRF, bounds): one outlier, two losses")
	ax.SetXLabel(`$t$`)
	ax.Legend()
	return fig
}