import numpy as np
import matplotlib.pyplot as plt

"""
Explanation: Optimization in 1-D and N-D. Debugging Functions
Unit 11, Lecture 1
Numerical Methods and Statistics

Prof. Andrew White,  April 19, 2020
Goals:

Learn the meaning of root-finding and minimization, the two types of optimization
Understand the iterative nature of these methods
Debug common problems when defining functions, which is essential for optimization
Be able to identify convex problems and understand their complexties
Learn the two standard methods for minimize and root-finding and how to call them in Python
End of explanation
"""


from scipy.optimize import newton

root =  newton(lambda x: np.cos(x) - x, x0=0)
print(root, np.cos(root) - root)

"""
Explanation: Numerical Optimization - Root Finding
What is $x$ in
$$\cos (x) = x$$
To use a root-finding method, we must make our equation have one side be 0.
$$\cos (x) - x = 0$$
Newton's Method - For finding Roots
End of explanation
"""


root =  newton(lambda x: np.cos(x) - x, x0=0, tol=1e-3)
print(root, np.cos(root) - root)

root =  newton(lambda x: np.cos(x) - x, x0=0, tol=1e-3, maxiter=5)
print(root, np.cos(root) - root)

"""
Explanation: $$x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}$$
This method, like all this unit are iterative. This is modifiable by you, either through choosing tolerance, or maximum iterations. The functions have sensible defaults
End of explanation
"""


#Note: this code is intentionally really bad
from scipy.integrate import quad
def equation(x):
    return quad(lambda s: np.exp(-s**2), 0, x)
root = newton(equation, x0=0)

"""
Explanation: Newton's Method
Use root instead, this is just for instructional purposes. root is better
Type: Root Finding
Discrete/Continuous: Continuous
Dimensions: 1
Derivative: optional
Non-Convex: not recommended
Python: newton
Let's try to solve a more interesting equation with Newton's method!
$$ \int_0^x e^{-s^2}\,ds = \frac{1}{2}$$
Find $x$
End of explanation
"""


def version_1(x):
    return x**2 - 3 *x + 2

version_2 = lambda x: x ** 2 - 3 *x + 2

np_version = np.vectorize(version_1)

print(version_1(3.))
print(version_2(3.))

some_threes = np.zeros(10)
some_threes[:] = 3.0 # -> Notice how we don't replace the numpy array with 3, but instead make all elements in it equal to three
#some_threes = 3 -> This would delete the numpy array and now some_threes would be a single 3

print(np_version(some_threes))

"""
Explanation: Steps for Debugging Code

Correct any Python errors
Restart the Kernel
Correct any Python errors
Print to discover logical errors

Most Common Mistakes

Copy-Pasta leaving in old variables
Mixing Numpy arrays and for loops
Forgetting order of function args and return values

Review of Functions
For functions like quad, newton, and minimize the function should take in 1 argument and return 1 argument.
Correct Examples
$$x^2 - 3x + 2$$
End of explanation
"""


from scipy.integrate import quad
import numpy as np
from math import pi

def integrate_sin2(x):
    ans, err = quad(lambda s: np.sin(s) ** 2, pi, x)
    return ans
    
print(integrate_sin2(2 * pi))

"""
Explanation: $$\int_{\pi}^x \sin^2(s)\,ds$$
End of explanation
"""


from math import sqrt

def distance(x, y):
    sum_sq = np.sum((x - y)**2)
    return sqrt(sum_sq)

def my_distance(s):
    x = np.zeros(3)
    y = np.zeros(3)
    x[0] = 2.0
    x[1] = s
    x[2] = -3.5
    y[0] = 1.0
    y[1] = -3.0
    y[2] = 0.0
    return distance(x, y)

print(my_distance(1))

"""
Explanation: 2 Vectors, $\vec{x}$ and $\vec{y}$, where one component of $\vec{x}$ is changing and we want to know the distance between the two vectors.
End of explanation
"""


def version_1(x):
    x**2 - 3 *x + 2

print(version_1(3.))

"""
Explanation: Incorrect Examples
No Return Value
End of explanation
"""


def integrate_sin2(x):
    return quad(lambda s: np.sin(s) ** 2, pi, x)

print(integrate_sin2(2 * pi))

"""
Explanation: Bad Return Value
End of explanation
"""


def distance(x, y):
    sum_sq = np.sum((x - y)**2)
    return sqrt(sum_sq)

print(distance(1))

"""
Explanation: Too many arguments
End of explanation
"""


#note still wrong
def equation(x):
    ans, err = quad(lambda s: np.exp(-s**2), 0, x)
    return ans
root = newton(equation, x0=0)
print(root)

equation(root)

"""
Explanation: Let's return to our example from above:
$$ \int_0^x e^{-s^2}\,ds = \frac{1}{2}$$
Find $x$
End of explanation
"""


root = newton(lambda x: equation(x) - 0.5, x0=0)
print(root, equation(root))

"""
Explanation: We forgot to rearrange the equation to be equal to $0$
End of explanation
"""


x = 4
y = 2

#Right now, there is x,y and all other functions/variables defined or imported above in the scope

"""
Explanation: Scope
Scope means the set of variables and functions which are defined and accessible in your code
End of explanation
"""


x = 4
y = 2
#Here, I have x and y in scope
def scope_example():
    z = 2
    #Here, I have x,y and z in scope
    
#Here I again have only x and y in scope

x = 4 
y = 2
print(y,"Before function")
def scope_example():
    y = 25 #This is a new version of y that exists only in this scope
    print(y, "Inside function")
scope_example()
print(y, "After Function")

x = 4 
y = [2]
print(y,"Before function")
def scope_example():
    y[0] = 25 #Here I'm not creating a y, but modifying y
    print(y, "Inside function")
scope_example()
print(y, "After Function")

"""
Explanation: Scopes nest
End of explanation
"""


def two_well(x):
    if x < 0.125:
        return (x + 2) ** 2
    if x >= 0.125:
        return (x - 2) ** 2 + 1
    
np_two_well = np.vectorize(two_well)

x = np.linspace(-4, 4, 1000)
plt.plot(x, np_two_well(x))
plt.show()

"""
Explanation: Things to remember about scope:

Scopes nest, so that you can access things above your scope
You can modify variables from any scope you can see, but ones you create disappear outside of the scope

Returning to Optimization
Applications of Optimization

Solving non-linear equations
Solving systems of equations
Optimizing equations with or without constraints
Fitting models to data

Choosing which method to use
There are five things to consider when doing optimization:

Is it 1 or N dimensions?
Are you minimizing or root-finding?
Are there constraints?
Are there bounds?
Is it convex or non-convex?

Identifying Convexity
If a problem has more than one minimum (derivative is 0), then the problem is non-convex. The opposite of convex is non-convex. A concave function can be made convex with a negative sign. A non-convex function cannot be made convex.
Knowing about convexity can come from:

Plots
Running convex optimization in two starting positions and getting different minimums
Knowing something specific about the problem

Consider a function with two minimums:
End of explanation
"""


def obj(x):
    return x**2

x = np.linspace(-1,1,100)
plt.plot(x, obj(x))
plt.show()

from scipy.optimize import minimize

minimize(obj, x0=3)

"""
Explanation: BFGS Optimization - Minimization
Broyden–Fletcher–Goldfarb–Shanno
BFGS
Type: Minimization
Discrete/Continuous: Continuous
Dimensions: N
Derivative: optional
Non-Convex: not recommended
Python: minimize if no constraints or bounds are given
Nomenclature
In optimization, you have a function called the objective function. That's what you're minimizing. It always returns a single value. It is sometimes called the fit, the error, the residual, or the penalty.
End of explanation
"""


def f(x):
    return ((x-4)**2)/2+((x-2)**2)/4
x = np.linspace(-10, 10, 100)
plt.plot(x, f(x))
plt.show()

result = minimize(f, x0=0)
print(f'the minimum value occurs at {result.x} and is {result.fun}')

"""
Explanation: Minimize Return value
     fun: The value of the function at the minimum
hess_inv: The inverse of the the Hessian
     jac: The value of the Jacobian
 message: A string describing what happened
    nfev: Number of function evaluations
     nit: Number of iterations of the x point
    njev: Number of times it computed the Jacobian
  status: The single digit message (0 = success, != 0 some error)
 success: Boolean indicating success
       x: The minimum x

Objective Functions
Be careful that your objective function is convex and it's minimum isn't at $\infty$ or $-\infty$. 
A good objective function
Minimize the following:
$$f(x) = \frac{(x - 4)^2}{2} + \frac{(x - 2)^2}{4}$$
End of explanation
"""


result = minimize(lambda x: (x - 4) / 2, x0=0)
print(f'the minimum value occurs at {result.x} and is {result.fun}')

"""
Explanation: A bad objective function
Minimize the following:
$$
f(x) = \frac{(x - 4)} { 2 }
$$
The minimum is at $-\infty$!
End of explanation
"""


result = minimize(lambda r: 4 * (r**(-12) - r**(-6)), x0=0)
print(f'the minimum is at {result.x}')

"""
Explanation: A good objective function but bad x0
Minimize the following function
$$4 \left[ r^{-12} - r^{-6}\right]$$
End of explanation
"""


r = np.linspace(0.9, 2, 100)
y = 4 * (r**(-12) - r**(-6))
plt.plot(r, y)
plt.show()

result = minimize(lambda r: 4 * (r**(-12) - r**(-6)), x0=1)
print(f'the minimum is at {result.x} and its value is {result.fun}')

"""
Explanation: Our initial value was not in the domain!
End of explanation
"""


#place - sign to make it a maxmimization problem
result = minimize(lambda x: (x - np.cos(x))**2, x0=1)
print(f'the maximum is at {result.x}')

"""
Explanation: Maximizing
In order to maximzie a function, just add a minus sign
Maximize the following function:
$$
-\left[x - \cos(x)\right]^2
$$
End of explanation
"""


result = minimize(lambda x: (x[0] - 4)**2 / 2 + (x[1] + 3)**2 / 5, x0=[0,0])
print(f'maximum occurs when x = {result.x[0]} and y = {result.x[1]}')

"""
Explanation: Multiple Dimensions
Just indicate multiple dimensions by using a multidimensional x0. Note that your $x$ becomes a vector!
Minimize the following:
$$
f(x, y) = \frac{(x - 4)^2} { 2 } + \frac{(y + 3)^2} { 5 }
$$
End of explanation
"""


x = np.linspace(-1, 1, 1000)
plt.plot(x,  2 * x ** 3 - 0 * x **2  - x)
plt.show()

from scipy.optimize import minimize
minimize(lambda x:  2 * x ** 3 - 0 * x **2  - x, x0=0.05)

x = np.linspace(-1, 1, 1000)
plt.plot(x,  2 * x ** 3 - 0 * x **2  - x)
plt.show()

from scipy.optimize import minimize
minimize(lambda x:  2 * x ** 3 - 0 * x **2  - x, x0=-0.5)

"""
Explanation: Identifying Convexity for Minimization - Example
The best ways to identify convexity are:

plot it
try optimizing in multiple starting positions
End of explanation
"""


from scipy.optimize import root
#rearranged equation so all terms on one side
result = root(lambda x: np.cos(x) + np.sin(x) - x, x0=1)
print(result)

"""
Explanation: Powell hybrid method
Type: Root finding
Discrete/Continuous: Continuous
Dimensions: N
Derivative: optional
Non-Convex: not recommended
Python: root unless method argument specifies a different method
1D Example
This is exactly like newton from above. Solve this equation:
$$
\cos x + \sin x = x
$$
End of explanation
"""


x = result.x 
print(np.cos(x) + np.sin(x) - x)

"""
Explanation: The result type is like what we saw for minimize. Similar terms are here, including the root and the value of the function at the root. Notice it's not exactly $0$ at the root.
End of explanation
"""


def sys(v):
    #I'm using v here to distinguish from the x in the equations
    #extract x and y
    x = v[0]
    y = v[1]
    #compute equations
    eq1 = 3 * x ** 2 - 2 * y - 4
    eq2 = x - 4 * y**2 + 2
    #pack into list
    sys_eq = [eq1, eq2]
    return sys_eq
root(sys, x0=[1,1])

"""
Explanation: Solve the following system of equations:
$$ 3 x^2 - 2 y = 4$$
$$ x - 4 y ^ 2 = -2$$
End of explanation
"""


def sys2(v):
    #I'm using v here to distinguish from the x in the equations
    #extract x and y
    x = v[0]
    y = v[1]
    #compute equations
    eq1 = 3 * x ** 2 - 2 * y**2 - 4
    eq2 = x**3 - 4 * y + 2
    #pack into list
    sys_eq = [eq1, eq2]
    return sys_eq
root(sys2, x0=[1,1])

"""
Explanation: So the answer is $x = 1.40,\: y = 0.921$
You should not ignore the information in the output of root. Imagine this small modification:
$$ 3 x^2 - 2 y^2 = 4$$
$$ x^3 - 4 y = -2$$
End of explanation
"""


root(lambda x: [3 * x[0] ** 2 - 2 * x[1] - 4, x[0] - 4 * x[1] ** 2 + 2], x0=[1,1])

"""
Explanation: If you did not read the message or status, you might believe the method succeeded. 
More compact version:
End of explanation
"""


root(lambda x: [3 * x[0] ** 2 - 2 * x[1] - 4, x[0] - 4 * x[1] ** 2 + 2], x0=[0, 0])

"""
Explanation: Importance of starting position
By the way, there are two roots to this function!
End of explanation
"""


