---
jupyter:
  jupytext:
    notebook_metadata_filter: all,-language_info
    split_at_heading: true
    text_representation:
      extension: .Rmd
      format_name: rmarkdown
      format_version: '1.2'
      jupytext_version: 1.14.1
  kernelspec:
    display_name: Python 3 (ipykernel)
    language: python
    name: python3
---

# Optimization

In [The Mean and Slopes](mean_and_slopes), we used a simple but slow way to
find the slope that best predicted one vector of values from another vector of
values.

First we go back to find that slope.

```{python}
import numpy as np
import matplotlib.pyplot as plt
# Make plots look a little bit more fancy
plt.style.use('fivethirtyeight')
# Print to 2 decimal places, show tiny values as 0
np.set_printoptions(precision=2, suppress=True)
import pandas as pd
pd.set_option('mode.copy_on_write', True)
```

Download the {download}`ckd_clean.csv <../data/ckd_clean.csv>` file to the
same directory as this notebook, if you are running on your own computer.

We fetch and process the data. See [mean and slopes](mean_and_slopes) for a slower description of this processing.

```{python}
# Load the data file
ckd = pd.read_csv('ckd_clean.csv')
pcv = np.array(ckd['Packed Cell Volume'])
hgb = np.array(ckd['Hemoglobin'])
```

Our criterion is the root mean squared error (RMSE).

```{python}
def calc_rmse(slope):
    predicted_pcv = hgb * slope  # 'hgb' comes from the top level
    errors = pcv - predicted_pcv # 'pcv' comes from the top level
    return np.sqrt(np.mean(errors ** 2))
```

We found the best slope by trying a very large number of slopes, and recording,
for each slope, the root mean squared error.  We chose the slope from the
slopes that we tried, that gave us the lowest root mean squared error.

```{python}
# Slopes to try
some_slopes = np.arange(2, 4, 0.01)
n_slopes = len(some_slopes)
# Try all these slopes, calculate and record RMSE.
rmse_errors = np.zeros(n_slopes)
best_slope = None
least_error = np.inf
for i in np.arange(n_slopes):
    # Get the slope.
    slope = some_slopes[i]
    # Calculate the error measure.
    this_error = calc_rmse(slope)
    # Record this slope if this is the smallest error.
    if this_error < least_error:
        least_error = this_error
        best_slope = slope
    # Put the error into the array of errors
    rmse_errors[i] = this_error

# The slope minimizing the RMSE.
best_slope
```

At the end, of [the mean and slopes](mean_and_slopes) notebook, we saw that
a function in Scipy called `minimize` can do this work for us, relatively
quickly.

```{python}
from scipy.optimize import minimize
minimize(calc_rmse, 3)
```

What we are doing, with our slow dumb technique, and with the `minimize`
function, is something called [mathematical
optimization](https://en.wikipedia.org/wiki/Mathematical_optimization).  We use
*optimization* when we have some *function* that takes one or more
*parameters*.  We want to chose, or *optimize* the parameters to give us some
desired output from the function.

In our case our *function* is the root mean squared error, `calc_rmse`.  The
*parameter* is the slope.  We are trying to find the value for the parameter
that *minimizes* the result of calling the function `calc_rmse`.

One way of doing this minimization, is the slow dumb way.  We just try a huge
number of values for the parameter (the slope), and chose the value that gives
us the lowest output value (the root mean squared error).

This is such a common problem, that there has been an enormous amount of
theoretical and practical work on building algorithms to make process of
searching for the minimum value more efficient.

This notebook is to give you an idea of how you might do this, and therefore,
the kind of things that `minimize` can do, to search quickly for the best
parameter.

Let's look again at the shape of the curve relating the slope to the root mean
squared error:

```{python}
plt.plot(some_slopes, rmse_errors)
plt.xlabel('Candidate slopes')
plt.ylabel('Root mean squared error')
plt.title('RMSE as a function of slope')
```

This is the *function* we are trying minimize.  Specifically, we are trying to
optimize the function that gives the *SSE* as a function of the *slope*
parameter.

We want to avoid trying every possible value for the slope.

To do this, we are going to start with one value for the slope, say 100, then
see if there is a good way to chose the next value to try.

Looking at the graph, we see that, when the slope is far away from the minimum,
the RMSE (on the y axis) changes very quickly as the slope changes.  That is,
the function has a steep *gradient*.

Maybe we could check what the gradient is, at our starting value of 100, by
calculating the RMSE (y) value, and then calculating the RMSE (y) value when we
increase the slope by a tiny amount.  This is the change in y for a very small
change in x.  We divide the change in y by the change in x, to get the
gradient:

```{python}
def calc_rmse_gradient(x, dx=0.0001):
    # Gradient of the RMSE at this value of x
    # RMSE at this x value.
    sos_0 = calc_rmse(x)
    # RMSE a tiny bit to the right on the x axis.
    sos_1 = calc_rmse(x + dx)
    # gradient is y difference divided by x difference.
    return (sos_1 - sos_0) / dx
```

```{python}
# The y value of the function.
calc_rmse(4)
```

```{python}
# The gradient of the function at this point.
calc_rmse_gradient(4)
```

A large positive gradient means the x value (slope) that we tried is still far
to the *right* of the minimum. This might encourage us to try an x value that
is well to the left.  We could call this a large step in x, and therefore
a large *step size*.

Let's try another value:

```{python}
# The y value of the function.
calc_rmse(2)
```

```{python}
# The gradient of the function at this point.
calc_rmse_gradient(2)
```

A large negative gradient means the x value (slope) that we tried is still far
to the *left* of the minimum.   This might encourage us to try an x value that
is well to the right.

As the gradients get small, we want to take smaller steps, so we don't miss the
minimum.

The general idea then, is to chose our step sizes in proportion to the gradient
of the function.

This is the optimization technique known as [gradient
descent](https://en.wikipedia.org/wiki/Gradient_descent).

Here it is in action. We try new x-axis values by making big jumps when the
gradient is steep, and small jumps when the gradient is shallow.  For each new
value of x, we check the gradient, and choose our new jump size.

Remember each value on the x-axis is another value we want to try, in order to
find the best-fitting slope for our original problem.

```{python}
next_x = 4 # We start the search at x=4
gamma = 0.00001 # Step size multiplier
precision = 0.00001 # Desired precision of result
max_iters = 1000 # Maximum number of iterations

for i in np.arange(max_iters):
    # Go to the next x value
    current_x = next_x
    # Estimate the gradient
    gradient = calc_rmse_gradient(current_x)
    # Use gradient to choose the next x value to try.
    # This takes negative steps when the gradient is positive
    # and positive steps when the gradient is negative.
    next_x = current_x - gamma * gradient
    step = next_x - current_x
    print('x: {:0.5f}; step {:0.5f}; gradient {:0.2f}'.format(
        current_x, step, gradient))
    # When the step size is equal to or smaller than the desired
    # precision, we are near enough.
    if abs(step) <= precision:
        # Break out of the loop.
        break

print("Minimum at", next_x)
```

As you can see, by doing this, we save ourselves from trying a very large
number of other potential x values, and we focus in on the minimum
very quickly.

This is just one method among many for optimizing our search for the minimum of
a function.  Now you know what kind of thing it is doing, we will just let
`minimize` do its job for us:

```{python}
# Use float to show us the final result in higher precision
result = minimize(calc_rmse, 100)
float(result.x)
```
