straw/pm21-dragon

Fork 23

Andrew Straw b6b4acaa36 lectures as finished

2024-12-02 10:47:50 +01:00

66 KiB

Raw Blame History

None <html> <head> </head>

Optimization¶

Let's consider the following code:

In [4]:

import numpy as np
import matplotlib.pyplot as plt

# position of flower in 2D
flower = np.array([29.1, 19.1])

def make_bee_track(t):
    """Find bee position at time t"""
    pos0 = (13.21, 4.56)
    velocity = (3.1, 0.25)
    pos_x = pos0[0] + t*velocity[0]
    pos_y = pos0[1] + t*velocity[1]
    return np.array([pos_x,pos_y])

t = np.linspace(0,15,10)
bee_track = make_bee_track(t)

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4,3))
ax.plot( [flower[0]], [flower[1]], 'or', label='flower' )
ax.plot( bee_track[0], bee_track[1], '.-k', label='bee')
ax.axis('equal')
ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('y')

Out[4]:

Text(0, 0.5, 'y')

In the above code, we parameterized the bee trajectory by the variable t in the function make_bee_track(). This means we could get a new point on the track by choosing a new value of t. For example:

In [5]:

print(make_bee_track(0.0))
print(make_bee_track(0.1))
print(make_bee_track(0.2))
print(make_bee_track(1.0))

[13.21  4.56]
[13.52   4.585]
[13.83  4.61]
[16.31  4.81]

Now let's measure the distance between the bee and the flower.

In [8]:

def compute_distance(a,b):
    a = np.array(a)
    b = np.array(b)
    return np.sqrt(np.sum((a-b)**2))

n_time_steps = bee_track.shape[1]
distance = np.zeros(n_time_steps)
for i in range(n_time_steps):
    bee_pos = bee_track[:,i]
    distance[i] = compute_distance(bee_pos, flower)
display(distance)

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4,3))
ax.plot( t, distance )
ax.set_xlabel('t')
ax.set_ylabel('distance')

array([21.5384238 , 17.73297556, 14.79017429, 13.29572112, 13.73095975,
       15.93858282, 19.32029244, 23.37188672, 27.80202491, 32.45606569])

Out[8]:

Text(0, 0.5, 'distance')

Given the plot of distance versus t above, we can see the distance is minimized when t is near 5. What is the bee's position when t is 5?

In [9]:

print(make_bee_track(5))

[28.71  5.81]

We can check back to the xy plot to see, indeed, this point is pretty close to the flower.

What if we want to know, however, exactly the closest point? There are several ways to find this. Here we are going to use a "brute force" approach which will work on many different problems. Specifically, we will use scipy.optimize.minimize_scalar. The overall idea of this kind of numerical optimization is that we find the best fitting parameter to minimize our "error".

In this example, we are not so much concerned with the exact algorithm being used, but with the way we call this algorithm.

Let's create a function which uses the flower location to compute distance:

In [10]:

def calc_distance_func(t):
    # assume variable 'flower' in our global scope
    # calculate bee position (also assuming 'make_bee_track' in scope)
    dist = compute_distance(flower, make_bee_track(t))
    print(f't: {t} -> dist: {dist}')
    return dist

In [11]:

calc_distance_func(0)

t: 0 -> dist: 21.538423804912004

Out[11]:

21.538423804912004

In [12]:

calc_distance_func(5)

t: 5 -> dist: 13.295721116208782

Out[12]:

13.295721116208782

In [15]:

import scipy.optimize
result = scipy.optimize.minimize_scalar(calc_distance_func)

t: 0.0 -> dist: 21.538423804912004
t: 1.0 -> dist: 19.17780487959975
t: 2.6180339999999998 -> dist: 15.913623755962107
t: 9.502603885705089 -> dist: 18.222622377854467
t: 5.24770562096323 -> dist: 13.233470063982764
t: 6.8729320915536185 -> dist: 13.918754747176704
t: 5.437736318938615 -> dist: 13.215989304828
t: 5.47123644934229 -> dist: 13.215645882055249
t: 5.468503780784455 -> dist: 13.215643128079853
t: 5.468493134264144 -> dist: 13.2156431280385
t: 5.468493215207842 -> dist: 13.215643128038504
t: 5.468493053320446 -> dist: 13.215643128038506

In [17]:

result

Out[17]:

 message: 
          Optimization terminated successfully;
          The returned value satisfies the termination criteria
          (using xtol = 1.48e-08 )
 success: True
     fun: 13.2156431280385
       x: 5.468493134264144
     nit: 8
    nfev: 12

In [18]:

scipy.optimize.minimize_scalar?

Signature:
scipy.optimize.minimize_scalar(
    fun,
    bracket=None,
    bounds=None,
    args=(),
    method=None,
    tol=None,
    options=None,
)
Docstring:
Local minimization of scalar function of one variable.

Parameters
----------
fun : callable
    Objective function.
    Scalar function, must return a scalar.
bracket : sequence, optional
    For methods 'brent' and 'golden', `bracket` defines the bracketing
    interval and is required.
    Either a triple ``(xa, xb, xc)`` satisfying ``xa < xb < xc`` and
    ``func(xb) < func(xa) and  func(xb) < func(xc)``, or a pair
    ``(xa, xb)`` to be used as initial points for a downhill bracket search
    (see `scipy.optimize.bracket`).
    The minimizer ``res.x`` will not necessarily satisfy
    ``xa <= res.x <= xb``.
bounds : sequence, optional
    For method 'bounded', `bounds` is mandatory and must have two finite
    items corresponding to the optimization bounds.
args : tuple, optional
    Extra arguments passed to the objective function.
method : str or callable, optional
    Type of solver.  Should be one of:

        - :ref:`Brent <optimize.minimize_scalar-brent>`
        - :ref:`Bounded <optimize.minimize_scalar-bounded>`
        - :ref:`Golden <optimize.minimize_scalar-golden>`
        - custom - a callable object (added in version 0.14.0), see below

    Default is "Bounded" if bounds are provided and "Brent" otherwise.
    See the 'Notes' section for details of each solver.

tol : float, optional
    Tolerance for termination. For detailed control, use solver-specific
    options.
options : dict, optional
    A dictionary of solver options.

        maxiter : int
            Maximum number of iterations to perform.
        disp : bool
            Set to True to print convergence messages.

    See :func:`show_options()` for solver-specific options.

Returns
-------
res : OptimizeResult
    The optimization result represented as a ``OptimizeResult`` object.
    Important attributes are: ``x`` the solution array, ``success`` a
    Boolean flag indicating if the optimizer exited successfully and
    ``message`` which describes the cause of the termination. See
    `OptimizeResult` for a description of other attributes.

See also
--------
minimize : Interface to minimization algorithms for scalar multivariate
    functions
show_options : Additional options accepted by the solvers

Notes
-----
This section describes the available solvers that can be selected by the
'method' parameter. The default method is the ``"Bounded"`` Brent method if
`bounds` are passed and unbounded ``"Brent"`` otherwise.

Method :ref:`Brent <optimize.minimize_scalar-brent>` uses Brent's
algorithm [1]_ to find a local minimum.  The algorithm uses inverse
parabolic interpolation when possible to speed up convergence of
the golden section method.

Method :ref:`Golden <optimize.minimize_scalar-golden>` uses the
golden section search technique [1]_. It uses analog of the bisection
method to decrease the bracketed interval. It is usually
preferable to use the *Brent* method.

Method :ref:`Bounded <optimize.minimize_scalar-bounded>` can
perform bounded minimization [2]_ [3]_. It uses the Brent method to find a
local minimum in the interval x1 < xopt < x2.

Note that the Brent and Golden methods do not guarantee success unless a
valid ``bracket`` triple is provided. If a three-point bracket cannot be
found, consider `scipy.optimize.minimize`. Also, all methods are intended
only for local minimization. When the function of interest has more than
one local minimum, consider :ref:`global_optimization`.

**Custom minimizers**

It may be useful to pass a custom minimization method, for example
when using some library frontend to minimize_scalar. You can simply
pass a callable as the ``method`` parameter.

The callable is called as ``method(fun, args, **kwargs, **options)``
where ``kwargs`` corresponds to any other parameters passed to `minimize`
(such as `bracket`, `tol`, etc.), except the `options` dict, which has
its contents also passed as `method` parameters pair by pair.  The method
shall return an `OptimizeResult` object.

The provided `method` callable must be able to accept (and possibly ignore)
arbitrary parameters; the set of parameters accepted by `minimize` may
expand in future versions and then these parameters will be passed to
the method. You can find an example in the scipy.optimize tutorial.

.. versionadded:: 0.11.0

References
----------
.. [1] Press, W., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery.
       Numerical Recipes in C. Cambridge University Press.
.. [2] Forsythe, G.E., M. A. Malcolm, and C. B. Moler. "Computer Methods
       for Mathematical Computations." Prentice-Hall Series in Automatic
       Computation 259 (1977).
.. [3] Brent, Richard P. Algorithms for Minimization Without Derivatives.
       Courier Corporation, 2013.

Examples
--------
Consider the problem of minimizing the following function.

>>> def f(x):
...     return (x - 2) * x * (x + 2)**2

Using the *Brent* method, we find the local minimum as:

>>> from scipy.optimize import minimize_scalar
>>> res = minimize_scalar(f)
>>> res.fun
-9.9149495908

The minimizer is:

>>> res.x
1.28077640403

Using the *Bounded* method, we find a local minimum with specified
bounds as:

>>> res = minimize_scalar(f, bounds=(-3, -1), method='bounded')
>>> res.fun  # minimum
3.28365179850e-13
>>> res.x  # minimizer
-2.0000002026
File:      ~/anaconda3/envs/pm21-dragon/lib/python3.11/site-packages/scipy/optimize/_minimize.py
Type:      function

In [19]:

result

Out[19]:

 message: 
          Optimization terminated successfully;
          The returned value satisfies the termination criteria
          (using xtol = 1.48e-08 )
 success: True
     fun: 13.2156431280385
       x: 5.468493134264144
     nit: 8
    nfev: 12

In [20]:

type(result)

Out[20]:

scipy.optimize._optimize.OptimizeResult

In [21]:

result.x

Out[21]:

5.468493134264144

In [22]:

print(calc_distance_func(5))

t: 5 -> dist: 13.295721116208782
13.295721116208782

In [23]:

calc_distance_func(5.468493134264144)

t: 5.468493134264144 -> dist: 13.2156431280385

Out[23]:

13.2156431280385

In [24]:

# Where is the bee for this value of `t`?
make_bee_track(5.468493134264144)

Out[24]:

array([30.16232872,  5.92712328])

In [25]:

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4,3))
ax.plot( [flower[0]], [flower[1]], 'or', label='flower' )
ax.plot( bee_track[0], bee_track[1], '.-k', label='bee')

x,y = make_bee_track(5.468493134264144)
ax.plot( [x], [y] ,'x', label='closest')
ax.axis('equal')
ax.legend()
ax.set_xlabel('x')
ax.set_ylabel('y')

Out[25]:

Text(0, 0.5, 'y')

</html>

66 KiB Raw Blame History

Optimization¶

66 KiB

Raw Blame History