# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# + [markdown] toc=true
# <h1>Table of Contents<span class="tocSkip"></span></h1>
# <div class="toc" style="margin-top: 1em;"><ul class="toc-item"><li><ul class="toc-item"><li><span><a href="#(a)" data-toc-modified-id="(a)-0.1"><span class="toc-item-num">0.1&nbsp;&nbsp;</span>(a)</a></span></li><li><span><a href="#(b)" data-toc-modified-id="(b)-0.2"><span class="toc-item-num">0.2&nbsp;&nbsp;</span>(b)</a></span></li><li><span><a href="#(c)" data-toc-modified-id="(c)-0.3"><span class="toc-item-num">0.3&nbsp;&nbsp;</span>(c)</a></span></li><li><span><a href="#(d)" data-toc-modified-id="(d)-0.4"><span class="toc-item-num">0.4&nbsp;&nbsp;</span>(d)</a></span></li></ul></li></ul></div>
# -

# Use the method of Steepest Descent with $TOL = 0.05$ to approximate the solutions of the following nonlinear systems.

# The nonlinear system continuous with ex10.3.3

# +
import numpy as np
from numpy import linalg
from abc import abstractmethod
import pandas as pd
import math

pd.options.display.float_format = '{:,.8f}'.format
np.set_printoptions(suppress=True, precision=8)

TOR = pow(10.0, -6)
MAX_ITR = 150


# -

class SteepestDescentMethod(object):

    def __init__(self):
        return

    @abstractmethod
    def f(self, x):
        return NotImplementedError('Implement f()!')

    @abstractmethod
    def g(self, x):
        return NotImplementedError('Implement g()!')

    @abstractmethod
    def grad_g(self, x):
        return NotImplementedError('Implement grad_g()!')

    @abstractmethod
    def jacobian(self, x):
        return NotImplementedError('Implement jacobian()!')

    @abstractmethod
    def run(self, x):
        return NotImplementedError('Implement run()!')


# ## (a) 
# $$\begin{align*}
# 4x_1^2 - 20 x_1 + \frac{1}{4} x_2^2 + 8 &= 0 \\
# \frac{1}{2}x_1x_2^2+2x_1-5x_2+8 &= 0
# \end{align*}$$

class SteepestDescent(SteepestDescentMethod):

    def __init__(self):
        super(SteepestDescentMethod, self).__init__()

    def f(self, x):
        sol = np.zeros(len(x))
        sol[0] = 4 * pow(x[0], 2) - 20 * x[0] + pow(x[1], 2) / 4 + 8
        sol[1] = x[0] * pow(x[1], 2) / 2 + 2 * x[0] - 5 * x[1] + 8
        return sol

    def jacobian(self, x):
        jac = np.zeros(shape=(2, 2))
        jac[0][0] = 8 * x[0] - 20
        jac[0][1] = x[1] / 2
        jac[1][0] = pow(x[1], 2) / 2
        jac[1][1] = x[0] * x[1] - 5
        return jac

    def g(self, x):
        sol = self.f(x)
        return sum([e * e for e in sol])

    def grad_g(self, x):
        return 2 * self.jacobian(x).transpose().dot(self.f(x))

    def run(self, x):
        df = pd.DataFrame(columns=['x' + str(i + 1) for i in range(len(x))] + ['g', 'residual', 'actual-residual'])

        row = len(df)
        df.loc[row] = [xe for xe in x] + [self.g(x), np.nan, np.nan]

        while True:
            prev_x = x
            g1 = self.g(x)
            z = self.grad_g(x)
            z0 = linalg.norm(z, 2)
            if z0 == 0.0:
                print('Zero gradient')
                return x

            z /= z0
            alpha3 = 1
            g3 = self.g(x - alpha3 * z)
            while g3 >= g1:
                alpha3 /= 2.0
                g3 = self.g(x - alpha3 * z)
                if alpha3 < TOR / 2.0:
                    print('No likely improvement')
                    return x

            alpha2 = alpha3 / 2.0
            g2 = self.g(x - alpha2 * z)

            h1 = (g2 - g1) / alpha2
            h2 = (g3 - g2) / (alpha3 - alpha2)
            h3 = (h2 - h1) / alpha3

            alpha0 = (alpha2 - h1 / h3) / 2.0
            g0 = self.g(x - alpha0 * z)

            alpha = alpha0
            g = g0
            if g3 < g:
                alpha = alpha3
                g = g3

            x = x - alpha * z
            residual = linalg.norm(x - prev_x, np.inf)
            row = len(df)
            df.loc[row] = [nxe for nxe in x] + [g, residual, np.nan]
            if math.fabs(g - g1) < TOR:
                break

        for i in range(len(df)):
            xk = np.array([df.loc[i][j] for j in range(len(x))])
            df.loc[i][4] = linalg.norm(xk - x, np.inf)

        print(self.f(x))
        return df


x0 = np.array([0, 0])
SteepestDescent().run(x0).astype(np.float64)

x0 = np.array([5, 5])
SteepestDescent().run(x0).astype(np.float64)


# ## (b) 
# $$\begin{align*}
# \sin(4\pi x_1 x_2)-2x_2-x_1 &= 0 \\
# \left( \frac{4\pi-1}{4\pi} \right)(\exp^{2x_1}-\exp)+4\exp x_2^2 -2\exp x_1 &= 0
# \end{align*}$$

class SteepestDescent(SteepestDescentMethod):

    def __init__(self):
        super(SteepestDescentMethod, self).__init__()

    def f(self, x):
        sol = np.zeros(len(x))
        sol[0] = math.sin(4 * math.pi * x[0] * x[1]) - 2 * x[1] - x[0]
        sol[1] = ((4 * math.pi - 1) / (4 * math.pi)) * (math.exp(2 * x[0]) - math.e) + 4 * math.e * pow(x[1], 2) - 2 * math.e * x[0]
        return sol

    def jacobian(self, x):
        jac = np.zeros(shape=(2, 2))
        jac[0][0] = 4 * math.pi * x[1] * math.cos(4 * math.pi * x[0] * x[1]) - 1
        jac[0][1] = 4 * math.pi * x[0] * math.cos(4 * math.pi * x[0] * x[1]) - 2
        jac[1][0] = 2 * (4 * math.pi - 1) / (4 * math.pi) * math.exp(2 * x[0]) - 2 * math.e
        jac[1][1] = 8 * math.e * x[1]
        return jac

    def g(self, x):
        sol = self.f(x)
        return sum([e * e for e in sol])

    def grad_g(self, x):
        return 2 * self.jacobian(x).transpose().dot(self.f(x))

    def run(self, x):
        df = pd.DataFrame(columns=['x' + str(i + 1) for i in range(len(x))] + ['g', 'residual', 'actual-residual'])

        row = len(df)
        df.loc[row] = [xe for xe in x] + [self.g(x), np.nan, np.nan]

        while True:
            prev_x = x
            g1 = self.g(x)
            z = self.grad_g(x)
            z0 = linalg.norm(z, 2)
            if z0 == 0.0:
                print('Zero gradient')
                return x

            z /= z0
            alpha3 = 1
            g3 = self.g(x - alpha3 * z)
            while g3 >= g1:
                alpha3 /= 2.0
                g3 = self.g(x - alpha3 * z)
                if alpha3 < TOR / 2.0:
                    print('No likely improvement')
                    return x

            alpha2 = alpha3 / 2.0
            g2 = self.g(x - alpha2 * z)

            h1 = (g2 - g1) / alpha2
            h2 = (g3 - g2) / (alpha3 - alpha2)
            h3 = (h2 - h1) / alpha3

            alpha0 = (alpha2 - h1 / h3) / 2.0
            g0 = self.g(x - alpha0 * z)

            alpha = alpha0
            g = g0
            if g3 < g:
                alpha = alpha3
                g = g3

            x = x - alpha * z
            residual = linalg.norm(x - prev_x, np.inf)
            row = len(df)
            df.loc[row] = [nxe for nxe in x] + [g, residual, np.nan]
            if math.fabs(g - g1) < TOR:
                break

        for i in range(len(df)):
            xk = np.array([df.loc[i][j] for j in range(len(x))])
            df.loc[i][4] = linalg.norm(xk - x, np.inf)

        print(self.f(x))
        return df


x0 = np.array([0, 0])
SteepestDescent().run(x0).astype(np.float64)

x0 = np.array([-10, -10])
SteepestDescent().run(x0).astype(np.float64)

x0 = np.array([1, -10])
SteepestDescent().run(x0).astype(np.float64)

x0 = np.array([-1, 1])
SteepestDescent().run(x0).astype(np.float64)


# ## (c) 
# $$\begin{align*}
# x_1(1-x_1)+4x_2 &= 12 \\
# (x_1 - x_2)^2 + (2x_2-3)^2 &= 25
# \end{align*}$$

class SteepestDescent(SteepestDescentMethod):

    def __init__(self):
        super(SteepestDescentMethod, self).__init__()

    def f(self, x):
        sol = np.zeros(len(x))
        sol[0] = x[0] * (1 - x[0]) + 4 * x[1] - 12
        sol[1] = pow(x[0] - x[1], 2) + pow(2 * x[1] - 3, 2) - 25
        return sol

    def jacobian(self, x):
        jac = np.zeros(shape=(2, 2))
        jac[0][0] = 1 - 2 * x[0]
        jac[0][1] = 4
        jac[1][0] = 2 * (x[0] - x[1])
        jac[1][1] = -2 * (x[0] * x[1]) + 4 * (2 * x[1] - 3)
        return jac

    def g(self, x):
        sol = self.f(x)
        return sum([e * e for e in sol])

    def grad_g(self, x):
        return 2 * self.jacobian(x).transpose().dot(self.f(x))

    def run(self, x):
        df = pd.DataFrame(columns=['x' + str(i + 1) for i in range(len(x))] + ['g', 'residual', 'actual-residual'])

        row = len(df)
        df.loc[row] = [xe for xe in x] + [self.g(x), np.nan, np.nan]

        while True:
            prev_x = x
            g1 = self.g(x)
            z = self.grad_g(x)
            z0 = linalg.norm(z, 2)
            if z0 == 0.0:
                print('Zero gradient')
                return x

            z /= z0
            alpha3 = 1
            g3 = self.g(x - alpha3 * z)
            while g3 >= g1:
                alpha3 /= 2.0
                g3 = self.g(x - alpha3 * z)
                if alpha3 < TOR / 2.0:
                    print('No likely improvement')
                    return x

            alpha2 = alpha3 / 2.0
            g2 = self.g(x - alpha2 * z)

            h1 = (g2 - g1) / alpha2
            h2 = (g3 - g2) / (alpha3 - alpha2)
            h3 = (h2 - h1) / alpha3

            alpha0 = (alpha2 - h1 / h3) / 2.0
            g0 = self.g(x - alpha0 * z)

            alpha = alpha0
            g = g0
            if g3 < g:
                alpha = alpha3
                g = g3

            x = x - alpha * z
            residual = linalg.norm(x - prev_x, np.inf)
            row = len(df)
            df.loc[row] = [nxe for nxe in x] + [g, residual, np.nan]
            if math.fabs(g - g1) < TOR:
                break

        for i in range(len(df)):
            xk = np.array([df.loc[i][j] for j in range(len(x))])
            df.loc[i][4] = linalg.norm(xk - x, np.inf)

        print(self.f(x))
        return df


x0 = np.array([0, 0])
SteepestDescent().run(x0).astype(np.float64)

x0 = np.array([1, 5])
SteepestDescent().run(x0).astype(np.float64)

x0 = np.array([-1, 500])
SteepestDescent().run(x0).astype(np.float64)


# ## (d) 
# $$\begin{align*}
# 5x_1^2-x_2^2 &= 0 \\
# x_2 - 0.25(\sin x_1 + \cos x_2) &= 0
# \end{align*}$$

class SteepestDescent(SteepestDescentMethod):

    def __init__(self):
        super(SteepestDescentMethod, self).__init__()

    def f(self, x):
        sol = np.zeros(len(x))
        sol[0] = 5 * pow(x[0], 2) - pow(x[1], 2)
        sol[1] = x[1] - 0.25 * (math.sin(x[0]) + math.cos(x[1]))
        return sol

    def jacobian(self, x):
        jac = np.zeros(shape=(2, 2))
        jac[0][0] = 10 * x[0]
        jac[0][1] = -2 * x[1]
        jac[1][0] = -0.25 * math.cos(x[0])
        jac[1][1] = 1 + 0.25 * math.sin(x[1])
        return jac

    def g(self, x):
        sol = self.f(x)
        return sum([e * e for e in sol])

    def grad_g(self, x):
        return 2 * self.jacobian(x).transpose().dot(self.f(x))

    def run(self, x):
        df = pd.DataFrame(columns=['x' + str(i + 1) for i in range(len(x))] + ['g', 'residual', 'actual-residual'])

        row = len(df)
        df.loc[row] = [xe for xe in x] + [self.g(x), np.nan, np.nan]

        while True:
            prev_x = x
            g1 = self.g(x)
            z = self.grad_g(x)
            z0 = linalg.norm(z, 2)
            if z0 == 0.0:
                print('Zero gradient')
                return x

            z /= z0
            alpha3 = 1
            g3 = self.g(x - alpha3 * z)
            while g3 >= g1:
                alpha3 /= 2.0
                g3 = self.g(x - alpha3 * z)
                if alpha3 < TOR / 2.0:
                    print('No likely improvement')
                    return x

            alpha2 = alpha3 / 2.0
            g2 = self.g(x - alpha2 * z)

            h1 = (g2 - g1) / alpha2
            h2 = (g3 - g2) / (alpha3 - alpha2)
            h3 = (h2 - h1) / alpha3

            alpha0 = (alpha2 - h1 / h3) / 2.0
            g0 = self.g(x - alpha0 * z)

            alpha = alpha0
            g = g0
            if g3 < g:
                alpha = alpha3
                g = g3

            x = x - alpha * z
            residual = linalg.norm(x - prev_x, np.inf)
            row = len(df)
            df.loc[row] = [nxe for nxe in x] + [g, residual, np.nan]
            if math.fabs(g - g1) < TOR:
                break

        for i in range(len(df)):
            xk = np.array([df.loc[i][j] for j in range(len(x))])
            df.loc[i][4] = linalg.norm(xk - x, np.inf)

        print(self.f(x))
        return df


x0 = np.array([0, 0])
SteepestDescent().run(x0).astype(np.float64)

x0 = np.array([1, 1])
SteepestDescent().run(x0).astype(np.float64)

x0 = np.array([-1, -1])
SteepestDescent().run(x0).astype(np.float64)

x0 = np.array([100, -100])
SteepestDescent().run(x0).astype(np.float64)

StarcoderdataJupyterScriptsDedupFiltered