diff --git a/docs/source/api/api-python/optim/index.rst b/docs/source/api/api-python/optim/index.rst index 7e847d53..b13d7d21 100644 --- a/docs/source/api/api-python/optim/index.rst +++ b/docs/source/api/api-python/optim/index.rst @@ -27,22 +27,12 @@ Module Contents Signature of the objective function :math:`f(, y)`. - !!!!!!!!! The function accepts four arguments: - !!!!!!!!! - `t`: current time, - !!!!!!!!! - `y`: state vector at time :math:`t`, - !!!!!!!!! - `ydot`: output array to which the result of function evalutation is stored, - !!!!!!!!! - `user_data`: additional context (user-defined data) that - !!!!!!!!! must be passed to the function (e.g., parameters of the system). - -.. py:class:: OptimResult - - .. py:attribute:: status - :type: int - - - .. py:attribute:: x - :type: numpy.ndarray - + !!!!!!!!! The function accepts four arguments: + !!!!!!!!! - `t`: current time, + !!!!!!!!! - `y`: state vector at time :math:`t`, + !!!!!!!!! - `ydot`: output array to which the result + !!!!!!!!! - `user_data`: additional context (user-defined data) that + !!!!!!!!! must be passed to the function (e.g., parameters of the system). .. py:class:: Optim(impl: str) @@ -56,54 +46,40 @@ Module Contents .. rubric:: Examples - Let's solve the following initial value problem: + Let's solve the following convex optimization problem: .. math:: - y'(t) = -y(t), \quad y(0) = 1. + minimize \sum_{i = 1}^N x_i^2 + + where the solution is :math:`[0, ..., 0]` as the problem is convex. First, import the necessary modules: >>> import numpy as np - >>> from oif.interfaces.ivp import IVP + >>> from oif.interfaces.optim import Optim + + Define the objective function: - Define the right-hand side function: + >>> def objective_fn(x): + ... return np.sum(x**2) - >>> def rhs(t, y, ydot, user_data): - ... ydot[0] = -y[0] - ... return 0 # No errors, optional + Create an instance of the optim solver using the implementation "scipy_optimize", + which is an adapter to the `scipy.optimize` Python package: - Now define the initial condition: + >>> s = Optim("scipy_optimize") - >>> y0, t0 = np.array([1.0]), 0.0 + We set the initial value, the right-hand side function, and the tolerance: - Create an instance of the IVP solver using the implementation "jl_diffeq", - which is an adapter to the `OrdinaryDiffeq.jl` Julia package: + >>> s.set_initial_guess([2.718, 3.142]) + >>> s.set_objective_fn(objective_fn) - >>> s = IVP("jl_diffeq") + Now we solve the minimization problem and print the return status and message: - We set the initial value, the right-hand side function, and the tolerance: + >>> status, message = s.minimize() + + We can print the resultant minimizer by retrieving it from the solver: + + >>> print(f"Minimizer is {s.x}") - >>> s.set_initial_value(y0, t0) - >>> s.set_rhs_fn(rhs) - >>> s.set_tolerances(1e-6, 1e-12) - - Now we integrate to time `t = 1.0` in a loop, outputting the current value - of `y` with time step `0.1`: - - >>> t = t0 - >>> times = np.linspace(t0, t0 + 1.0, num=11) - >>> for t in times[1:]: - ... s.integrate(t) - ... print(f"{t:.1f} {s.y[0]:.6f}") - 0.1 0.904837 - 0.2 0.818731 - 0.3 0.740818 - 0.4 0.670320 - 0.5 0.606531 - 0.6 0.548812 - 0.7 0.496585 - 0.8 0.449329 - 0.9 0.406570 - 1.0 0.367879 .. py:attribute:: x0 @@ -112,13 +88,16 @@ Module Contents Current value of the state vector. - .. py:attribute:: status - :value: -1 + .. py:attribute:: x + :type: numpy.ndarray + + .. py:attribute:: user_data + :type: object - .. py:attribute:: x - :type: numpy.ndarray + .. py:attribute:: oif_user_data + :type: openinterfaces.core.OIFUserData .. py:method:: set_initial_guess(x0: numpy.ndarray) @@ -130,6 +109,21 @@ Module Contents .. py:method:: set_objective_fn(objective_fn: ObjectiveFn) + .. py:method:: set_grad_fn(grad_fn) + + + .. py:method:: set_user_data(user_data: object) + + Specify additional data that will be used for right-hand side function. + + + + .. py:method:: set_method(method_name: str, method_params: dict = {}) + + Set integrator, if the name is recognizable. + + + .. py:method:: minimize() Integrate to time `t` and write solution to `y`. diff --git a/docs/source/getting-started/example-optim-rosenbrock-julia.md b/docs/source/getting-started/example-optim-rosenbrock-julia.md new file mode 100644 index 00000000..07cc93ae --- /dev/null +++ b/docs/source/getting-started/example-optim-rosenbrock-julia.md @@ -0,0 +1,155 @@ +# Example: Minimize the Rosenbrock function using Julia + +This tutorial shows how to use _Open Interfaces_ to +minimize the +[Rosenbrock function][rosenbrock-wiki] using Julia `optim` gateway. + +The problem is mathematically defined as +finding the minimizing vector $x \in \mathbb{R}^n$ +for the objective function + +$$ + f(x) = \sum_{i=1}^{n-1} \left(a (x_{i+1} - x_i^2)^2 + + (1 - x_i)^2\right), +$$ + +where $a = 10$. The global minimizer is +$x = [1, 1, \ldots, 1]$, where $f(x) = 0$. + +## Loading the interface + +The example begins by importing the optimization interface: + +```julia +using OpenInterfaces.Interfaces.Optim +``` + +The first call obtains an optimizer instance for a specific backend implementation: + +```julia +s = Optim.Self(impl) +``` + +The `impl` string selects the optimization backend. +The supported values for the `optim` interface are: + +- `optim_jl` +- `scipy_optimize` + +## Defining the objective and gradient + +Let's define the objective and gradient functions: + +```julia +function rosenbrock_objective_fn(x, a) + return sum(a * (x[2:end] - x[1:(end-1)] .^ 2.0) .^ 2.0 + (1 .- x[1:(end-1)]) .^ 2.0) +end + +function rosenbrock_grad_fn(x, grad_f, a) + xi = @view x[1:(end-1)] + xip1 = @view x[2:end] + + grad_f[1:(end-1)] .= -4.0 .* a .* xi .* (xip1 .- xi .^ 2.0) .- 2.0 .* (1.0 .- xi) + grad_f[end] = 0.0 + grad_f[2:end] .+= 2.0 .* a .* (xip1 .- xi .^ 2.0) + return 0 +end +``` + + +## Configuring the optimizer + +First of all we instantiate an instance of the `optim` interface. +As the Julia components use the object-oriented approach, +we instantiate the provided structure `Self`: + +```julia +s = Optim.Self("scipy_optimize") +``` +where we use `scipy_optimize` (Python implementation). + +Next we need to set an initial guess for the nonlinear optimization process: + +```julia +x0 = [3.14, 2.72, 6.18, 9.81, 8.31] + +Optim.set_initial_guess(s, x0) +``` + +Our objective function takes an additional parameter $a$, +so we need to pass it as _user data_: +```julia +user_data = 10 +Optim.set_user_data(s, user_data) +``` + +Now we can select the optimization method, +that we want to use and their respective parameters. +These parameters are always method-specific, and should be +looked for in the documentation for the given method. + +For example, we can select the `NelderMead` algorithm, +which in SciPy allows to set the termination criterion +to be the norm of the objective function: + +```julia +Optim.set_method(s, "nelder-mead", Dict("fatol" => 1e-11)) +``` + +or we can select the `BFGS` method, +with the termination criterion defined as the inf-norm +of the gradient function: +```julia +Optim.set_method(s, "BFGS", Dict("gtol" => 1e-8)) +``` + +The objective and gradient callbacks are then passed +to the implementations (note that the `Nelder-Mead` does not require +`set_grad_fn`, so it can be omitted): + +```julia +Optim.set_objective_fn(s, rosenbrock_objective_fn) +Optim.set_grad_fn(s, rosenbrock_grad_fn) +``` + +## Running the minimization + +Once the problem is fully configured, the minimization can be invoked: + +```julia +status, message = Optim.minimize(s) +x = s.x + +println("Message: ", message) +@assert status == 0 +println("x = ", x) +``` + +Note that we use the convention that the `status` is equal to zero, +when the used solver reports successful termination. +Additional return argument `message` is used to pass a free-form message, +that different solver use to define more details about the minimization +(for example, they can report only local convergence, +not the global one). + +A successful run finds a vector whose entries are all near `1.0`, which is the global minimum of the Rosenbrock function. + +## Running the full script + +The complete example is available at `examples/lang_julia/call_optim_rosenbrock.jl`. +Assuming the project has been built and the environment is loaded via `source env.sh`, run it from the project root with: + +```shell +julia examples/lang_julia/call_optim_rosenbrock.jl [implementation] [method] [linesearch] +``` + +Supported combinations of the command-line arguments for this script are: + +- `optim_jl NelderMead` +- `optim_jl BFGS StrongWolfe` +- `scipy_optimize NelderMead` +- `scipy_optimize BFGS` + +The line-search argument is only used for the Julia `BFGS` implementation. The script prints the chosen backend, the method, the solver message, and the optimizer output, and then verifies that the result is close to the true minimizer. + +[rosenbrock-wiki]: https://en.wikipedia.org/wiki/Rosenbrock_function diff --git a/docs/source/getting-started/example-optim-rosenbrock-python.md b/docs/source/getting-started/example-optim-rosenbrock-python.md new file mode 100644 index 00000000..73788287 --- /dev/null +++ b/docs/source/getting-started/example-optim-rosenbrock-python.md @@ -0,0 +1,175 @@ +# Example: Minimize the Rosenbrock function using Python + +In this example, we use _Open Interfaces_ from Python to minimize the +[Rosenbrock function][rosenbrock-wiki]. For a vector $x \in \mathbb{R}^n$, +the objective function used by the example is + +$$ + f(x) = \sum_{i=1}^{n-1} \left(a (x_{i+1} - x_i^2)^2 + + (1 - x_i)^2\right), +$$ + +where $a = 10$. The global minimizer is +$x = [1, 1, \ldots, 1]$, where $f(x) = 0$. + +## Objective and gradient functions + +The optimization interface calls the objective function with the current +vector and the user data. The example uses the user data to pass the parameter +$a$ without hard-coding it in the callback: + +```python +def rosenbrock_objective_fn(x, a): + return np.sum(a * (x[1:] - x[:-1] ** 2.0) ** 2.0 + + (1.0 - x[:-1]) ** 2.0) +``` + +The gradient callback receives the current vector, an output array, and the +same user data. It fills the output array in place and returns `0` to indicate +successful computation: + +```python +def rosenbrock_grad_fn(x, grad_f, a): + xi = x[:-1] + xip1 = x[1:] + + grad_f[:-1] = -4.0 * a * xi * (xip1 - xi**2.0) - 2.0 * (1.0 - xi) + grad_f[-1] = 0.0 + grad_f[1:] += 2.0 * a * (xip1 - xi**2.0) + return 0 +``` + +The callback signatures are defined by the interface so that the same +optimization problem can be used by implementations written in different +languages. + +## Configuring the optimizer + +After importing NumPy and the optimization interface, define an initial guess +and create an optimizer. The example starts from some random initial guess +("random", but consisting of approximations of some constants): + +```python +import numpy as np +from openinterfaces.interfaces.optim import Optim + +x0 = np.array([3.14, 2.72, 6.18, 9.81, 8.31]) +# Instantiate a Gateway to the `optim` interface +s = Optim("optim_jl") +s.set_initial_guess(x0) +``` + +Because our objective function has the user-defined parameter `a`, +we need to pass it as well: +```python +user_data = 10.0 + +s.set_user_data(user_data) +``` + + +## Passing callbacks to the objective and gradient functions + +```python +s.set_objective_fn(rosenbrock_objective_fn) +s.set_grad_fn(rosenbrock_grad_fn) +``` + +The `optim_jl` implementation uses optimization algorithms from the Julia +ecosystem. +For Nelder-Mead, we can find in the documentation +at , +that the absolute function tolerance is controlled +by the parameter `g_abstol` (or `g_tol`), +hence we define the termination criterion of the optimization process +as follows: + +```python +s.set_method("NelderMead", {"g_abstol": 1e-11}) +``` + +To use the BFGS method, +select the method, setting the termination criterion +to be the inf-norm of the gradient, +also setting the line-search strategy to `StrongWolfe`: + +```python +s.set_method("BFGS", {"g_abstol": 1e-8, "linesearch": "StrongWolfe"}) +``` + +Nelder-Mead is a gradient-free method, +so one does not need to pass the gradient callback. + +For the BFGS, we pass the gradient callback +to improve the stability by using analytic gradient +instead of finite differences: +```python +s.set_grad_fn(rosenbrock_grad_fn) +``` + +So the full set of required calls to the `Optim` +interface to define an optimization problem +for, say, the gradient-based BFGS method +from the `Optim.jl` package is: + +```python +s = Optim("optim_jl") +s.set_initial_guess(x0) +s.set_user_data(user_data) +s.set_method("BFGS", {"g_abstol": 1e-8, "linesearch": "StrongWolfe"}) +s.set_objective_fn(rosenbrock_objective_fn) +s.set_grad_fn(rosenbrock_grad_fn) +``` + +## Minimizing + +Now as the optimization problem is fully defined, +we are ready to call `minimize`. +Calling `minimize()` returns a status code and a message. +By convention, status code is equal to zero, +when the implementation reports a successful termination +of the optimization process. +The message is a string, which can be useful for the user +to see, as different solvers use free-form messages +to describe the status of optimization process +(for example, they can report that only local convergence is achieved, +but not the global one). + +The resultant minimizer is available as the property `x` +of the `Optim` gateway: +```python +status, message = s.minimize() +assert status == 0 + +print(f"Message: {message}") +print(f"x = {s.x}") +``` + +For a successful run, every component of `s.x` is close to `1` +for the Rosenbrock function used in this example. + +## Running the complete example + +The full example is available in +`examples/lang_python/call_optim_rosenbrock.py`. +Assuming that _Open Interfaces_ are successfully built, +and the environment variables are set via `source env.sh`, +we can call the example from the root directory: +```shell +python examples/lang_python/call_optim_rosenbrock.py [implementation] [method] +[linesearch] +``` + +Supported combinations include: + + - `optim_jl NelderMead` + - `optim_jl BFGS StrongWolfe` + - `scipy_optimize NelderMead` + - `scipy_optimize BFGS` + +The line-search argument is used only by `optim_jl` with BFGS. The script +prints the implementation, method, solver message, and computed minimizer, +then checks that the result is close to the exact solution. + + +[rosenbrock-wiki]: https://en.wikipedia.org/wiki/Rosenbrock_function diff --git a/docs/source/getting-started/index.rst b/docs/source/getting-started/index.rst index 7fe85dfc..e60a79dc 100644 --- a/docs/source/getting-started/index.rst +++ b/docs/source/getting-started/index.rst @@ -9,3 +9,5 @@ Getting Started example-solve-vdp-python example-solve-vdp-c example-ivp-burgers-equation-in-python + example-optim-rosenbrock-python + example-optim-rosenbrock-julia