# TODO:

* It is incorrect to state that the x axis of the left column is "iteration" because those are only the accepted points in both algorithms. I don't know what else to call it. The assignment states to plot the estimate of the maximum as a function of the iteration number. I'm filtering my domain, so my plot will look like it converges much faster.

  Mention this in the report, but don't worry about fixing it.


```python
%config InlineBackend.figure_format = 'svg'
%matplotlib inline

import itertools
import multiprocessing
from multiprocessing.dummy import Pool as ThreadPool
from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
# Apparently, SNS stands for "Samuel Norman Seaborn", a fictional
# character from The West Wing
import seaborn as sns
import sympy

sns.set()
sympy.init_printing()
# Make the figures directory if it doesn't exist.
Path('figures/').mkdir(exist_ok=True)
```


```python
def f(x):
    """The function to evaluate.
    
    This function returns a sympy symbolic function, float, or np.ndarray
    depending on the type of the input.
    """
    # Use symbolic sine, pi if necessary.
    sin, pi = (sympy.sin, sympy.pi) if isinstance(x, sympy.Symbol) else (np.sin, np.pi)

    return 2 ** (-2 * ((x - 0.1) / 0.9)**2) * sin(5 * pi * x) ** 6
```


```python
x = np.linspace(0, 1, 200)
plt.plot(x, f(x))
plt.title('The objective function $f(x)$')
plt.xlabel('$x$')
plt.ylabel('$f(x)$')
plt.savefig('figures/prob1-function.pdf')
plt.show()
```

We can visually pick out the periodic extremals occuring at $\frac{n}{5} - \frac{1}{10}$, with the global optimum over $[0, 1]$ occuring at $\frac{1}{10}$.

# Symbolic and Gradient Methods

yuck


```python
x_ = sympy.Symbol('x')
f(x_)
```


```python
fp = sympy.diff(f(x_), x_)
fp
```


```python
# solve f'(x) = 0
sympy.solveset(fp, x_)
```


```python
fp = sympy.lambdify([x_], fp, 'numpy')
```


```python
plt.plot(x, fp(x))
plt.title('The derivative of $f(x)$')
plt.xlabel('$x$')
plt.ylabel('$f\'(x)$')
plt.savefig('figures/prob1-derivative.pdf')
plt.show()
```


```python
result = sp.optimize.minimize(lambda x: -f(x), 0.19, bounds=[(0, 1)])
result
```


```python
result = sp.optimize.minimize(lambda x: -f(x), 0.2, bounds=[(0, 1)])
result
```


```python
result = sp.optimize.minimize(lambda x: -f(x), 0.21, bounds=[(0, 1)])
result
```

Observations:

The (default) gradient method is just as susceptible to local solutions as the hill climbing algorithm. This makes sense.

# Hill Climbing


```python
def perturb(x, bounds, sigma):
    """Perturb the given value by adding zero-mean white noise.

    :param x: The value to perturb.
    :type x: float
    :param bounds: The lower and upper bounds on the feasible region.
    :type bounds: a (lower, upper) tuple.
    :param sigma: The standard deviation to use when adding white noise.
    """
    m, M = bounds
    xp = x + np.random.normal(scale=sigma)
    while xp >= M or xp <= m:
        xp = x + np.random.normal(scale=sigma)
    return xp
```


```python
def hill_climbing(func, bounds, sigma, iters):
    """Minimize the given function using the Hill Climbing algorithm.

    Each run generates a random initial guess.

    :param func: The function to minimize.
    :type func: a function f: float -> float
    :param bounds: The lower and upper bounds on the feasible region.
    :type bounds: a (lower, upper) tuple.
    :param sigma: The standard deviation to use when perturbing the current guess.
    :param iters: The number of iterations to run the hill climbing algorithm for.
    :returns: The path of points visited from the initial guess to the final solution.
    :rtype: An np.ndarray with shape (iters,)
    """
    path = np.zeros(iters)
    x0 = np.random.uniform(*bounds)
    for i in range(iters):
        xp = perturb(x0, bounds, sigma)
        if func(xp) < func(x0):
            x0 = xp
        path[i] = x0
    return path
```


```python
rows = 3
fig, axes = plt.subplots(rows, 2, figsize=(9, 9), sharey=True)
axes = iter(axes.flatten())
for _ in range(rows):
    path = hill_climbing(lambda x: -f(x), bounds=(0, 1), sigma=0.1, iters=100)
    
    ax = next(axes)
    
    ax.plot(path)
    ax.set_title('Hill Climbing Path')
    ax.set_xlabel('iteration')
    ax.set_ylabel('$x$')
    
    ax = next(axes)
    
    # Remove duplicate entries
    path = np.unique(path)
    ax.plot(x, f(x), label='$f(x)$')
    ax.plot(path, f(path), 'o', label='Points accepted')

    ax.set_title('Hill Climbing Path')
    ax.set_xlabel('$x$')
    ax.set_ylabel('$f(x)$')
    ax.legend()

plt.tight_layout()
plt.savefig('figures/prob1-hill-climbing-results.pdf')
plt.show()
```


```python
def iterated_hill_climbing(func, bounds, sigma, inner_iters, iters):
    """Repeatedly climb the hill to find the less-local extremum via PTSD.
    
    Runs the iterations in parallel with as many threads as there are processors.
    
    :param func: The function to minimize.
    :type func: a function f: float -> float
    :param bounds: The lower and upper bounds on the feasible region.
    :type bounds: a (lower, upper) tuple.
    :param sigma: The standard deviation to use when perturbing the current guess.
    :param inner_iters: The number of iterations to use for each run of the algorithm.
    :param iters: The number of times to run the algorithm.
    :returns: An array of solutions from each run, sorted by their fitness.
    """
    pool = ThreadPool(multiprocessing.cpu_count())
    # starmap consumes the given iterable in parallel until it is exhausted, collecting the results.
    results = pool.starmap(hill_climbing, itertools.repeat((func, bounds, sigma, inner_iters), times=iters))
    # Each result is a full path, not the optimal value
    optimums = [r[-1] for r in results]
    # Sort the optimums by their fitness.
    optimums.sort(key=func)
    return np.array(optimums)
```


```python
optimums = iterated_hill_climbing(lambda x: -f(x), bounds=(0, 1), sigma=0.1, inner_iters=50, iters=20)
solution = optimums[0]
print(solution)

plt.plot(x, f(x), label='$f(x)$')
plt.plot(solution, f(solution), 'go', label='Best solution')
plt.plot(optimums, f(optimums), '.', label='Found optimums')

plt.title('Iterated Hill Climbing Optimums')
plt.xlabel('$x$')
plt.ylabel('$f(x)$')
plt.legend()
plt.savefig('figures/prob1-hill-climbing-solution.pdf')
plt.show()
```

# Simulated Annealing


```python
def simulated_annealing(func, bounds, sigma, temp, cooling_factor):
    """Use simulated annealing to optimize the given function.

    :param func: The function to minimize.
    :type func: a function f: float -> float
    :param bounds: The lower and upper bounds on the feasible region.
    :type bounds: a (lower, upper) tuple.
    :param sigma: The standard deviation to use when perturbing the current guess.
    :param temp: The initial temperature of the system.
    :param cooling_factor: How quickly the system should cool.
    :returns: An array of points *accepted* by the random condition.
    :rtype: a 1D np.ndarray
    """
    # Pick a random starting point somewhere in the domain.
    x0 = np.random.uniform(*bounds)
    current = func(x0)
    path = []
    while temp > 0.001:
        xp = perturb(x0, bounds, sigma)
        potential = func(xp)
        if potential < current or np.random.random() < np.exp((current - potential) / temp):
            x0 = xp
            path.append(x0)
            current = potential
        temp *= 1 - cooling_factor

    return np.array(path)
```


```python
rows = 3
fig, axes = plt.subplots(rows, 2, figsize=(9, 9), sharey=True)
axes = iter(axes.flatten())
for _ in range(rows):
    path = simulated_annealing(lambda x: -f(x), bounds=(0, 1), sigma=0.1, temp=0.01, cooling_factor=0.001)
    solution = path[-1]

    ax = next(axes)

    ax.plot(path)
    ax.set_title('Simulated Annealing Path')
    ax.set_xlabel('iteration')
    ax.set_ylabel('$x$')

    ax = next(axes)

    ax.plot(x, f(x), label='$f(x)$')
    ax.plot(path, f(path), '.', label='Points accepted')
    ax.plot(solution, f(solution), 'o', label='Optimum')

    ax.set_title('Simulated Annealing Path')
    ax.set_xlabel('$x$')
    ax.set_ylabel('$f(x)$')
    ax.legend()

plt.tight_layout()
plt.savefig('figures/prob1-simulated-annealing-results.pdf')
plt.show()
```


```python
def iterated_simulated_annealing(func, bounds, sigma, temp, cooling_factor, iters):
    """Repeatedly run simulated annealing to improve the quality of the results.

    Runs the iterations in parallel with as many threads as there are processors.

    :param func: The function to minimize.
    :type func: a function f: float -> float
    :param bounds: The lower and upper bounds on the feasible region.
    :type bounds: a (lower, upper) tuple.
    :param sigma: The standard deviation to use when perturbing the current guess.
    :param temp: The initial temperature of the system.
    :param cooling_factor: How quickly the system should cool.
    :param iters: The number of times to run the algorithm.
    :returns: An array of solutions from each run, sorted by their fitness.
    """
    pool = ThreadPool(multiprocessing.cpu_count())
    # starmap consumes the given iterable in parallel until it is exhausted, collecting the results.
    results = pool.starmap(simulated_annealing, itertools.repeat((func, bounds, sigma, temp, cooling_factor), times=iters))
    # Each result is a full path, not the optimal value
    optimums = [r[-1] for r in results]
    # Sort the optimums by their fitness.
    optimums.sort(key=func)
    return np.array(optimums)
```


```python
optimums = iterated_simulated_annealing(lambda x: -f(x), bounds=(0, 1), sigma=0.1, temp=0.01, cooling_factor=0.001, iters=20)
solution = optimums[0]
print(solution)

plt.plot(x, f(x), label='$f(x)$')
plt.plot(solution, f(solution), 'go', label='Best solution')
plt.plot(optimums, f(optimums), '.', alpha=0.5, label='Found optimums')

plt.title('Iterated Simulated Annealing Optimums')
plt.xlabel('$x$')
plt.ylabel('$f(x)$')
plt.legend()
plt.savefig('figures/prob1-simulated-annealing-solutions.pdf')
plt.show()
```

Observations:

* When the simulated annealing algorithm converges, it converges closer to one of the actual peaks. The hill climbing algorithm has more noise in its convergence.
* Simulated annealing required a lot of tuning to work, and even more tuning to work well.
* The tunable parameters that affected the correctness of the solution the most were the standard deviation of the random noise and the cooling factor. The parameter that affected speed of convergence the most was the initial temperature.
* The solution clusters for the simulated annealing algorithm are much tighter than those for the hill climbing algorithm.

You want enough noise to explore, yet not so much that you jump all over the place. Temperatures ranging from 1000 to 0.002 (because 0.001 is the convergence criterion) worked well, but smaller temperatures converged much faster.

67b457bbcdab86a241af4248363e2cd5619620fc
ipynb
Jupyter Notebook
homework/hw1/prob1/prob1.ipynb
Notgnoshi/natural-computing
af6e10aaa4efcc6787bd5f71fb43fb66e67154ce
MIT
homework/hw1/prob1/prob1.ipynb
Notgnoshi/natural-computing
af6e10aaa4efcc6787bd5f71fb43fb66e67154ce
MIT
homework/hw1/prob1/prob1.ipynb
Notgnoshi/natural-computing
af6e10aaa4efcc6787bd5f71fb43fb66e67154ce
MIT
Qwen/Qwen-72B
1. YES
2. YES