diff --git a/doc/pub/week39/html/week39-bs.html b/doc/pub/week39/html/week39-bs.html index c7edb4c96..165880adb 100644 --- a/doc/pub/week39/html/week39-bs.html +++ b/doc/pub/week39/html/week39-bs.html @@ -232,10 +232,18 @@ doconce format html week39.do.txt --html_style=bootstrap --pygments_html_style=d 'the-syntax-a-dot-b-when-finding-the-dot-product'), ('Recommended to avoid', 2, None, 'recommended-to-avoid'), ('Using Autograd with OLS', 2, None, 'using-autograd-with-ols'), + ('Same code but now with momentum gradient descent', + 2, + None, + 'same-code-but-now-with-momentum-gradient-descent'), ('Including Stochastic Gradient Descent with Autograd', 2, None, 'including-stochastic-gradient-descent-with-autograd'), + ('Same code but now with momentum gradient descent', + 2, + None, + 'same-code-but-now-with-momentum-gradient-descent'), ('And Logistic Regression', 2, None, 'and-logistic-regression'), ('Introducing "JAX":"https://jax.readthedocs.io/en/latest/"', 2, @@ -319,7 +327,7 @@ MathJax.Hub.Config({
  • Program example for gradient descent with Ridge Regression
  • Using gradient descent methods, limitations
  • Improving gradient descent with momentum
  • -
  • Same code but now with momentum gradient descent
  • +
  • Same code but now with momentum gradient descent
  • Overview video on Stochastic Gradient Descent
  • Batches and mini-batches
  • Stochastic Gradient Descent (SGD)
  • @@ -352,9 +360,11 @@ MathJax.Hub.Config({
  • The syntax a.dot(b) when finding the dot product
  • Recommended to avoid
  • Using Autograd with OLS
  • -
  • Including Stochastic Gradient Descent with Autograd
  • -
  • And Logistic Regression
  • -
  • Introducing "JAX":"https://jax.readthedocs.io/en/latest/"
  • +
  • Same code but now with momentum gradient descent
  • +
  • Including Stochastic Gradient Descent with Autograd
  • +
  • Same code but now with momentum gradient descent
  • +
  • And Logistic Regression
  • +
  • Introducing "JAX":"https://jax.readthedocs.io/en/latest/"
  • @@ -384,7 +394,7 @@ MathJax.Hub.Config({
    -

    Sep 29, 2022

    +

    Sep 30, 2022


    @@ -409,7 +419,7 @@ MathJax.Hub.Config({
  • 9
  • 10
  • ...
  • -
  • 81
  • +
  • 83
  • »
  • diff --git a/doc/pub/week39/html/week39-reveal.html b/doc/pub/week39/html/week39-reveal.html index 620723633..0a0f5455a 100644 --- a/doc/pub/week39/html/week39-reveal.html +++ b/doc/pub/week39/html/week39-reveal.html @@ -184,7 +184,7 @@ MathJax.Hub.Config({
    -

    Sep 29, 2022

    +

    Sep 30, 2022


    @@ -3151,6 +3151,84 @@ plt.show() +
    +

    Same code but now with momentum gradient descent

    + + +
    +
    +
    +
    +
    +
    # Using Autograd to calculate gradients for OLS
    +from random import random, seed
    +import numpy as np
    +import autograd.numpy as np
    +import matplotlib.pyplot as plt
    +from autograd import grad
    +
    +def CostOLS(beta):
    +    return (1.0/n)*np.sum((y-X @ beta)**2)
    +
    +n = 100
    +x = 2*np.random.rand(n,1)
    +y = 4+3*x#+np.random.randn(n,1)
    +
    +X = np.c_[np.ones((n,1)), x]
    +XT_X = X.T @ X
    +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
    +print("Own inversion")
    +print(theta_linreg)
    +# Hessian matrix
    +H = (2.0/n)* XT_X
    +EigValues, EigVectors = np.linalg.eig(H)
    +print(f"Eigenvalues of Hessian Matrix:{EigValues}")
    +
    +theta = np.random.randn(2,1)
    +eta = 1.0/np.max(EigValues)
    +Niterations = 30
    +
    +# define the gradient
    +training_gradient = grad(CostOLS)
    +
    +for iter in range(Niterations):
    +    gradients = training_gradient(theta)
    +    theta -= eta*gradients
    +    print(iter,gradients[0],gradients[1])
    +print("theta from own gd")
    +print(theta)
    +
    +# Now improve with momentum gradient descent
    +change = 0.0
    +delta_momentum = 0.3
    +for iter in range(Niterations):
    +    # calculate gradient
    +    gradients = training_gradient(theta)
    +    # calculate update
    +    new_change = eta*gradients+delta_momentum*change
    +    # take a step
    +    theta -= new_change
    +    # save the change
    +    change = new_change
    +    print(iter,gradients[0],gradients[1])
    +print("theta from own gd wth momentum")
    +print(theta)
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +

    Including Stochastic Gradient Descent with Autograd

    In this code we include the stochastic gradient descent approach discussed above. Note here that we specify which argument we are taking the derivative with respect to when using autograd.

    @@ -3251,6 +3329,98 @@ theta = np.random.randn(2,Same code but now with momentum gradient descent + + +
    +
    +
    +
    +
    +
    # Using Autograd to calculate gradients using SGD
    +# OLS example
    +from random import random, seed
    +import numpy as np
    +import autograd.numpy as np
    +import matplotlib.pyplot as plt
    +from autograd import grad
    +
    +# Note change from previous example
    +def CostOLS(y,X,theta):
    +    return np.sum((y-X @ theta)**2)
    +
    +n = 100
    +x = 2*np.random.rand(n,1)
    +y = 4+3*x+np.random.randn(n,1)
    +
    +X = np.c_[np.ones((n,1)), x]
    +XT_X = X.T @ X
    +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
    +print("Own inversion")
    +print(theta_linreg)
    +# Hessian matrix
    +H = (2.0/n)* XT_X
    +EigValues, EigVectors = np.linalg.eig(H)
    +print(f"Eigenvalues of Hessian Matrix:{EigValues}")
    +
    +theta = np.random.randn(2,1)
    +eta = 1.0/np.max(EigValues)
    +Niterations = 100
    +
    +# Note that we request the derivative wrt third argument (theta, 2 here)
    +training_gradient = grad(CostOLS,2)
    +
    +for iter in range(Niterations):
    +    gradients = (1.0/n)*training_gradient(y, X, theta)
    +    theta -= eta*gradients
    +print("theta from own gd")
    +print(theta)
    +
    +
    +n_epochs = 50
    +M = 5   #size of each minibatch
    +m = int(n/M) #number of minibatches
    +t0, t1 = 5, 50
    +def learning_schedule(t):
    +    return t0/(t+t1)
    +
    +theta = np.random.randn(2,1)
    +
    +change = 0.0
    +delta_momentum = 0.3
    +
    +for epoch in range(n_epochs):
    +    for i in range(m):
    +        random_index = M*np.random.randint(m)
    +        xi = X[random_index:random_index+M]
    +        yi = y[random_index:random_index+M]
    +        gradients = (1.0/M)*training_gradient(yi, xi, theta)
    +        eta = learning_schedule(epoch*m+i)
    +        # calculate update
    +        new_change = eta*gradients+delta_momentum*change
    +        # take a step
    +        theta -= new_change
    +        # save the change
    +        change = new_change
    +print("theta from own sdg with momentum")
    +print(theta)
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +

    And Logistic Regression

    diff --git a/doc/pub/week39/html/week39-solarized.html b/doc/pub/week39/html/week39-solarized.html index 6914a3f5e..0cdc418fa 100644 --- a/doc/pub/week39/html/week39-solarized.html +++ b/doc/pub/week39/html/week39-solarized.html @@ -259,10 +259,18 @@ div.toc p,a { 'the-syntax-a-dot-b-when-finding-the-dot-product'), ('Recommended to avoid', 2, None, 'recommended-to-avoid'), ('Using Autograd with OLS', 2, None, 'using-autograd-with-ols'), + ('Same code but now with momentum gradient descent', + 2, + None, + 'same-code-but-now-with-momentum-gradient-descent'), ('Including Stochastic Gradient Descent with Autograd', 2, None, 'including-stochastic-gradient-descent-with-autograd'), + ('Same code but now with momentum gradient descent', + 2, + None, + 'same-code-but-now-with-momentum-gradient-descent'), ('And Logistic Regression', 2, None, 'and-logistic-regression'), ('Introducing "JAX":"https://jax.readthedocs.io/en/latest/"', 2, @@ -305,7 +313,7 @@ MathJax.Hub.Config({
    -

    Sep 29, 2022

    +

    Sep 30, 2022


    @@ -3061,6 +3069,84 @@ plt.show() +









    +

    Same code but now with momentum gradient descent

    + + +
    +
    +
    +
    +
    +
    # Using Autograd to calculate gradients for OLS
    +from random import random, seed
    +import numpy as np
    +import autograd.numpy as np
    +import matplotlib.pyplot as plt
    +from autograd import grad
    +
    +def CostOLS(beta):
    +    return (1.0/n)*np.sum((y-X @ beta)**2)
    +
    +n = 100
    +x = 2*np.random.rand(n,1)
    +y = 4+3*x#+np.random.randn(n,1)
    +
    +X = np.c_[np.ones((n,1)), x]
    +XT_X = X.T @ X
    +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
    +print("Own inversion")
    +print(theta_linreg)
    +# Hessian matrix
    +H = (2.0/n)* XT_X
    +EigValues, EigVectors = np.linalg.eig(H)
    +print(f"Eigenvalues of Hessian Matrix:{EigValues}")
    +
    +theta = np.random.randn(2,1)
    +eta = 1.0/np.max(EigValues)
    +Niterations = 30
    +
    +# define the gradient
    +training_gradient = grad(CostOLS)
    +
    +for iter in range(Niterations):
    +    gradients = training_gradient(theta)
    +    theta -= eta*gradients
    +    print(iter,gradients[0],gradients[1])
    +print("theta from own gd")
    +print(theta)
    +
    +# Now improve with momentum gradient descent
    +change = 0.0
    +delta_momentum = 0.3
    +for iter in range(Niterations):
    +    # calculate gradient
    +    gradients = training_gradient(theta)
    +    # calculate update
    +    new_change = eta*gradients+delta_momentum*change
    +    # take a step
    +    theta -= new_change
    +    # save the change
    +    change = new_change
    +    print(iter,gradients[0],gradients[1])
    +print("theta from own gd wth momentum")
    +print(theta)
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    + +









    Including Stochastic Gradient Descent with Autograd

    In this code we include the stochastic gradient descent approach discussed above. Note here that we specify which argument we are taking the derivative with respect to when using autograd.

    @@ -3161,6 +3247,98 @@ theta = np.random.randn(2,Same code but now with momentum gradient descent + + +
    +
    +
    +
    +
    +
    # Using Autograd to calculate gradients using SGD
    +# OLS example
    +from random import random, seed
    +import numpy as np
    +import autograd.numpy as np
    +import matplotlib.pyplot as plt
    +from autograd import grad
    +
    +# Note change from previous example
    +def CostOLS(y,X,theta):
    +    return np.sum((y-X @ theta)**2)
    +
    +n = 100
    +x = 2*np.random.rand(n,1)
    +y = 4+3*x+np.random.randn(n,1)
    +
    +X = np.c_[np.ones((n,1)), x]
    +XT_X = X.T @ X
    +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
    +print("Own inversion")
    +print(theta_linreg)
    +# Hessian matrix
    +H = (2.0/n)* XT_X
    +EigValues, EigVectors = np.linalg.eig(H)
    +print(f"Eigenvalues of Hessian Matrix:{EigValues}")
    +
    +theta = np.random.randn(2,1)
    +eta = 1.0/np.max(EigValues)
    +Niterations = 100
    +
    +# Note that we request the derivative wrt third argument (theta, 2 here)
    +training_gradient = grad(CostOLS,2)
    +
    +for iter in range(Niterations):
    +    gradients = (1.0/n)*training_gradient(y, X, theta)
    +    theta -= eta*gradients
    +print("theta from own gd")
    +print(theta)
    +
    +
    +n_epochs = 50
    +M = 5   #size of each minibatch
    +m = int(n/M) #number of minibatches
    +t0, t1 = 5, 50
    +def learning_schedule(t):
    +    return t0/(t+t1)
    +
    +theta = np.random.randn(2,1)
    +
    +change = 0.0
    +delta_momentum = 0.3
    +
    +for epoch in range(n_epochs):
    +    for i in range(m):
    +        random_index = M*np.random.randint(m)
    +        xi = X[random_index:random_index+M]
    +        yi = y[random_index:random_index+M]
    +        gradients = (1.0/M)*training_gradient(yi, xi, theta)
    +        eta = learning_schedule(epoch*m+i)
    +        # calculate update
    +        new_change = eta*gradients+delta_momentum*change
    +        # take a step
    +        theta -= new_change
    +        # save the change
    +        change = new_change
    +print("theta from own sdg with momentum")
    +print(theta)
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    + +









    And Logistic Regression

    diff --git a/doc/pub/week39/html/week39.html b/doc/pub/week39/html/week39.html index 85922cc9a..e2a59f217 100644 --- a/doc/pub/week39/html/week39.html +++ b/doc/pub/week39/html/week39.html @@ -336,10 +336,18 @@ div.toc p,a { 'the-syntax-a-dot-b-when-finding-the-dot-product'), ('Recommended to avoid', 2, None, 'recommended-to-avoid'), ('Using Autograd with OLS', 2, None, 'using-autograd-with-ols'), + ('Same code but now with momentum gradient descent', + 2, + None, + 'same-code-but-now-with-momentum-gradient-descent'), ('Including Stochastic Gradient Descent with Autograd', 2, None, 'including-stochastic-gradient-descent-with-autograd'), + ('Same code but now with momentum gradient descent', + 2, + None, + 'same-code-but-now-with-momentum-gradient-descent'), ('And Logistic Regression', 2, None, 'and-logistic-regression'), ('Introducing "JAX":"https://jax.readthedocs.io/en/latest/"', 2, @@ -382,7 +390,7 @@ MathJax.Hub.Config({
    -

    Sep 29, 2022

    +

    Sep 30, 2022


    @@ -3138,6 +3146,84 @@ plt.show() +









    +

    Same code but now with momentum gradient descent

    + + +
    +
    +
    +
    +
    +
    # Using Autograd to calculate gradients for OLS
    +from random import random, seed
    +import numpy as np
    +import autograd.numpy as np
    +import matplotlib.pyplot as plt
    +from autograd import grad
    +
    +def CostOLS(beta):
    +    return (1.0/n)*np.sum((y-X @ beta)**2)
    +
    +n = 100
    +x = 2*np.random.rand(n,1)
    +y = 4+3*x#+np.random.randn(n,1)
    +
    +X = np.c_[np.ones((n,1)), x]
    +XT_X = X.T @ X
    +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
    +print("Own inversion")
    +print(theta_linreg)
    +# Hessian matrix
    +H = (2.0/n)* XT_X
    +EigValues, EigVectors = np.linalg.eig(H)
    +print(f"Eigenvalues of Hessian Matrix:{EigValues}")
    +
    +theta = np.random.randn(2,1)
    +eta = 1.0/np.max(EigValues)
    +Niterations = 30
    +
    +# define the gradient
    +training_gradient = grad(CostOLS)
    +
    +for iter in range(Niterations):
    +    gradients = training_gradient(theta)
    +    theta -= eta*gradients
    +    print(iter,gradients[0],gradients[1])
    +print("theta from own gd")
    +print(theta)
    +
    +# Now improve with momentum gradient descent
    +change = 0.0
    +delta_momentum = 0.3
    +for iter in range(Niterations):
    +    # calculate gradient
    +    gradients = training_gradient(theta)
    +    # calculate update
    +    new_change = eta*gradients+delta_momentum*change
    +    # take a step
    +    theta -= new_change
    +    # save the change
    +    change = new_change
    +    print(iter,gradients[0],gradients[1])
    +print("theta from own gd wth momentum")
    +print(theta)
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    + +









    Including Stochastic Gradient Descent with Autograd

    In this code we include the stochastic gradient descent approach discussed above. Note here that we specify which argument we are taking the derivative with respect to when using autograd.

    @@ -3238,6 +3324,98 @@ theta = np. +









    +

    Same code but now with momentum gradient descent

    + + +
    +
    +
    +
    +
    +
    # Using Autograd to calculate gradients using SGD
    +# OLS example
    +from random import random, seed
    +import numpy as np
    +import autograd.numpy as np
    +import matplotlib.pyplot as plt
    +from autograd import grad
    +
    +# Note change from previous example
    +def CostOLS(y,X,theta):
    +    return np.sum((y-X @ theta)**2)
    +
    +n = 100
    +x = 2*np.random.rand(n,1)
    +y = 4+3*x+np.random.randn(n,1)
    +
    +X = np.c_[np.ones((n,1)), x]
    +XT_X = X.T @ X
    +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)
    +print("Own inversion")
    +print(theta_linreg)
    +# Hessian matrix
    +H = (2.0/n)* XT_X
    +EigValues, EigVectors = np.linalg.eig(H)
    +print(f"Eigenvalues of Hessian Matrix:{EigValues}")
    +
    +theta = np.random.randn(2,1)
    +eta = 1.0/np.max(EigValues)
    +Niterations = 100
    +
    +# Note that we request the derivative wrt third argument (theta, 2 here)
    +training_gradient = grad(CostOLS,2)
    +
    +for iter in range(Niterations):
    +    gradients = (1.0/n)*training_gradient(y, X, theta)
    +    theta -= eta*gradients
    +print("theta from own gd")
    +print(theta)
    +
    +
    +n_epochs = 50
    +M = 5   #size of each minibatch
    +m = int(n/M) #number of minibatches
    +t0, t1 = 5, 50
    +def learning_schedule(t):
    +    return t0/(t+t1)
    +
    +theta = np.random.randn(2,1)
    +
    +change = 0.0
    +delta_momentum = 0.3
    +
    +for epoch in range(n_epochs):
    +    for i in range(m):
    +        random_index = M*np.random.randint(m)
    +        xi = X[random_index:random_index+M]
    +        yi = y[random_index:random_index+M]
    +        gradients = (1.0/M)*training_gradient(yi, xi, theta)
    +        eta = learning_schedule(epoch*m+i)
    +        # calculate update
    +        new_change = eta*gradients+delta_momentum*change
    +        # take a step
    +        theta -= new_change
    +        # save the change
    +        change = new_change
    +print("theta from own sdg with momentum")
    +print(theta)
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    +
    + +









    And Logistic Regression

    diff --git a/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz b/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz index 75e5c1cb8..3b8729668 100644 Binary files a/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz and b/doc/pub/week39/ipynb/ipynb-week39-src.tar.gz differ diff --git a/doc/pub/week39/ipynb/week39.ipynb b/doc/pub/week39/ipynb/week39.ipynb index be6bd4b83..181035900 100644 --- a/doc/pub/week39/ipynb/week39.ipynb +++ b/doc/pub/week39/ipynb/week39.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "266a989e", + "id": "a3f2a49b", "metadata": { "editable": true }, @@ -14,7 +14,7 @@ }, { "cell_type": "markdown", - "id": "0556a7de", + "id": "ad1e01ff", "metadata": { "editable": true }, @@ -22,14 +22,14 @@ "# Week 39: Optimization and Gradient Methods\n", "**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State University\n", "\n", - "Date: **Sep 29, 2022**\n", + "Date: **Sep 30, 2022**\n", "\n", "Copyright 1999-2022, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license" ] }, { "cell_type": "markdown", - "id": "fba6c7a0", + "id": "1ce0a810", "metadata": { "editable": true }, @@ -55,7 +55,7 @@ }, { "cell_type": "markdown", - "id": "b647f340", + "id": "8ea463cd", "metadata": { "editable": true }, @@ -76,7 +76,7 @@ }, { "cell_type": "markdown", - "id": "579fc60e", + "id": "7bcc3cee", "metadata": { "editable": true }, @@ -93,7 +93,7 @@ }, { "cell_type": "markdown", - "id": "ab60fa1c", + "id": "8258a91e", "metadata": { "editable": true }, @@ -108,7 +108,7 @@ }, { "cell_type": "markdown", - "id": "86e72472", + "id": "fc4d239c", "metadata": { "editable": true }, @@ -118,7 +118,7 @@ }, { "cell_type": "markdown", - "id": "d313a983", + "id": "5e411ed9", "metadata": { "editable": true }, @@ -134,7 +134,7 @@ }, { "cell_type": "markdown", - "id": "a81835f9", + "id": "8559bad7", "metadata": { "editable": true }, @@ -146,7 +146,7 @@ }, { "cell_type": "markdown", - "id": "a6eb84ee", + "id": "8cb89fe1", "metadata": { "editable": true }, @@ -157,7 +157,7 @@ }, { "cell_type": "markdown", - "id": "63610cb8", + "id": "4e7f874e", "metadata": { "editable": true }, @@ -169,7 +169,7 @@ }, { "cell_type": "markdown", - "id": "b1049f68", + "id": "1f6c709e", "metadata": { "editable": true }, @@ -179,7 +179,7 @@ }, { "cell_type": "markdown", - "id": "61782910", + "id": "618e2feb", "metadata": { "editable": true }, @@ -193,7 +193,7 @@ }, { "cell_type": "markdown", - "id": "889f7cf7", + "id": "19b17dae", "metadata": { "editable": true }, @@ -205,7 +205,7 @@ }, { "cell_type": "markdown", - "id": "cf341d16", + "id": "9ad58b9c", "metadata": { "editable": true }, @@ -215,7 +215,7 @@ }, { "cell_type": "markdown", - "id": "99617477", + "id": "6763d42d", "metadata": { "editable": true }, @@ -227,7 +227,7 @@ }, { "cell_type": "markdown", - "id": "048efeba", + "id": "ce95b330", "metadata": { "editable": true }, @@ -239,7 +239,7 @@ }, { "cell_type": "markdown", - "id": "7eeb2a3f", + "id": "4c7a5561", "metadata": { "editable": true }, @@ -259,7 +259,7 @@ }, { "cell_type": "markdown", - "id": "4b243b9f", + "id": "b2f5f204", "metadata": { "editable": true }, @@ -275,7 +275,7 @@ }, { "cell_type": "markdown", - "id": "232e81e8", + "id": "25f1ec55", "metadata": { "editable": true }, @@ -291,7 +291,7 @@ }, { "cell_type": "markdown", - "id": "5048b0ce", + "id": "12800af1", "metadata": { "editable": true }, @@ -302,7 +302,7 @@ }, { "cell_type": "markdown", - "id": "26fc11aa", + "id": "24eac4c2", "metadata": { "editable": true }, @@ -314,7 +314,7 @@ }, { "cell_type": "markdown", - "id": "268b40ad", + "id": "fefe799b", "metadata": { "editable": true }, @@ -324,7 +324,7 @@ }, { "cell_type": "markdown", - "id": "1e8df82c", + "id": "ac788ad2", "metadata": { "editable": true }, @@ -336,7 +336,7 @@ }, { "cell_type": "markdown", - "id": "9dad7b22", + "id": "6837c77b", "metadata": { "editable": true }, @@ -346,7 +346,7 @@ }, { "cell_type": "markdown", - "id": "40fb9be6", + "id": "ea620a33", "metadata": { "editable": true }, @@ -358,7 +358,7 @@ }, { "cell_type": "markdown", - "id": "3d5b5285", + "id": "12439301", "metadata": { "editable": true }, @@ -380,7 +380,7 @@ }, { "cell_type": "markdown", - "id": "92ce5653", + "id": "2991e70c", "metadata": { "editable": true }, @@ -393,7 +393,7 @@ }, { "cell_type": "markdown", - "id": "0defe140", + "id": "0bb42d16", "metadata": { "editable": true }, @@ -406,7 +406,7 @@ }, { "cell_type": "markdown", - "id": "9c81d7b3", + "id": "cc166498", "metadata": { "editable": true }, @@ -416,7 +416,7 @@ }, { "cell_type": "markdown", - "id": "5ec8b48f", + "id": "941aae4f", "metadata": { "editable": true }, @@ -434,7 +434,7 @@ }, { "cell_type": "markdown", - "id": "39d6e38f", + "id": "957408a1", "metadata": { "editable": true }, @@ -444,7 +444,7 @@ }, { "cell_type": "markdown", - "id": "ae2035ba", + "id": "3d6b3d4a", "metadata": { "editable": true }, @@ -459,7 +459,7 @@ }, { "cell_type": "markdown", - "id": "96aca494", + "id": "17d6443c", "metadata": { "editable": true }, @@ -469,7 +469,7 @@ }, { "cell_type": "markdown", - "id": "a8df5a1b", + "id": "583a3f9c", "metadata": { "editable": true }, @@ -483,7 +483,7 @@ }, { "cell_type": "markdown", - "id": "1b245bd4", + "id": "55065669", "metadata": { "editable": true }, @@ -493,7 +493,7 @@ }, { "cell_type": "markdown", - "id": "2e07b2d6", + "id": "8ca8e144", "metadata": { "editable": true }, @@ -507,7 +507,7 @@ }, { "cell_type": "markdown", - "id": "1996803f", + "id": "9173c236", "metadata": { "editable": true }, @@ -522,7 +522,7 @@ }, { "cell_type": "markdown", - "id": "a0a6d7cd", + "id": "34850c9f", "metadata": { "editable": true }, @@ -539,7 +539,7 @@ }, { "cell_type": "markdown", - "id": "be745e5e", + "id": "524b6a67", "metadata": { "editable": true }, @@ -551,7 +551,7 @@ }, { "cell_type": "markdown", - "id": "1fb290b3", + "id": "b1d69a19", "metadata": { "editable": true }, @@ -565,7 +565,7 @@ }, { "cell_type": "markdown", - "id": "5daf5176", + "id": "5c48c377", "metadata": { "editable": true }, @@ -580,7 +580,7 @@ }, { "cell_type": "markdown", - "id": "01556473", + "id": "d847b012", "metadata": { "editable": true }, @@ -592,7 +592,7 @@ }, { "cell_type": "markdown", - "id": "bb911199", + "id": "430c26cb", "metadata": { "editable": true }, @@ -603,7 +603,7 @@ }, { "cell_type": "markdown", - "id": "8c872bd4", + "id": "988b26d1", "metadata": { "editable": true }, @@ -631,7 +631,7 @@ }, { "cell_type": "markdown", - "id": "e89e8886", + "id": "5527bcc6", "metadata": { "editable": true }, @@ -653,7 +653,7 @@ }, { "cell_type": "markdown", - "id": "b975f76b", + "id": "85a3c5fc", "metadata": { "editable": true }, @@ -675,7 +675,7 @@ }, { "cell_type": "markdown", - "id": "5c99441c", + "id": "103ca9e0", "metadata": { "editable": true }, @@ -687,7 +687,7 @@ }, { "cell_type": "markdown", - "id": "04eb49ef", + "id": "5f043b19", "metadata": { "editable": true }, @@ -724,7 +724,7 @@ }, { "cell_type": "markdown", - "id": "5fdf3577", + "id": "7c5ac593", "metadata": { "editable": true }, @@ -752,7 +752,7 @@ }, { "cell_type": "markdown", - "id": "49b00751", + "id": "6bd0a3f3", "metadata": { "editable": true }, @@ -782,7 +782,7 @@ }, { "cell_type": "markdown", - "id": "eda76339", + "id": "3e40ce1f", "metadata": { "editable": true }, @@ -802,7 +802,7 @@ }, { "cell_type": "markdown", - "id": "853dc57b", + "id": "ade236c4", "metadata": { "editable": true }, @@ -814,7 +814,7 @@ }, { "cell_type": "markdown", - "id": "5379fa5a", + "id": "e1c7967b", "metadata": { "editable": true }, @@ -824,7 +824,7 @@ }, { "cell_type": "markdown", - "id": "a49c9a24", + "id": "ad63044a", "metadata": { "editable": true }, @@ -836,7 +836,7 @@ }, { "cell_type": "markdown", - "id": "42e20b26", + "id": "2151a6e5", "metadata": { "editable": true }, @@ -848,7 +848,7 @@ }, { "cell_type": "markdown", - "id": "c0af34f2", + "id": "08d1dcfc", "metadata": { "editable": true }, @@ -860,7 +860,7 @@ }, { "cell_type": "markdown", - "id": "93e49cad", + "id": "714a13d7", "metadata": { "editable": true }, @@ -872,7 +872,7 @@ }, { "cell_type": "markdown", - "id": "991756d0", + "id": "5447991e", "metadata": { "editable": true }, @@ -883,7 +883,7 @@ }, { "cell_type": "markdown", - "id": "9e0bab1f", + "id": "9d8d5636", "metadata": { "editable": true }, @@ -896,7 +896,7 @@ }, { "cell_type": "markdown", - "id": "b7e002d6", + "id": "c4c62bf8", "metadata": { "editable": true }, @@ -908,7 +908,7 @@ }, { "cell_type": "markdown", - "id": "41a88a3a", + "id": "e75ef87c", "metadata": { "editable": true }, @@ -918,7 +918,7 @@ }, { "cell_type": "markdown", - "id": "bf7d650e", + "id": "967df7e5", "metadata": { "editable": true }, @@ -930,7 +930,7 @@ }, { "cell_type": "markdown", - "id": "21b306ec", + "id": "9770216b", "metadata": { "editable": true }, @@ -940,7 +940,7 @@ }, { "cell_type": "markdown", - "id": "c2595c2c", + "id": "600bb9bf", "metadata": { "editable": true }, @@ -951,7 +951,7 @@ }, { "cell_type": "markdown", - "id": "d9631d23", + "id": "1172ae79", "metadata": { "editable": true }, @@ -963,7 +963,7 @@ }, { "cell_type": "markdown", - "id": "e9520154", + "id": "3760450e", "metadata": { "editable": true }, @@ -975,7 +975,7 @@ }, { "cell_type": "markdown", - "id": "fe8ab6ff", + "id": "b0441a59", "metadata": { "editable": true }, @@ -987,7 +987,7 @@ }, { "cell_type": "markdown", - "id": "48b2df0a", + "id": "88c2cdab", "metadata": { "editable": true }, @@ -998,7 +998,7 @@ }, { "cell_type": "markdown", - "id": "3120af28", + "id": "df8ddc4b", "metadata": { "editable": true }, @@ -1009,7 +1009,7 @@ }, { "cell_type": "markdown", - "id": "349f81a7", + "id": "0f297c3d", "metadata": { "editable": true }, @@ -1021,7 +1021,7 @@ }, { "cell_type": "markdown", - "id": "2fa0028e", + "id": "0d0d0aa2", "metadata": { "editable": true }, @@ -1031,7 +1031,7 @@ }, { "cell_type": "markdown", - "id": "f9ec53a2", + "id": "c4d38149", "metadata": { "editable": true }, @@ -1043,7 +1043,7 @@ }, { "cell_type": "markdown", - "id": "e881ad44", + "id": "d4372560", "metadata": { "editable": true }, @@ -1053,7 +1053,7 @@ }, { "cell_type": "markdown", - "id": "9e652153", + "id": "e718ad83", "metadata": { "editable": true }, @@ -1065,7 +1065,7 @@ }, { "cell_type": "markdown", - "id": "235b588a", + "id": "457ea769", "metadata": { "editable": true }, @@ -1075,7 +1075,7 @@ }, { "cell_type": "markdown", - "id": "1b1cf967", + "id": "e5677842", "metadata": { "editable": true }, @@ -1087,7 +1087,7 @@ }, { "cell_type": "markdown", - "id": "998bb78c", + "id": "f41dcf2f", "metadata": { "editable": true }, @@ -1097,7 +1097,7 @@ }, { "cell_type": "markdown", - "id": "59cccb52", + "id": "c3a84d60", "metadata": { "editable": true }, @@ -1109,7 +1109,7 @@ }, { "cell_type": "markdown", - "id": "8e031fcc", + "id": "54eec9d8", "metadata": { "editable": true }, @@ -1120,7 +1120,7 @@ { "cell_type": "code", "execution_count": 1, - "id": "42d5c457", + "id": "e65d59f3", "metadata": { "collapsed": false, "editable": true @@ -1153,7 +1153,7 @@ }, { "cell_type": "markdown", - "id": "3d0317f0", + "id": "a1b1ec51", "metadata": { "editable": true }, @@ -1164,7 +1164,7 @@ { "cell_type": "code", "execution_count": 2, - "id": "23fce077", + "id": "d509c489", "metadata": { "collapsed": false, "editable": true @@ -1178,7 +1178,7 @@ }, { "cell_type": "markdown", - "id": "19c01aa8", + "id": "0a234ac1", "metadata": { "editable": true }, @@ -1189,7 +1189,7 @@ { "cell_type": "code", "execution_count": 3, - "id": "3c018860", + "id": "6dee8ee2", "metadata": { "collapsed": false, "editable": true @@ -1202,7 +1202,7 @@ }, { "cell_type": "markdown", - "id": "07d57e1b", + "id": "6739f9d4", "metadata": { "editable": true }, @@ -1213,7 +1213,7 @@ { "cell_type": "code", "execution_count": 4, - "id": "5f77c7ff", + "id": "a635a93f", "metadata": { "collapsed": false, "editable": true @@ -1231,7 +1231,7 @@ }, { "cell_type": "markdown", - "id": "31dfadfb", + "id": "1738466f", "metadata": { "editable": true }, @@ -1242,7 +1242,7 @@ { "cell_type": "code", "execution_count": 5, - "id": "0cedd6a4", + "id": "c58d0115", "metadata": { "collapsed": false, "editable": true @@ -1257,7 +1257,7 @@ }, { "cell_type": "markdown", - "id": "c9a57db7", + "id": "7ea54653", "metadata": { "editable": true }, @@ -1267,7 +1267,7 @@ }, { "cell_type": "markdown", - "id": "720e9388", + "id": "dd7e00ed", "metadata": { "editable": true }, @@ -1281,7 +1281,7 @@ }, { "cell_type": "markdown", - "id": "ccdd73b8", + "id": "7f01fc09", "metadata": { "editable": true }, @@ -1293,7 +1293,7 @@ }, { "cell_type": "markdown", - "id": "76cd5686", + "id": "6205a0c0", "metadata": { "editable": true }, @@ -1304,7 +1304,7 @@ }, { "cell_type": "markdown", - "id": "001483f6", + "id": "0a474db3", "metadata": { "editable": true }, @@ -1316,7 +1316,7 @@ }, { "cell_type": "markdown", - "id": "29da6d30", + "id": "281c3e87", "metadata": { "editable": true }, @@ -1327,7 +1327,7 @@ }, { "cell_type": "markdown", - "id": "1cccf0ea", + "id": "72adadf6", "metadata": { "editable": true }, @@ -1338,7 +1338,7 @@ }, { "cell_type": "markdown", - "id": "ef9f4a74", + "id": "b1a60997", "metadata": { "editable": true }, @@ -1350,7 +1350,7 @@ }, { "cell_type": "markdown", - "id": "cf59c588", + "id": "4d655399", "metadata": { "editable": true }, @@ -1360,7 +1360,7 @@ }, { "cell_type": "markdown", - "id": "a2106bb7", + "id": "116ebbac", "metadata": { "editable": true }, @@ -1372,7 +1372,7 @@ }, { "cell_type": "markdown", - "id": "f3cd45cd", + "id": "55185d8f", "metadata": { "editable": true }, @@ -1384,7 +1384,7 @@ }, { "cell_type": "markdown", - "id": "ae979921", + "id": "afd187d3", "metadata": { "editable": true }, @@ -1396,7 +1396,7 @@ }, { "cell_type": "markdown", - "id": "ad4149ba", + "id": "06b67a8c", "metadata": { "editable": true }, @@ -1408,7 +1408,7 @@ }, { "cell_type": "markdown", - "id": "9d4a2bb5", + "id": "a6c814bc", "metadata": { "editable": true }, @@ -1419,7 +1419,7 @@ }, { "cell_type": "markdown", - "id": "876225f8", + "id": "8e4c87d2", "metadata": { "editable": true }, @@ -1431,7 +1431,7 @@ }, { "cell_type": "markdown", - "id": "fd69b613", + "id": "71ac6a87", "metadata": { "editable": true }, @@ -1441,7 +1441,7 @@ }, { "cell_type": "markdown", - "id": "f2de9fd2", + "id": "2a5a1639", "metadata": { "editable": true }, @@ -1453,7 +1453,7 @@ }, { "cell_type": "markdown", - "id": "31e2d00e", + "id": "a21020d2", "metadata": { "editable": true }, @@ -1463,7 +1463,7 @@ }, { "cell_type": "markdown", - "id": "88c8b3ca", + "id": "fc5b548c", "metadata": { "editable": true }, @@ -1475,7 +1475,7 @@ }, { "cell_type": "markdown", - "id": "af1050c2", + "id": "52b40bde", "metadata": { "editable": true }, @@ -1495,7 +1495,7 @@ }, { "cell_type": "markdown", - "id": "8f3d0c5e", + "id": "361660ca", "metadata": { "editable": true }, @@ -1507,7 +1507,7 @@ }, { "cell_type": "markdown", - "id": "e8d03d21", + "id": "519d134c", "metadata": { "editable": true }, @@ -1517,7 +1517,7 @@ }, { "cell_type": "markdown", - "id": "0c4d6c4b", + "id": "0a7faed4", "metadata": { "editable": true }, @@ -1529,7 +1529,7 @@ }, { "cell_type": "markdown", - "id": "22c9d676", + "id": "e7a1c68e", "metadata": { "editable": true }, @@ -1539,7 +1539,7 @@ }, { "cell_type": "markdown", - "id": "33f135cc", + "id": "6226cc6a", "metadata": { "editable": true }, @@ -1550,7 +1550,7 @@ }, { "cell_type": "markdown", - "id": "4638cec5", + "id": "2e7286a9", "metadata": { "editable": true }, @@ -1562,7 +1562,7 @@ }, { "cell_type": "markdown", - "id": "78176236", + "id": "356d6972", "metadata": { "editable": true }, @@ -1574,7 +1574,7 @@ }, { "cell_type": "markdown", - "id": "c71e2985", + "id": "b8d473e6", "metadata": { "editable": true }, @@ -1586,7 +1586,7 @@ }, { "cell_type": "markdown", - "id": "a82f2c13", + "id": "2d22f619", "metadata": { "editable": true }, @@ -1599,7 +1599,7 @@ }, { "cell_type": "markdown", - "id": "f1726a68", + "id": "16102292", "metadata": { "editable": true }, @@ -1610,7 +1610,7 @@ }, { "cell_type": "markdown", - "id": "d10d83de", + "id": "8cac320f", "metadata": { "editable": true }, @@ -1622,7 +1622,7 @@ }, { "cell_type": "markdown", - "id": "0785af64", + "id": "7d17a555", "metadata": { "editable": true }, @@ -1638,7 +1638,7 @@ }, { "cell_type": "markdown", - "id": "bcbabbe5", + "id": "dfc11703", "metadata": { "editable": true }, @@ -1650,7 +1650,7 @@ }, { "cell_type": "markdown", - "id": "e177513d", + "id": "511d4f6c", "metadata": { "editable": true }, @@ -1661,7 +1661,7 @@ }, { "cell_type": "markdown", - "id": "f62ef65f", + "id": "167ad0dd", "metadata": { "editable": true }, @@ -1673,7 +1673,7 @@ }, { "cell_type": "markdown", - "id": "51c56e08", + "id": "5e3e8f17", "metadata": { "editable": true }, @@ -1683,7 +1683,7 @@ }, { "cell_type": "markdown", - "id": "4b4b7571", + "id": "d17a807f", "metadata": { "editable": true }, @@ -1695,7 +1695,7 @@ }, { "cell_type": "markdown", - "id": "f1d91da4", + "id": "c6545761", "metadata": { "editable": true }, @@ -1705,7 +1705,7 @@ }, { "cell_type": "markdown", - "id": "ca59c2b0", + "id": "48fa9c57", "metadata": { "editable": true }, @@ -1717,7 +1717,7 @@ }, { "cell_type": "markdown", - "id": "35b934d5", + "id": "14ca3319", "metadata": { "editable": true }, @@ -1727,7 +1727,7 @@ }, { "cell_type": "markdown", - "id": "7b81ff02", + "id": "96480b6a", "metadata": { "editable": true }, @@ -1739,7 +1739,7 @@ }, { "cell_type": "markdown", - "id": "cc532b40", + "id": "f5f3cbcf", "metadata": { "editable": true }, @@ -1763,7 +1763,7 @@ { "cell_type": "code", "execution_count": 6, - "id": "fd40460a", + "id": "58aecda1", "metadata": { "collapsed": false, "editable": true @@ -1776,7 +1776,7 @@ }, { "cell_type": "markdown", - "id": "e82a0302", + "id": "307bbb57", "metadata": { "editable": true }, @@ -1787,7 +1787,7 @@ }, { "cell_type": "markdown", - "id": "68a4daa0", + "id": "7a110c79", "metadata": { "editable": true }, @@ -1799,7 +1799,7 @@ }, { "cell_type": "markdown", - "id": "1e3644b8", + "id": "28b08d11", "metadata": { "editable": true }, @@ -1809,7 +1809,7 @@ }, { "cell_type": "markdown", - "id": "b52d0e59", + "id": "c28cd1d6", "metadata": { "editable": true }, @@ -1821,7 +1821,7 @@ }, { "cell_type": "markdown", - "id": "78ad760b", + "id": "4e123dce", "metadata": { "editable": true }, @@ -1835,7 +1835,7 @@ }, { "cell_type": "markdown", - "id": "d71f7c0d", + "id": "d4dae5a6", "metadata": { "editable": true }, @@ -1851,7 +1851,7 @@ }, { "cell_type": "markdown", - "id": "ce9d921f", + "id": "f02bdc7a", "metadata": { "editable": true }, @@ -1861,7 +1861,7 @@ }, { "cell_type": "markdown", - "id": "27aa96d7", + "id": "8a27caf3", "metadata": { "editable": true }, @@ -1873,7 +1873,7 @@ }, { "cell_type": "markdown", - "id": "9a8d11f0", + "id": "750e00b1", "metadata": { "editable": true }, @@ -1883,7 +1883,7 @@ }, { "cell_type": "markdown", - "id": "f3e390cf", + "id": "24b9613d", "metadata": { "editable": true }, @@ -1895,7 +1895,7 @@ }, { "cell_type": "markdown", - "id": "b7228342", + "id": "68208c7e", "metadata": { "editable": true }, @@ -1909,7 +1909,7 @@ }, { "cell_type": "markdown", - "id": "62c62291", + "id": "456d3441", "metadata": { "editable": true }, @@ -1919,7 +1919,7 @@ }, { "cell_type": "markdown", - "id": "1564d6a0", + "id": "5e9d90e1", "metadata": { "editable": true }, @@ -1930,7 +1930,7 @@ }, { "cell_type": "markdown", - "id": "8cc61784", + "id": "99072c22", "metadata": { "editable": true }, @@ -1945,7 +1945,7 @@ }, { "cell_type": "markdown", - "id": "54ba3f55", + "id": "339105e4", "metadata": { "editable": true }, @@ -1955,7 +1955,7 @@ }, { "cell_type": "markdown", - "id": "c150ee31", + "id": "fa3aeb52", "metadata": { "editable": true }, @@ -1967,7 +1967,7 @@ }, { "cell_type": "markdown", - "id": "5bec697f", + "id": "868d976a", "metadata": { "editable": true }, @@ -1979,7 +1979,7 @@ }, { "cell_type": "markdown", - "id": "639805a1", + "id": "34207df7", "metadata": { "editable": true }, @@ -1994,7 +1994,7 @@ }, { "cell_type": "markdown", - "id": "734657b8", + "id": "38a9d087", "metadata": { "editable": true }, @@ -2007,7 +2007,7 @@ { "cell_type": "code", "execution_count": 7, - "id": "0aec41a8", + "id": "dc7badb6", "metadata": { "collapsed": false, "editable": true @@ -2064,7 +2064,7 @@ }, { "cell_type": "markdown", - "id": "a6d8295f", + "id": "59d83633", "metadata": { "editable": true }, @@ -2075,7 +2075,7 @@ { "cell_type": "code", "execution_count": 8, - "id": "a903b53b", + "id": "ccaa3e15", "metadata": { "collapsed": false, "editable": true @@ -2102,7 +2102,7 @@ }, { "cell_type": "markdown", - "id": "c6eb1b40", + "id": "b6683bd2", "metadata": { "editable": true }, @@ -2114,7 +2114,7 @@ }, { "cell_type": "markdown", - "id": "8a8cd057", + "id": "2bb6c0fe", "metadata": { "editable": true }, @@ -2126,7 +2126,7 @@ }, { "cell_type": "markdown", - "id": "f2983a73", + "id": "e9d915fa", "metadata": { "editable": true }, @@ -2136,7 +2136,7 @@ }, { "cell_type": "markdown", - "id": "e352badb", + "id": "f05f103f", "metadata": { "editable": true }, @@ -2150,7 +2150,7 @@ }, { "cell_type": "markdown", - "id": "aa884fde", + "id": "f1006d59", "metadata": { "editable": true }, @@ -2160,7 +2160,7 @@ }, { "cell_type": "markdown", - "id": "f91ba4e0", + "id": "64482546", "metadata": { "editable": true }, @@ -2172,7 +2172,7 @@ }, { "cell_type": "markdown", - "id": "9a18134e", + "id": "1080bf06", "metadata": { "editable": true }, @@ -2183,7 +2183,7 @@ }, { "cell_type": "markdown", - "id": "d3ec0891", + "id": "030ff44d", "metadata": { "editable": true }, @@ -2198,7 +2198,7 @@ }, { "cell_type": "markdown", - "id": "8a042993", + "id": "69f5722e", "metadata": { "editable": true }, @@ -2212,7 +2212,7 @@ }, { "cell_type": "markdown", - "id": "148f61dd", + "id": "1c05f1a1", "metadata": { "editable": true }, @@ -2223,7 +2223,7 @@ { "cell_type": "code", "execution_count": 9, - "id": "707deca6", + "id": "7f2203f8", "metadata": { "collapsed": false, "editable": true @@ -2284,7 +2284,7 @@ }, { "cell_type": "markdown", - "id": "f39bd794", + "id": "e20918f8", "metadata": { "editable": true }, @@ -2306,7 +2306,7 @@ }, { "cell_type": "markdown", - "id": "41eae562", + "id": "c4cd2dc7", "metadata": { "editable": true }, @@ -2319,7 +2319,7 @@ { "cell_type": "code", "execution_count": 10, - "id": "2f65ff78", + "id": "4d62141d", "metadata": { "collapsed": false, "editable": true @@ -2385,7 +2385,7 @@ }, { "cell_type": "markdown", - "id": "185959cb", + "id": "f6af1edb", "metadata": { "editable": true }, @@ -2396,7 +2396,7 @@ { "cell_type": "code", "execution_count": 11, - "id": "84299e90", + "id": "69121b13", "metadata": { "collapsed": false, "editable": true @@ -2470,7 +2470,7 @@ }, { "cell_type": "markdown", - "id": "cce1ff04", + "id": "f7ace6c9", "metadata": { "editable": true }, @@ -2482,7 +2482,7 @@ }, { "cell_type": "markdown", - "id": "ec171de6", + "id": "9ec31a1d", "metadata": { "editable": true }, @@ -2503,7 +2503,7 @@ }, { "cell_type": "markdown", - "id": "61f771e2", + "id": "397dc9e2", "metadata": { "editable": true }, @@ -2535,7 +2535,7 @@ }, { "cell_type": "markdown", - "id": "167c0617", + "id": "9e34ba26", "metadata": { "editable": true }, @@ -2552,7 +2552,7 @@ }, { "cell_type": "markdown", - "id": "90e37c2b", + "id": "3192866d", "metadata": { "editable": true }, @@ -2565,7 +2565,7 @@ }, { "cell_type": "markdown", - "id": "514288e1", + "id": "d0eacb6e", "metadata": { "editable": true }, @@ -2578,7 +2578,7 @@ }, { "cell_type": "markdown", - "id": "bdb9af89", + "id": "7f71e422", "metadata": { "editable": true }, @@ -2591,7 +2591,7 @@ }, { "cell_type": "markdown", - "id": "3d7c4c2f", + "id": "68fe139a", "metadata": { "editable": true }, @@ -2605,7 +2605,7 @@ }, { "cell_type": "markdown", - "id": "8a5b43e2", + "id": "1b27dbae", "metadata": { "editable": true }, @@ -2627,7 +2627,7 @@ }, { "cell_type": "markdown", - "id": "0d4cc828", + "id": "d3aabc29", "metadata": { "editable": true }, @@ -2642,7 +2642,7 @@ }, { "cell_type": "markdown", - "id": "00ecaf91", + "id": "a393e4dc", "metadata": { "editable": true }, @@ -2654,7 +2654,7 @@ }, { "cell_type": "markdown", - "id": "f0e3f62a", + "id": "454663c9", "metadata": { "editable": true }, @@ -2667,7 +2667,7 @@ }, { "cell_type": "markdown", - "id": "45411341", + "id": "2a2a9b25", "metadata": { "editable": true }, @@ -2681,7 +2681,7 @@ }, { "cell_type": "markdown", - "id": "58c3927e", + "id": "dd7bc0c2", "metadata": { "editable": true }, @@ -2692,7 +2692,7 @@ { "cell_type": "code", "execution_count": 12, - "id": "eb9b91fb", + "id": "7fc19ed3", "metadata": { "collapsed": false, "editable": true @@ -2717,7 +2717,7 @@ }, { "cell_type": "markdown", - "id": "9eba2696", + "id": "b33dfafd", "metadata": { "editable": true }, @@ -2733,7 +2733,7 @@ }, { "cell_type": "markdown", - "id": "9b74dee2", + "id": "224cb6c6", "metadata": { "editable": true }, @@ -2754,7 +2754,7 @@ }, { "cell_type": "markdown", - "id": "1b0f3bea", + "id": "eea96e98", "metadata": { "editable": true }, @@ -2774,7 +2774,7 @@ }, { "cell_type": "markdown", - "id": "d21909a4", + "id": "91dcd74e", "metadata": { "editable": true }, @@ -2793,7 +2793,7 @@ { "cell_type": "code", "execution_count": 13, - "id": "c48e042b", + "id": "0de8404e", "metadata": { "collapsed": false, "editable": true @@ -2828,7 +2828,7 @@ }, { "cell_type": "markdown", - "id": "c3e56e42", + "id": "17e6deb7", "metadata": { "editable": true }, @@ -2841,7 +2841,7 @@ { "cell_type": "code", "execution_count": 14, - "id": "27c60383", + "id": "21917f1d", "metadata": { "collapsed": false, "editable": true @@ -2918,7 +2918,7 @@ }, { "cell_type": "markdown", - "id": "8eb3e862", + "id": "18e21a9b", "metadata": { "editable": true }, @@ -2933,7 +2933,7 @@ }, { "cell_type": "markdown", - "id": "34c63f75", + "id": "90ea9013", "metadata": { "editable": true }, @@ -2948,7 +2948,7 @@ }, { "cell_type": "markdown", - "id": "b526eb97", + "id": "f9fa1bee", "metadata": { "editable": true }, @@ -2960,7 +2960,7 @@ }, { "cell_type": "markdown", - "id": "a65b9f7c", + "id": "66018d26", "metadata": { "editable": true }, @@ -2978,7 +2978,7 @@ }, { "cell_type": "markdown", - "id": "501f2ee6", + "id": "4e034652", "metadata": { "editable": true }, @@ -2997,7 +2997,7 @@ }, { "cell_type": "markdown", - "id": "d7f9d1c1", + "id": "41863ff4", "metadata": { "editable": true }, @@ -3009,7 +3009,7 @@ }, { "cell_type": "markdown", - "id": "2f6d7af0", + "id": "dd86f43f", "metadata": { "editable": true }, @@ -3019,7 +3019,7 @@ }, { "cell_type": "markdown", - "id": "894ce2fd", + "id": "be1bec74", "metadata": { "editable": true }, @@ -3035,7 +3035,7 @@ }, { "cell_type": "markdown", - "id": "11ddd872", + "id": "f67f13cf", "metadata": { "editable": true }, @@ -3047,7 +3047,7 @@ }, { "cell_type": "markdown", - "id": "c0cadfff", + "id": "b00e4d9b", "metadata": { "editable": true }, @@ -3057,7 +3057,7 @@ }, { "cell_type": "markdown", - "id": "428438f5", + "id": "645ea43b", "metadata": { "editable": true }, @@ -3069,7 +3069,7 @@ }, { "cell_type": "markdown", - "id": "92e9357c", + "id": "f3155674", "metadata": { "editable": true }, @@ -3079,7 +3079,7 @@ }, { "cell_type": "markdown", - "id": "5c0e73c6", + "id": "38c39b06", "metadata": { "editable": true }, @@ -3091,7 +3091,7 @@ }, { "cell_type": "markdown", - "id": "0ea1cb78", + "id": "d09fe7eb", "metadata": { "editable": true }, @@ -3107,7 +3107,7 @@ }, { "cell_type": "markdown", - "id": "89f597f5", + "id": "66fd53f5", "metadata": { "editable": true }, @@ -3119,7 +3119,7 @@ }, { "cell_type": "markdown", - "id": "954d762d", + "id": "1b552b2f", "metadata": { "editable": true }, @@ -3152,7 +3152,7 @@ }, { "cell_type": "markdown", - "id": "38641257", + "id": "bec97cf5", "metadata": { "editable": true }, @@ -3164,7 +3164,7 @@ }, { "cell_type": "markdown", - "id": "22c5b321", + "id": "d482e2fe", "metadata": { "editable": true }, @@ -3182,7 +3182,7 @@ }, { "cell_type": "markdown", - "id": "c7309191", + "id": "d0b6d3dc", "metadata": { "editable": true }, @@ -3192,7 +3192,7 @@ }, { "cell_type": "markdown", - "id": "56a16b68", + "id": "0510a924", "metadata": { "editable": true }, @@ -3223,7 +3223,7 @@ }, { "cell_type": "markdown", - "id": "b7c6d361", + "id": "05cebe24", "metadata": { "editable": true }, @@ -3238,7 +3238,7 @@ }, { "cell_type": "markdown", - "id": "78d79a03", + "id": "f6feb0dd", "metadata": { "editable": true }, @@ -3256,7 +3256,7 @@ }, { "cell_type": "markdown", - "id": "ad384bea", + "id": "5c3d2f41", "metadata": { "editable": true }, @@ -3268,7 +3268,7 @@ }, { "cell_type": "markdown", - "id": "a2b9882c", + "id": "718af0d0", "metadata": { "editable": true }, @@ -3280,7 +3280,7 @@ }, { "cell_type": "markdown", - "id": "70d9b90b", + "id": "34265c75", "metadata": { "editable": true }, @@ -3298,7 +3298,7 @@ }, { "cell_type": "markdown", - "id": "6b4c7e03", + "id": "20ed6d82", "metadata": { "editable": true }, @@ -3321,7 +3321,7 @@ }, { "cell_type": "markdown", - "id": "4fc5c09c", + "id": "89caf637", "metadata": { "editable": true }, @@ -3339,7 +3339,7 @@ }, { "cell_type": "markdown", - "id": "6dac702f", + "id": "f5e7a1ac", "metadata": { "editable": true }, @@ -3351,7 +3351,7 @@ }, { "cell_type": "markdown", - "id": "38fb872a", + "id": "5337c5ea", "metadata": { "editable": true }, @@ -3363,7 +3363,7 @@ }, { "cell_type": "markdown", - "id": "48436706", + "id": "9fa7d7c4", "metadata": { "editable": true }, @@ -3375,7 +3375,7 @@ }, { "cell_type": "markdown", - "id": "711b05c2", + "id": "475f8da0", "metadata": { "editable": true }, @@ -3387,7 +3387,7 @@ }, { "cell_type": "markdown", - "id": "7c5be7c4", + "id": "8b337b55", "metadata": { "editable": true }, @@ -3399,7 +3399,7 @@ }, { "cell_type": "markdown", - "id": "91fbb6a2", + "id": "71e429cc", "metadata": { "editable": true }, @@ -3416,7 +3416,7 @@ }, { "cell_type": "markdown", - "id": "e25b8b4b", + "id": "391659f3", "metadata": { "editable": true }, @@ -3435,7 +3435,7 @@ }, { "cell_type": "markdown", - "id": "bc510d5a", + "id": "e1d8467c", "metadata": { "editable": true }, @@ -3447,7 +3447,7 @@ }, { "cell_type": "markdown", - "id": "0957ade8", + "id": "75e2c931", "metadata": { "editable": true }, @@ -3467,7 +3467,7 @@ }, { "cell_type": "markdown", - "id": "c0c4e70c", + "id": "09303519", "metadata": { "editable": true }, @@ -3505,7 +3505,7 @@ }, { "cell_type": "markdown", - "id": "2479e6ec", + "id": "f1be6f6b", "metadata": { "editable": true }, @@ -3517,7 +3517,7 @@ }, { "cell_type": "markdown", - "id": "25573059", + "id": "992f4a02", "metadata": { "editable": true }, @@ -3527,7 +3527,7 @@ }, { "cell_type": "markdown", - "id": "99867c27", + "id": "a4b6d08e", "metadata": { "editable": true }, @@ -3539,7 +3539,7 @@ }, { "cell_type": "markdown", - "id": "ab095717", + "id": "c7d831ef", "metadata": { "editable": true }, @@ -3550,7 +3550,7 @@ { "cell_type": "code", "execution_count": 15, - "id": "957de0bd", + "id": "54acdf56", "metadata": { "collapsed": false, "editable": true @@ -3595,7 +3595,7 @@ }, { "cell_type": "markdown", - "id": "6b45687d", + "id": "9a753f05", "metadata": { "editable": true }, @@ -3612,7 +3612,7 @@ { "cell_type": "code", "execution_count": 16, - "id": "3250aae4", + "id": "b8d502a7", "metadata": { "collapsed": false, "editable": true @@ -3640,7 +3640,7 @@ }, { "cell_type": "markdown", - "id": "4bb0ad31", + "id": "db184594", "metadata": { "editable": true }, @@ -3655,7 +3655,7 @@ { "cell_type": "code", "execution_count": 17, - "id": "515ff31b", + "id": "61613ed6", "metadata": { "collapsed": false, "editable": true @@ -3699,7 +3699,7 @@ }, { "cell_type": "markdown", - "id": "774f5eda", + "id": "7543cea7", "metadata": { "editable": true }, @@ -3709,7 +3709,7 @@ }, { "cell_type": "markdown", - "id": "68b8b8f7", + "id": "a6c094d3", "metadata": { "editable": true }, @@ -3720,7 +3720,7 @@ { "cell_type": "code", "execution_count": 18, - "id": "bcf539e8", + "id": "e57abf0c", "metadata": { "collapsed": false, "editable": true @@ -3748,7 +3748,7 @@ }, { "cell_type": "markdown", - "id": "cb8ac207", + "id": "f63bb9f5", "metadata": { "editable": true }, @@ -3763,7 +3763,7 @@ }, { "cell_type": "markdown", - "id": "c5f87a92", + "id": "11623240", "metadata": { "editable": true }, @@ -3774,7 +3774,7 @@ { "cell_type": "code", "execution_count": 19, - "id": "7d46685e", + "id": "148dce0d", "metadata": { "collapsed": false, "editable": true @@ -3802,7 +3802,7 @@ }, { "cell_type": "markdown", - "id": "1f7f96c0", + "id": "7890ae18", "metadata": { "editable": true }, @@ -3813,7 +3813,7 @@ { "cell_type": "code", "execution_count": 20, - "id": "e49128c7", + "id": "1ec6199c", "metadata": { "collapsed": false, "editable": true @@ -3838,7 +3838,7 @@ }, { "cell_type": "markdown", - "id": "03135da4", + "id": "e87650cc", "metadata": { "editable": true }, @@ -3849,7 +3849,7 @@ { "cell_type": "code", "execution_count": 21, - "id": "968ba627", + "id": "57766f6a", "metadata": { "collapsed": false, "editable": true @@ -3885,7 +3885,7 @@ { "cell_type": "code", "execution_count": 22, - "id": "ebaba4fb", + "id": "43f2b332", "metadata": { "collapsed": false, "editable": true @@ -3905,7 +3905,7 @@ }, { "cell_type": "markdown", - "id": "680f026e", + "id": "b3defe0d", "metadata": { "editable": true }, @@ -3916,7 +3916,7 @@ { "cell_type": "code", "execution_count": 23, - "id": "092801c3", + "id": "666d028d", "metadata": { "collapsed": false, "editable": true @@ -3954,7 +3954,7 @@ }, { "cell_type": "markdown", - "id": "aaaf1ad3", + "id": "1e69976c", "metadata": { "editable": true }, @@ -3964,7 +3964,7 @@ }, { "cell_type": "markdown", - "id": "f5c6d17b", + "id": "f8ca253c", "metadata": { "editable": true }, @@ -3978,7 +3978,7 @@ { "cell_type": "code", "execution_count": 24, - "id": "1d32abfa", + "id": "6f717821", "metadata": { "collapsed": false, "editable": true @@ -4000,7 +4000,7 @@ }, { "cell_type": "markdown", - "id": "5f858e86", + "id": "06c934c9", "metadata": { "editable": true }, @@ -4010,7 +4010,7 @@ }, { "cell_type": "markdown", - "id": "1179e3dd", + "id": "d199bb84", "metadata": { "editable": true }, @@ -4021,7 +4021,7 @@ { "cell_type": "code", "execution_count": 25, - "id": "dba3ad31", + "id": "49073461", "metadata": { "collapsed": false, "editable": true @@ -4043,7 +4043,7 @@ }, { "cell_type": "markdown", - "id": "2d71a680", + "id": "673157aa", "metadata": { "editable": true }, @@ -4056,7 +4056,7 @@ { "cell_type": "code", "execution_count": 26, - "id": "6ac01a09", + "id": "3aecc3d4", "metadata": { "collapsed": false, "editable": true @@ -4081,7 +4081,7 @@ }, { "cell_type": "markdown", - "id": "5becf640", + "id": "5eab280c", "metadata": { "editable": true }, @@ -4093,7 +4093,7 @@ { "cell_type": "code", "execution_count": 27, - "id": "6511ed46", + "id": "ddff4eb6", "metadata": { "collapsed": false, "editable": true @@ -4108,7 +4108,7 @@ }, { "cell_type": "markdown", - "id": "3f49669a", + "id": "ebfad0f4", "metadata": { "editable": true }, @@ -4123,7 +4123,7 @@ { "cell_type": "code", "execution_count": 28, - "id": "7207a933", + "id": "d99a8175", "metadata": { "collapsed": false, "editable": true @@ -4183,7 +4183,82 @@ }, { "cell_type": "markdown", - "id": "eb71cd5b", + "id": "36efc0e8", + "metadata": { + "editable": true + }, + "source": [ + "## Same code but now with momentum gradient descent" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "7f807aeb", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "# Using Autograd to calculate gradients for OLS\n", + "from random import random, seed\n", + "import numpy as np\n", + "import autograd.numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from autograd import grad\n", + "\n", + "def CostOLS(beta):\n", + " return (1.0/n)*np.sum((y-X @ beta)**2)\n", + "\n", + "n = 100\n", + "x = 2*np.random.rand(n,1)\n", + "y = 4+3*x#+np.random.randn(n,1)\n", + "\n", + "X = np.c_[np.ones((n,1)), x]\n", + "XT_X = X.T @ X\n", + "theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)\n", + "print(\"Own inversion\")\n", + "print(theta_linreg)\n", + "# Hessian matrix\n", + "H = (2.0/n)* XT_X\n", + "EigValues, EigVectors = np.linalg.eig(H)\n", + "print(f\"Eigenvalues of Hessian Matrix:{EigValues}\")\n", + "\n", + "theta = np.random.randn(2,1)\n", + "eta = 1.0/np.max(EigValues)\n", + "Niterations = 30\n", + "\n", + "# define the gradient\n", + "training_gradient = grad(CostOLS)\n", + "\n", + "for iter in range(Niterations):\n", + " gradients = training_gradient(theta)\n", + " theta -= eta*gradients\n", + " print(iter,gradients[0],gradients[1])\n", + "print(\"theta from own gd\")\n", + "print(theta)\n", + "\n", + "# Now improve with momentum gradient descent\n", + "change = 0.0\n", + "delta_momentum = 0.3\n", + "for iter in range(Niterations):\n", + " # calculate gradient\n", + " gradients = training_gradient(theta)\n", + " # calculate update\n", + " new_change = eta*gradients+delta_momentum*change\n", + " # take a step\n", + " theta -= new_change\n", + " # save the change\n", + " change = new_change\n", + " print(iter,gradients[0],gradients[1])\n", + "print(\"theta from own gd wth momentum\")\n", + "print(theta)" + ] + }, + { + "cell_type": "markdown", + "id": "02751092", "metadata": { "editable": true }, @@ -4194,8 +4269,8 @@ }, { "cell_type": "code", - "execution_count": 29, - "id": "f5fade01", + "execution_count": 30, + "id": "c958f749", "metadata": { "collapsed": false, "editable": true @@ -4279,7 +4354,96 @@ }, { "cell_type": "markdown", - "id": "8a03d0c5", + "id": "177b5e84", + "metadata": { + "editable": true + }, + "source": [ + "## Same code but now with momentum gradient descent" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "755d4596", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ + "# Using Autograd to calculate gradients using SGD\n", + "# OLS example\n", + "from random import random, seed\n", + "import numpy as np\n", + "import autograd.numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from autograd import grad\n", + "\n", + "# Note change from previous example\n", + "def CostOLS(y,X,theta):\n", + " return np.sum((y-X @ theta)**2)\n", + "\n", + "n = 100\n", + "x = 2*np.random.rand(n,1)\n", + "y = 4+3*x+np.random.randn(n,1)\n", + "\n", + "X = np.c_[np.ones((n,1)), x]\n", + "XT_X = X.T @ X\n", + "theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)\n", + "print(\"Own inversion\")\n", + "print(theta_linreg)\n", + "# Hessian matrix\n", + "H = (2.0/n)* XT_X\n", + "EigValues, EigVectors = np.linalg.eig(H)\n", + "print(f\"Eigenvalues of Hessian Matrix:{EigValues}\")\n", + "\n", + "theta = np.random.randn(2,1)\n", + "eta = 1.0/np.max(EigValues)\n", + "Niterations = 100\n", + "\n", + "# Note that we request the derivative wrt third argument (theta, 2 here)\n", + "training_gradient = grad(CostOLS,2)\n", + "\n", + "for iter in range(Niterations):\n", + " gradients = (1.0/n)*training_gradient(y, X, theta)\n", + " theta -= eta*gradients\n", + "print(\"theta from own gd\")\n", + "print(theta)\n", + "\n", + "\n", + "n_epochs = 50\n", + "M = 5 #size of each minibatch\n", + "m = int(n/M) #number of minibatches\n", + "t0, t1 = 5, 50\n", + "def learning_schedule(t):\n", + " return t0/(t+t1)\n", + "\n", + "theta = np.random.randn(2,1)\n", + "\n", + "change = 0.0\n", + "delta_momentum = 0.3\n", + "\n", + "for epoch in range(n_epochs):\n", + " for i in range(m):\n", + " random_index = M*np.random.randint(m)\n", + " xi = X[random_index:random_index+M]\n", + " yi = y[random_index:random_index+M]\n", + " gradients = (1.0/M)*training_gradient(yi, xi, theta)\n", + " eta = learning_schedule(epoch*m+i)\n", + " # calculate update\n", + " new_change = eta*gradients+delta_momentum*change\n", + " # take a step\n", + " theta -= new_change\n", + " # save the change\n", + " change = new_change\n", + "print(\"theta from own sdg with momentum\")\n", + "print(theta)" + ] + }, + { + "cell_type": "markdown", + "id": "66103f4e", "metadata": { "editable": true }, @@ -4289,8 +4453,8 @@ }, { "cell_type": "code", - "execution_count": 30, - "id": "885dd6a7", + "execution_count": 32, + "id": "b1f5ba7c", "metadata": { "collapsed": false, "editable": true @@ -4334,7 +4498,7 @@ }, { "cell_type": "markdown", - "id": "3168e44a", + "id": "152c46ee", "metadata": { "editable": true }, @@ -4352,8 +4516,8 @@ }, { "cell_type": "code", - "execution_count": 31, - "id": "592c4321", + "execution_count": 33, + "id": "34a7c3a4", "metadata": { "collapsed": false, "editable": true diff --git a/doc/src/week39/.ipynb_checkpoints/test-checkpoint.ipynb b/doc/src/week39/.ipynb_checkpoints/test-checkpoint.ipynb new file mode 100644 index 000000000..24babda28 --- /dev/null +++ b/doc/src/week39/.ipynb_checkpoints/test-checkpoint.ipynb @@ -0,0 +1,532 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "14315bd2", + "metadata": {}, + "source": [ + "\n", + "" + ] + }, + { + "cell_type": "markdown", + "id": "fe75bd92", + "metadata": {}, + "source": [ + "# Codes\n", + "**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State University\n", + "\n", + "Date: **Sep 30, 2022**\n", + "\n", + "Copyright 1999-2022, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "62131698", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + ">0 f([0.74724774]) = 0.55838\n", + ">1 f([0.59779819]) = 0.35736\n", + ">2 f([0.47823856]) = 0.22871\n", + ">3 f([0.38259084]) = 0.14638\n", + ">4 f([0.30607268]) = 0.09368\n", + ">5 f([0.24485814]) = 0.05996\n", + ">6 f([0.19588651]) = 0.03837\n", + ">7 f([0.15670921]) = 0.02456\n", + ">8 f([0.12536737]) = 0.01572\n", + ">9 f([0.10029389]) = 0.01006\n", + ">10 f([0.08023512]) = 0.00644\n", + ">11 f([0.06418809]) = 0.00412\n", + ">12 f([0.05135047]) = 0.00264\n", + ">13 f([0.04108038]) = 0.00169\n", + ">14 f([0.0328643]) = 0.00108\n", + ">15 f([0.02629144]) = 0.00069\n", + ">16 f([0.02103315]) = 0.00044\n", + ">17 f([0.01682652]) = 0.00028\n", + ">18 f([0.01346122]) = 0.00018\n", + ">19 f([0.01076897]) = 0.00012\n", + ">20 f([0.00861518]) = 0.00007\n", + ">21 f([0.00689214]) = 0.00005\n", + ">22 f([0.00551372]) = 0.00003\n", + ">23 f([0.00441097]) = 0.00002\n", + ">24 f([0.00352878]) = 0.00001\n", + ">25 f([0.00282302]) = 0.00001\n", + ">26 f([0.00225842]) = 0.00001\n", + ">27 f([0.00180673]) = 0.00000\n", + ">28 f([0.00144539]) = 0.00000\n", + ">29 f([0.00115631]) = 0.00000\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtNUlEQVR4nO3dd3wUdf7H8ddnNz0hhJDQEtKooqBAaCqIHcuBFREbpyfq2fXuxDvP7u/0zn7qKWA5PQW7ouKhNAENJYDSISFACiWBhAAhPd/fH7vxVkzIArs7u5vP8/HgQXZmsvNmEt6ZzMx3RowxKKWUCnw2qwMopZTyDC10pZQKElroSikVJLTQlVIqSGihK6VUkAixasUJCQkmLS3NqtUrpVRAWr58+W5jTGJT8ywr9LS0NLKzs61avVJKBSQR2dbcPD3kopRSQUILXSmlgoQWulJKBQktdKWUChJa6EopFSRaLHQReUNEikVkTTPzRUReFJFcEVklIgM8H1MppVRL3NlDfwsYdZj55wE9nH8mAv869lhKKaWOVIuFboxZAJQeZpExwNvGYTEQJyKdPRXwUD8V7OXJrzegt/1VSgWahgbDE1+tY932fV55f08cQ08CClxeFzqn/YqITBSRbBHJLikpOaqVrSrcy6vfbWZVYflRfb5SSlklK28PUxZuYeMu/y10txljJhtjMo0xmYmJTY5cbdGY/klEhtp5b0m+h9MppZR3vbckn7ioUM47wTsHMTxR6EVAV5fXyc5pXhEbEcqYk7ow46ftlFfWems1SinlUcX7q5i1dieXDUgmItTulXV4otBnANc6r3YZCpQbY3Z44H2bddWQVCpr6/lspdd+biillEd9mF1IXYPhyiEpXluHO5ctTgOygF4iUigiN4jIzSJys3ORmUAekAtMAX7vtbROfZPb0i+5Le8u2aYnR5VSfq++wfDeknxO7taebokxXltPi3dbNMZc2cJ8A9zqsURuGj84hUmfrGb5tjIy0+J9vXqllHLbgpwSivZWcv/5vb26noAdKfqbE7vQJjyEd/XkqFLKz727OJ+EmDDO6dPJq+sJ2EKPDg/h4gFJfLV6B2UVNVbHUUqpJm3fW8ncDbsYm9mVsBDvVm7AFjrA+CEp1NQ18PGKQqujKKVUk95fVoABrhzsvZOhjQK60Ht3imVgajveXZKvJ0eVUn6nrr6B6cvyGdEjka7xUV5fX0AXOsBVQ1LYsruCrM17rI6ilFK/MGdDMbv2VXOVFy9VdBXwhX5+387ERYXqyVGllN95d0k+nWIjOKN3B5+sL+ALPSLUzqUDkpm1dicl+6utjqOUUgDk7znIwpwSrhjUlRC7b6o24AsdHCdH6xoMH2QXtLywUkr5wLRl+QgwbnDXFpf1lKAo9G6JMQzLaM+0pfk0NOjJUaWUtWrqGvgwu4Azj+tI57aRPltvUBQ6OPbSC8sqWZBzdLflVUopT/lm3U52H6hhvI9OhjYKmkI/9/hOtI8O05OjSinLvbs4n+R2kYzocXS3CT9aQVPoYSE2xg7qypz1u9hRXml1HKVUK5VbfICsvD1cOTgFu018uu6gKXSAKwelYIDpS/XkqFLKGtOW5hNiE8Zm+u5kaKOgKvSU9lEM75HI+8sKqKtvsDqOUqqVqaqt5+MVhZx7fCcS24T7fP1BVejgGDm6c18VczcUWx1FKdXKzFy9g70Ha302MvRQQVfoZ/buQMfYcD05qpTyuXeX5JOREM2wbu0tWX/QFXqI3ca4QSksyCmhoPSg1XGUUq3Ehp37WL6tjPFDUhDx7cnQRkFX6OAYmSU4Tk4opZQvvLckn7AQG5cOSLYsQ1AWeue2kZzRuyMfZBdQU6cnR5VS3lVRXccnK4q4oG9n2kWHWZYjKAsd4KqhKew+UMM363ZaHUUpFeS++Gk7B6rrLDsZ2ihoC31Ej0SS4iJ5T0+OKqW87L2l+fTsGMPA1HaW5gjaQrfbhPFDUvhh8x42lxywOo5SKkitKtzLqsJyrhqSatnJ0EZBW+gAl2cmE2ITpuleulLKS95bkk9kqJ2LByRZHSW4C71DmwjOPb4TH60opKq23uo4Sqkgs6+qlhk/bWf0iV2IjQi1Ok5wFzo4bqu792AtX6/ZYXUUpVSQ+XxlEQdr6n1+m9zmBH2hD8toT3pCNO8u1sMuSinPMcbw7pJ8TkiKpV9yW6vjAK2g0G02YfzgFLK3lbFx536r4yilgsSK/DI27NzvFydDGwV9oQNcOjCZMLuN95ZsszqKUipIvLskn5jwEEaf2MXqKD9rFYUeHx3G+X078cmKIg7W1FkdRykV4PYerOHLVTu4qH8XosNDrI7zs1ZR6ABXDU1lf3UdX/y03eooSqkA99HyQmrqGhg/ONXqKL/Qago9M7UdPTvG6G11lVLHxBjDe0vzGZASR58usVbH+QW3Cl1ERonIRhHJFZFJTcxPEZF5IrJSRFaJyPmej3psRBwnR1cVlrOqcK/VcZRSASpr8x7ySioYP8S/9s7BjUIXETvwMnAe0Ae4UkT6HLLYA8AHxpj+wDjgFU8H9YRLBiYTEx7C1IVbrI6ilApQUxdtoX10GBf262x1lF9xZw99MJBrjMkzxtQA04ExhyxjgMbfPdoCfnmgOjYilCsHd+Wr1TsoLNOHXyiljkzOrv3M3VDMtcPSiAi1Wx3nV9wp9CSgwOV1oXOaq4eBq0WkEJgJ3N7UG4nIRBHJFpHskpKSo4h77H57SjoCvPn9VkvWr5QKXFMXbiE8xMY1w/zvcAt47qTolcBbxphk4HzgHRH51XsbYyYbYzKNMZmJiYkeWvWR6RIXyYX9OjN9aT7llbWWZFBKBZ7i/VV8urKIyzOTibfwIRaH406hFwFdXV4nO6e5ugH4AMAYkwVEAAmeCOgNvxueQUVNvT6iTinltrd/2EZtQwM3nJphdZRmuVPoy4AeIpIuImE4TnrOOGSZfOBMABE5DkehW3NMxQ0nJLXllO7tefP7LfqIOqVUiw7W1PHO4m2c06cj6QnRVsdpVouFboypA24DZgHrcVzNslZEHhWR0c7F7gVuFJGfgGnABGOM8VZoT7hxeAa79lXrQCOlVIs+zC6kvLKWiSP8d+8cwK0xq8aYmThOdrpOe9Dl43XAKZ6N5l2n9UykV8c2TFmYxyUDkvzm5jpKKf9S32CYuiiPASlxDEyNtzrOYbWakaKHEhF+NzydDTv3szBnt9VxlFJ+atbanRSUVvr93jm04kIHGH1SFzq0CWfKwjyroyil/JAxhtcW5JHaPoqz+3SyOk6LWnWhh4fYmXBKGgtzdrNu+z6r4yil/Ez2tjJ+KtjL705Nx27z/8OyrbrQAa4anEpUmJ2pi3QvXSn1S1MW5NEuKpTLBnZteWE/0OoLvW1UKFcM6sqMH7ezo7zS6jhKKT+RV3KAb9fv4pqhqUSG+d8w/6a0+kIHuP6UdBqM4a0ftlodRSnlJ15ftIVQu41rhqVZHcVtWuhA1/gozu/bmfcW57O/Sm8HoFRrt+dANR8tL+TSAUkktgm3Oo7btNCdJo7IYH91He8vK2h5YaVUUHtn8Taq6/x7mH9TtNCd+iXHMTg9nje/30ptvd4OQKnWqqq2nreztnFm7w507xBjdZwjooXuYuLwDIr2VjJz9Q6royilLPLxikJKK2q4MQAGEh1KC93FGb07kJEYzZSFefj5rWiUUl7Q0GCYunAL/ZLbMiTdv4f5N0UL3YXNJtw4PIM1RfvIyttjdRyllI/NXr+LLbsruHF4RkDe30kL/RAX908iISaMKQt0oJFSrc2UhXkkxUVy3gn+P8y/KVroh4gItXPtsDTmbSxh0679VsdRSvnIivwylm0t44ZT0wmxB2Y1BmZqL7t6aCoRoTam6k27lGo1pi7MIzYihLGDAmOYf1O00JsQHx3G5QO78tnK7RTvq7I6jlLKy7btqeC/a3Zy1dBUYsLdekyEX9JCb8YNp6ZT29DAv7O2Wh1FKeVlbyzagt0mTDg5zeoox0QLvRlpCdGc26cT/1mcT0V1ndVxlFJeUlZRwwfZhYw5KYmOsRFWxzkmWuiHceOIDMora/kwW28HoFSwenfJNipr67lxeOANJDqUFvphDExtx8DUdrz+/Rbq9HYASgWdqtp63vphm+MZw53aWB3nmGmht+DG4RkUlFYya+0uq6MopTzs8x+L2H2gOiCeF+oOLfQWnN2nI2nto5i8YLPeDkCpINLQYJiycAt9Osdycrf2VsfxCC30Fthtwg3DM/ipsJxlW8usjqOU8pD5m4rJLT7AxBGBOcy/KVrobrhsQDLtokKZrLcDUCpoTF6QR+e2EVzQr7PVUTxGC90NkWF2rhmWxpwNu9hccsDqOEqpY7S6sJzFeaVcf0o6oQE6zL8pwfMv8bJrh6USZrfxyrzNVkdRSh2jl+fl0iY8hCsGB+4w/6ZoobspISacq4em8unKQt1LVyqArSkq579rd3L9qenERoRaHcejtNCPwC0juxEeYueF2TlWR1FKHaVnv91E28hQbhiebnUUj9NCPwIJMeFMOCWNL1ZtZ+NOvbWuUoFmRX4ZczcUM3FERtDtnYMW+hGbODyD6LAQnvt2k9VRlFJH6NlvNtE+Oizgb8LVHLcKXURGichGEckVkUnNLDNWRNaJyFoRec+zMf1Hu+gwbjg1nf+u3cmaonKr4yil3LQkbw+Lcndzy8huRAfwLXIPp8VCFxE78DJwHtAHuFJE+hyyTA/gfuAUY8zxwF2ej+o/bhieTtvIUN1LVypAGGN45ptNdGjjuLghWLmzhz4YyDXG5BljaoDpwJhDlrkReNkYUwZgjCn2bEz/EhsRysQRGczZUMyKfB09qpS/W5S7m6VbS7n19O5EhNqtjuM17hR6EuB6/9hC5zRXPYGeIvK9iCwWkVFNvZGITBSRbBHJLikpObrEfmLCyWnER4fpXrpSfq5x77xL2wjGBdl154fy1EnREKAHMBK4EpgiInGHLmSMmWyMyTTGZCYmJnpo1daIDg/hltO6sTBnN0vy9lgdRynVjLkbivmxYC+3n9mD8JDg3TsH9wq9CHD9sZbsnOaqEJhhjKk1xmwBNuEo+KB29dBUOrQJ55lvN+mdGJXyQ8YYnv12EynxUVw2MNnqOF7nTqEvA3qISLqIhAHjgBmHLPMZjr1zRCQBxyGYoL+TVWSYnVtP787SLaV8n6t76Ur5m1lrd7J2+z7uPLNHUN2zpTkt/guNMXXAbcAsYD3wgTFmrYg8KiKjnYvNAvaIyDpgHvBHY0yraLhxg7vSpW0ET3+zUffSlfIj9Q2OvfOMxGgu6n/oab/g5NaPLGPMTGNMT2NMN2PME85pDxpjZjg/NsaYe4wxfYwxfY0x070Z2p+Eh9i5/cwe/Fiwl3kbg/riHqUCypertrNp1wHuPqsndltw3O+8JcH/O4gPXDYwmZT4KJ75Ro+lK+UP6uobeH52Dr07teGCvsFzv/OWaKF7QKjdxh1n9mDt9n3MWrvT6jhKtXqLn3uTq6c/x8MJ+7C1kr1z0EL3mItO6kJGYjTPfruJ+gbdS1fKKrULv2fYnyZy/fIZDJk4FrKyrI7kM1roHhJit3HXWT3ZtOsAX67abnUcpVqtrc++gt00IIDU1MD8+VZH8hktdA+6sG9nenVswwuzc6irb7A6jlKtTlV1LSGLFmIAY7dDWBiMHGl1LJ/RQvcgm024++ye5O2u4LMfdS9dKV9b9tgLpO8uoOC2PyCPPQZz5sCwYVbH8hmx6qqMzMxMk52dbcm6vckYw29eWkR5ZS1z7x3ZKgYzKOUPKkvLOZCWwd72neieuwqxB+cwfxFZbozJbGqeto2HiQj3nt2LgtJKPswutDqOUq3G+rv+TOL+UmqeeTZoy7wlWuheMLJXIv1T4vjn3ByqauutjqNU0KvYtJnjp03hhyHncPwl51odxzJa6F7QuJe+o7yK6UvzrY6jVNDbftPtNGCj7QvPWB3FUlroXnJK9/YMSY/n5fmbqazRvXSlvOXAnPn0mP81315wDccPOcHqOJbSQvcSEeHec3pRsr+adxZvtTqOUsGpoYEDt9zG9jYJdH/mUavTWE4L3YsGp8czvEcCr36Xx4HqOqvjKBV0Kqa8Qaectfz3qjvp06313LOlOVroXnbP2T0prajh3z9stTqKUsFl/37M/fezoksvTn3oDqvT+AUtdC/rn9KOM3t34LXvNlNeWWt1HKWCRsWjjxNTtpv5t/yZnp1irY7jF7TQfeDus3uyr6qOqQuD/iFOSvnGli2Ev/A8nx5/OhfddInVafyGFroPnJDUlgv7dWbygjwKyw5aHUepgLf/rnupQVh/x/1kJMZYHcdvaKH7yP3nH4cIPPHVequjKBXQzHff0WbGp7x5ylgmjh9hdRy/ooXuI0lxkdx2ene+XrOThTklVsdRKjDV17PvltsoapNI7AOTSIgJtzqRX9FC96HfDc8gtX0UD89YS02d3l5XqSNV8/qbtF2/hncuupkrR/a2Oo7f0UL3oYhQOw9e2IfNJRV6GaNSR2r/fmon3c/yLr0Z+dAdhOidTH9Ft4iPnXlcR87o3YHnZ2+ieF+V1XGUChjlf32E6LLdzLv5zwztlmB1HL+khW6BBy/sQ2294cmvN1gdRanAsGULUS+/yOd9z+TqOy63Oo3f0kK3QFpCNDeOSOeTlUUs21pqdRyl/F7xzXdQKzbK//oIndpGWB3Hb2mhW+TW07vTuW0ED32+lvoGa54apVQgqJk7lw7ffMn008dzxcVDrY7j17TQLRIVFsJfLjiOdTv28Z7eM12pptXXU36T4zLF7n9/iPCQ1vkkIndpoVvogr6dGZbRnqdnbaS0osbqOEr5nfJXJpOYu56vrr6LESemWh3H72mhW0hEeGTM8RyoruPpbzZaHUcp/7JvH/LAAyxPPo7zHrvT6jQBQQvdYj07tuG6YWlMW5rP6sJyq+Mo5Te2//EBYveVsuG+R+naPtrqOAFBC90P3HV2D9pHh/HgjDU06AlSpajblEPC668yc8A5XDLxYqvjBAwtdD8QGxHKfaN6szJ/L5+sLLI6jlKWK7rxNmptdiL/8SSRYXoi1F1uFbqIjBKRjSKSKyKTDrPcpSJiRCTTcxFbh0sHJNM/JY4nv17Pvip9EIZqvcpnziJ1wTfMvOA6Rp5+ktVxAkqLhS4iduBl4DygD3CliPRpYrk2wJ3AEk+HbA1sNuHR0Sewp6KGF2bnWB1HKWvU11Px+zspbNuBAc8/hohYnSiguLOHPhjINcbkGWNqgOnAmCaWewx4CtAblBylvsltGTcohbd+2MqmXfutjqOUz+U//RJdtm1k2S330S010eo4AcedQk8CClxeFzqn/UxEBgBdjTFfHe6NRGSiiGSLSHZJid4TvCl/PLcXMeEhPPT5WozRE6Sq9Wgo20vs4w/zY+oJnPXgbVbHCUjHfFJURGzAs8C9LS1rjJlsjMk0xmQmJupP36bER4fxh3N7kZW3h5mrd1odRymf2XjHJGIPlFP2xN9pExlmdZyA5E6hFwFdXV4nO6c1agOcAMwXka3AUGCGnhg9euMHp9CncyyPf7WOgzV1VsdRyuv2r15P92mv892w8xg5fpTVcQKWO4W+DOghIukiEgaMA2Y0zjTGlBtjEowxacaYNGAxMNoYk+2VxK2A3SY8OuZ4dpRX8cq8zVbHUcrrCm+8jRpbCF1eflZPhB6DFgvdGFMH3AbMAtYDHxhj1orIoyIy2tsBW6vMtHgu7p/E5AV5bN1dYXUcpbxm24dfcNySufwwdiK9+veyOk5AE6tOvGVmZprsbN2JP5zifVWc/vR8hmS0540Jg6yOo5THmbo68tN6E1pxgKjNm4iLj7U6kt8TkeXGmCYPaetIUT/WITaCu87qydwNxcxZv8vqOEp5VlYWO86+kNSizeT96SEtcw/QQvdz152cRrfEaB7+Yi0V1XqCVAWJrCzMGWfQef4sGkQYNryf1YmCgha6nwsLsfG3S/pRWFbJ41+tszqOUh5h5s2DqioEx22k7QsXWB0pKGihB4DB6fHcNKIb05YW8O06PfSiAt/6/N0IYESQ8HAYOdLqSEFBCz1A3HN2T/p0jmXSx6so2V9tdRyljtqOH9fR9a3X2JrUDfPoozBnDgwbZnWsoKCFHiDCQmw8P+4k9lfXcd/Hq/S2ACog1dfUsu/SKwAI+/ILbA88oGXuQVroAaRnxzbcf15v5m4o1gdLq4C0cuK99Mpbw9q/PkWXk46zOk7Q0UIPMNcNS2N4jwQe/3I9eSUHrI6jlNu2fP4N/d9+mcWnnM+QP99qdZygpIUeYGw24R+XnUhYiI27P/iJ2voGqyMp1aKq3aVE/vZadsR1pOeH/9bh/V6ihR6AOrWN4P8u7stPBXt5aW6u1XGUalHOZdeSsLeE4n9NJb5zgtVxgpYWeoC6oF9nLumfxEvzclmRX2Z1HKWatfHpf9H3u6/47srfM+CK862OE9S00APYw2OOp1NsBHe//6OOIlV+ad/aDSQ/8AdWpffj5KnPWB0n6GmhB7DYiFCeHXsi+aUHdRSp8jumpoY9F42lHiFs2rtE6kMrvE4LPcANyWivo0iVX9p4259Iz11N1p+eoPeQE6yO0ypooQcBHUWq/E3JzNn0mPpP5g09j7MevsPqOK2GFnoQcB1FOklHkSqLNZSWwTVXUxTXkR7T38Ru00sUfUULPUj07NiGSaN6M2dDMdOWFlgdR7VWxpB3+bW0Kyth07Ovkpza0epErYoWehCZcHIap3ZP4LEv1+koUmWJohcn033ul3x56c2ced1vrI7T6mihBxGbTXj6ch1FqqxRvX4j7e67m+Vp/Rjx+tM6GtQCWuhBplPbCJ64+AQdRap8q7aW3RddTi02qt98i/jYSKsTtUpa6EHown5dfh5FulJHkSofKLjzPpI2rea/dzzCySP7Wx2n1dJCD1I6ilT5yoFv5pL06vN8PWgUox+/y+o4rZoWepBqHEW6rfQgj3+13uo4KkiZ0lJqxo8nP64zKe9MJTLMbnWkVk0LPYgNyWjPxBEZTFuaz4fZeimj8jBj2Hr5tbQp283Sx1/k+F5JVidq9bTQg9wfzunFqd0T+POnq8navMfqOCqIrHvyRdLnfsWMS2/m0psvsTqOQgs96IXabbx81QDS2kdz83+Ws1mvT1cekLNoBWkPTWJV9/6c9+9ndDSon9BCbwXaRobyxoRBhNiE699aRmlFjdWRVADb8/VsYi84FyM2Os/4gKjIcKsjKSct9Faia3wUU67LZGd5FRPfzqaqtt7qSCoAVc5fSNyFo+i4bzeRpo7EvSVWR1IutNBbkQEp7Xh27ElkbyvjPr2JlzpC9Q2G3Nv/hL3BsTNga6iH+fOtDaV+QQu9lbmgX2f+eG4vPv9xO8/NzrE6jgogM/7wFH3XLKbBZgO7HcLCYORIq2MpF24VuoiMEpGNIpIrIpOamH+PiKwTkVUiMkdEUj0fVXnK70d2Y2xmMi/OyeGTFYVWx1EBYPZz/+bCFx4gr98QbHPnwmOPwZw5MGyY1dGUi5CWFhARO/AycDZQCCwTkRnGGNdnnq0EMo0xB0XkFuDvwBXeCKyOnYjw+EV9KSit5L6PV5EUF8mQjPZWx1J+atnH3zLsvpvZkZRB6nezIK4tnHaa1bFUE9zZQx8M5Bpj8owxNcB0YIzrAsaYecaYg86Xi4Fkz8ZUnhYWYuPVqweSEh/FTf9ZzpbdFVZHUn4oN+tH0q8by4E2cSQsmIM9rq3VkdRhuFPoSYDrMMNC57Tm3AB83dQMEZkoItkikl1SomfHrdY2KpQ3JwzGJsJv31xKmV7OqFyUbM4n/DcXEtpQj23WLKLSulodSbXAoydFReRqIBP4R1PzjTGTjTGZxpjMxMRET65aHaWU9lFMuXYg28uruOmd5VTX6eWMCg6W7qXsjFEklO9mz/SPSczsZ3Uk5QZ3Cr0IcP3RnOyc9gsichbwF2C0MUafVBxABqbG84/L+rF0aymTPl6tlzO2cvXVNWw+/QK6FWxk0/OTyRh9ttWRlJvcKfRlQA8RSReRMGAcMMN1ARHpD7yGo8yLPR9TeduYk5K45+yefLqyiBfn6IMxWi1jWHPhFfRd9QNZf3icE2+91upE6gi0eJWLMaZORG4DZgF24A1jzFoReRTINsbMwHGIJQb40PnYqXxjzGgv5lZecPsZ3dm6p4LnZm8iLSGKMSfp3fNamzXX38GJsz9j3hW3MPKpX12hrPxci4UOYIyZCcw8ZNqDLh+f5eFcygIiwt8u6UthWSV//HAVXeIiGZQWb3Us5SObHv47J7z1EvNHjGH4f/6pzwQNQDpSVP1CeIid164eSFK7SCa+nc1WvZyxVSh64z26PzKJxcefTObM6YSE6IMqApEWuvqVdtFhvDFhEAb07oytQNmsOSTcNIG1XXuT+s0MYqIjrI6kjpIWumpSekI0k6/JpHBvJWNfy2L73kqrIykvKPo+m5CLL2J7bAfsX35J5y46YjiQaaGrZg1Oj+ft6wezq7yKS//1Azm79lsdSXnQ1tffI/bMkdQhVH4xkz79MqyOpI6RFro6rKEZ7Xn/pmHUNRguezWL5dvKrI6kPGDd5HdJ+d3VxFRX0Lahhj6i50qCgRa6alGfLrF8csvJtIsK5aqpi5m7YZfVkdQxmP/l93S582YEgwC2ujq9r3mQ0EJXbukaH8VHt5xMjw5tuPHt5Xy8XG+7G4hmvv45fceeT4gA4eF6X/Mgo4Wu3JYQE860iUMZmhHPvR/+xOQFm62OpNxkjOHTJ6Yw8pYrqI9pQ+jybGTePL2veZBxa2CRUo1iwkN4Y8Ig7vngJ/5v5gZ2H6hh0qje2PSp736rvsHw+R2PMfqVR9iZ1otOi+YQ0qWzY6YWeVDRQldHLDzEzj/H9SchOozJC/LYfaCapy7tR6hdf+HzN1U1dXwz/nYu+fhVtgw4hbR5M5HYWKtjKS/RQldHxWYTHh59PAkx4Tzz7SbKKmp4+aoBRIXpt5S/2FdRxZLzxzN6wafknnsx3WdMdxwvV0FLd6nUURMRbj+zB3+7pC/fbSrhqqlL9CEZfqK4uIw1w87h7AWfkjPh93T/+mMt81ZAC10dsysHp/DKVQNZu30fl+uoUssV5Bawc8gIhq5eRM4D/0ePN18GvdFWq6CFrjxi1AmddFSp1bKyKPn9XdgGDaJ3wUa2vfIGPR673+pUyoe00JXHuI4qvfw1HVXqU1lZNIw8nYR/vUCXvbsof+Ip0m+ZYHUq5WNa6MqjGkeVxkU6RpV+8dN2qyMFPdPQQNGkh5CaagTAZiexocrqWMoCWujK4xpHlfbuFMvt01Zy8zvLKd6nBeMNO3Lz+XHgSJIWfIsRG8ZuR8J15GdrpYWuvCIhJpyPbh7GpPN6M29jMWc9+x0fZBfoA6g9pKHBMPv5twnpfxJ9Vi9myZ0PwoIFiI78bNXEqv9gmZmZJjs725J1K9/KKznApI9Xs3RrKad2T+Bvl/Sla3yU1bECVl7hHtZNuJUL57xPYZd0QqZPo9PwIVbHUj4iIsuNMZlNzdM9dOV1GYkxTJ84lMcuOoGV+WWc89wC3vx+C/UNurd+JOrqG5j+1n+pHjiYC+e8T+7l15GUs0bLXP1MC135hM0mXDM0lW/uOY0hGfE88sU6Ln/1B3KL9fJGd+TN+IZl/U/j4t/9huTKUva+/zHdP3gLidLfdNT/aKErn0qKi+TNCYN47ooT2bK7gvNfWMRLc3OorW+wOppfqqqt57O/vkjqRecxbPUiwkwDbd59m7ixl1gdTfkhLXTlcyLCxf2T+fae0zjn+I48/c0mfvPPRawuLLc6ml9ZuXYbn4y6lgufuBubcfzAExFYs8biZMpfaaEryyTEhPPS+AFMvmYgpRU1XPTK9zz59QaqauutjmapisoaPrvzCZKG9Gf83PfYd+pIJCJCH0ahWqS3xlOWO+f4TgzJaM/fZq7n1e82M2vtTp68pC9DMlrRE+izsmD+fDbVhNAwZSoXFW2ioNeJVE75gvjhp/w8n5Ej9ZJE1Sy9bFH5le9zdzPpk1UUlFZycf8krh6awoCUdo5DDUHK/PAD5owzkGrHSM/S6Dj2PvEUGXfcqDfVUr9yuMsWdQ9d+ZVTuicw664RPD87h3cXb+PTlUV07xDDFZlduWRAEu1jwq2O6DHF+6pY+M4XDHjsj6Q5y7xBhDb33kn8nROtjqcCkO6hK79VUV3HV6t2MH1ZPivy9xJqF87u05GxmV0Z3iMRewA+9q5u0fds/egrftyxn7Qf5pJZuI6D4VFE1NUgGCQsTEd6qsM63B66FroKCDm79vP+sgI+WVlEaUUNXdpGcHlmVy7PTCa5nf9fi71tTwXLXn6HMY/cTkhDHQIcaJdA1R/+RMIdt8Dq1XqMXLlFC10Fjeq6emavK+b97AIW5pQAcGr3BMYNSuGsPh0ID7FbnPB/qmrryX7rE2TKZMK2FzGgaAN2HP/fjM2GPPoo/OUvFqdUgUaPoaugER5i54J+nbmgX2cKyw7yYXYhH2YXcOt7K4iPDuOS/klcMagrPTq2sSZgVhYl0z5m3c79hC9bwilbf0IAA1SfeTb27xdCba3j0MoZZ1iTUQUtt/bQRWQU8AJgB6YaY548ZH448DYwENgDXGGM2Xq499Q9dOUp9Q2GRbm7eX9ZPt+u20VtvaF3pzZ06xBD13ZRpMRH0TU+kq7tougSF0lYyDEMv3BePlgXF8++op1s6zuIsoIddPjqM2zFxfRevwwxBgGqwyIIq6ly3KPcbofHHnMcUtFDK+oYHNMhFxGxA5uAs4FCYBlwpTFmncsyvwf6GWNuFpFxwMXGmCsO975a6Mob9hyo5pMVRSzIKaGg9CCFZZXUudwEzCbQuW0kye0i6Rr/y7LvGh9FYkw4tox0GrZtw4hQHxtHTVQ0u1K6syu6HUPmf440NPy81y3OPzhf43xtbHbkxt/B229DTY1jQJCe7FQecKyFPgx42BhzrvP1/QDGmL+5LDPLuUyWiIQAO4FEc5g310JXvlDfYNi5r4qC0oPklx6ksPQgBWWVP78u3l/9i+UXvTKBpP27f/U+h15P01jmriWOzeb4Y8z/Chx0j1x51LEeQ08CClxeFwKH3q/z52WMMXUiUg60B37xP0NEJgITAVJSUtwKr9SxsNuEpLhIkuIiGdrEyNOq2noKyyopKHOUfZenHN+yLV0Q2VjgEhrqKPD6ekeJP/887NnzywLXIlc+4tOTosaYycBkcOyh+3LdSjUlItRO9w4xdO8Q45iQmgrbtjX/CSIQEgJ33w1xcf+7r4ruhSs/4E6hFwFdXV4nO6c1tUyh85BLWxwnR5UKLFu3Qlqao9RFoF07aNMG+veH88779d53Iy1y5QfcKfRlQA8RScdR3OOA8YcsMwO4DsgCLgPmHu74uVJ+betWqxModVRaLHTnMfHbgFk4Llt8wxizVkQeBbKNMTOA14F3RCQXKMVR+koppXzIrWPoxpiZwMxDpj3o8nEVcLlnoymllDoS+oALpZQKElroSikVJLTQlVIqSGihK6VUkLDs9rkiUgIcZgTHYSVwyChUP6G5jozmOnL+mk1zHZljyZVqjElsaoZlhX4sRCS7uXsZWElzHRnNdeT8NZvmOjLeyqWHXJRSKkhooSulVJAI1EKfbHWAZmiuI6O5jpy/ZtNcR8YruQLyGLpSSqlfC9Q9dKWUUofQQldKqSDht4UuIpeLyFoRaRCRZi/vEZFRIrJRRHJFZJLL9HQRWeKc/r6IhHkoV7yIfCsiOc6/2zWxzOki8qPLnyoRucg57y0R2eIy7yRf5XIuV++y7hku063cXieJSJbz671KRK5wmefR7dXc94vL/HDnvz/XuT3SXObd75y+UUTOPZYcR5HrHhFZ59w+c0Qk1WVek19TH+WaICIlLuv/ncu865xf9xwRuc7HuZ5zybRJRPa6zPPm9npDRIpFZE0z80VEXnTmXiUiA1zmHfv2Msb45R/gOKAXMB/IbGYZO7AZyADCgJ+APs55HwDjnB+/CtzioVx/ByY5P54EPNXC8vE4bikc5Xz9FnCZF7aXW7mAA81Mt2x7AT2BHs6PuwA7gDhPb6/Dfb+4LPN74FXnx+OA950f93EuHw6kO9/H7sNcp7t8D93SmOtwX1Mf5ZoAvNTE58YDec6/2zk/buerXIcsfzuO2357dXs533sEMABY08z884GvcTzFcCiwxJPby2/30I0x640xG1tYbDCQa4zJM8bUANOBMSIiwBnAR87l/g1c5KFoY5zv5+77XgZ8bYw56KH1N+dIc/3M6u1ljNlkjMlxfrwdKAaaHAl3jJr8fjlM3o+AM53bZwww3RhTbYzZAuQ6388nuYwx81y+hxbjeHKYt7mzvZpzLvCtMabUGFMGfAuMsijXlcA0D637sIwxC3DswDVnDPC2cVgMxIlIZzy0vfy20N3U1AOsk3A8oHqvMabukOme0NEYs8P58U6gYwvLj+PX30xPOH/dek5Ewn2cK0JEskVkceNhIPxoe4nIYBx7XZtdJntqezX3/dLkMs7t0fjAc3c+15u5XN2AYy+vUVNfU1/mutT59flIRBofV+kX28t5aCodmOsy2Vvbyx3NZffI9vLpQ6IPJSKzgU5NzPqLMeZzX+dpdLhcri+MMUZEmr3u0/mTty+Opz01uh9HsYXhuBb1PuBRH+ZKNcYUiUgGMFdEVuMoraPm4e31DnCdMabBOfmot1cwEpGrgUzgNJfJv/qaGmM2N/0OHvcFMM0YUy0iN+H47eYMH63bHeOAj4wx9S7TrNxeXmVpoRtjzjrGt2juAdZ7cPwqE+Lcy2rqwdZHlUtEdolIZ2PMDmcBFR/mrcYCnxpjal3eu3FvtVpE3gT+4Mtcxpgi5995IjIf6A98jMXbS0Riga9w/DBf7PLeR729mnAsDzx353O9mQsROQvHD8nTjDHVjdOb+Zp6oqBazGWMcX0Y/FQc50waP3fkIZ873wOZ3MrlYhxwq+sEL24vdzSX3SPbK9APufz8AGtxXJUxDphhHGcZ5uE4fg2OB1h7ao+/8YHY7rzvr47dOUut8bj1RUCTZ8O9kUtE2jUeshCRBOAUYJ3V28v5tfsUx7HFjw6Z58nt1eT3y2Hyuj7wfAYwThxXwaQDPYClx5DliHKJSH/gNWC0MabYZXqTX1Mf5urs8nI0sN758SzgHGe+dsA5/PI3Va/mcmbrjeMEY5bLNG9uL3fMAK51Xu0yFCh37rR4Znt562zvsf4BLsZxHKka2AXMck7vAsx0We58YBOOn7B/cZmegeM/XC7wIRDuoVztgTlADjAbiHdOzwSmuiyXhuOnru2Qz58LrMZRTP8BYnyVCzjZue6fnH/f4A/bC7gaqAV+dPlzkje2V1PfLzgO4Yx2fhzh/PfnOrdHhsvn/sX5eRuB8zz8/d5SrtnO/weN22dGS19TH+X6G7DWuf55QG+Xz73euR1zgd/6Mpfz9cPAk4d8nre31zQcV2nV4uivG4CbgZud8wV42Zl7NS5X8Hlie+nQf6WUChKBfshFKaWUkxa6UkoFCS10pZQKElroSikVJLTQlVIqSGihK6VUkNBCV0qpIPH/s/Ss5hTLXcMAAAAASUVORK5CYII=\n", + "text/plain": [ + "
    " + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "%matplotlib inline\n", + "\n", + "from numpy import asarray\n", + "from numpy import arange\n", + "from numpy.random import rand\n", + "from numpy.random import seed\n", + "from matplotlib import pyplot\n", + " \n", + "# objective function\n", + "def objective(x):\n", + "\treturn x**2.0\n", + " \n", + "# derivative of objective function\n", + "def derivative(x):\n", + "\treturn x * 2.0\n", + " \n", + "# gradient descent algorithm\n", + "def gradient_descent(objective, derivative, bounds, n_iter, step_size):\n", + "\t# track all solutions\n", + "\tsolutions, scores = list(), list()\n", + "\t# generate an initial point\n", + "\tsolution = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0])\n", + "\t# run the gradient descent\n", + "\tfor i in range(n_iter):\n", + "\t\t# calculate gradient\n", + "\t\tgradient = derivative(solution)\n", + "\t\t# take a step\n", + "\t\tsolution = solution - step_size * gradient\n", + "\t\t# evaluate candidate point\n", + "\t\tsolution_eval = objective(solution)\n", + "\t\t# store solution\n", + "\t\tsolutions.append(solution)\n", + "\t\tscores.append(solution_eval)\n", + "\t\t# report progress\n", + "\t\tprint('>%d f(%s) = %.5f' % (i, solution, solution_eval))\n", + "\treturn [solutions, scores]\n", + " \n", + "# seed the pseudo random number generator\n", + "seed(4)\n", + "# define range for input\n", + "bounds = asarray([[-1.0, 1.0]])\n", + "# define the total iterations\n", + "n_iter = 30\n", + "# define the step size\n", + "step_size = 0.1\n", + "# perform the gradient descent search\n", + "solutions, scores = gradient_descent(objective, derivative, bounds, n_iter, step_size)\n", + "# sample input range uniformly at 0.1 increments\n", + "inputs = arange(bounds[0,0], bounds[0,1]+0.1, 0.1)\n", + "# compute targets\n", + "results = objective(inputs)\n", + "# create a line plot of input vs result\n", + "pyplot.plot(inputs, results)\n", + "# plot the solutions found\n", + "pyplot.plot(solutions, scores, '.-', color='red')\n", + "# show the plot\n", + "pyplot.show()" + ] + }, + { + "cell_type": "markdown", + "id": "d27edab5", + "metadata": {}, + "source": [ + "## Same code but now with momentum gradient descent" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "4c2769f0", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + ">0 f([0.74724774]) = 0.55838\n", + ">1 f([0.54175461]) = 0.29350\n", + ">2 f([0.37175575]) = 0.13820\n", + ">3 f([0.24640494]) = 0.06072\n", + ">4 f([0.15951871]) = 0.02545\n", + ">5 f([0.1015491]) = 0.01031\n", + ">6 f([0.0638484]) = 0.00408\n", + ">7 f([0.03976851]) = 0.00158\n", + ">8 f([0.02459084]) = 0.00060\n", + ">9 f([0.01511937]) = 0.00023\n", + ">10 f([0.00925406]) = 0.00009\n", + ">11 f([0.00564365]) = 0.00003\n", + ">12 f([0.0034318]) = 0.00001\n", + ">13 f([0.00208188]) = 0.00000\n", + ">14 f([0.00126053]) = 0.00000\n", + ">15 f([0.00076202]) = 0.00000\n", + ">16 f([0.00046006]) = 0.00000\n", + ">17 f([0.00027746]) = 0.00000\n", + ">18 f([0.00016719]) = 0.00000\n", + ">19 f([0.00010067]) = 0.00000\n", + ">20 f([6.05804744e-05]) = 0.00000\n", + ">21 f([3.64373635e-05]) = 0.00000\n", + ">22 f([2.19069576e-05]) = 0.00000\n", + ">23 f([1.31664443e-05]) = 0.00000\n", + ">24 f([7.91100141e-06]) = 0.00000\n", + ">25 f([4.75216828e-06]) = 0.00000\n", + ">26 f([2.85408468e-06]) = 0.00000\n", + ">27 f([1.71384267e-06]) = 0.00000\n", + ">28 f([1.02900153e-06]) = 0.00000\n", + ">29 f([6.17748881e-07]) = 0.00000\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtxklEQVR4nO3dd3RUdf7/8ed7Jo0kBEihpRKq9BKaCqKggrhgoYPKVwTFsutaYXWVlbXvKrq6q4h9KdZdUXGRIgJKC106BEgCQgKEBALpn98fM/gbMSEBZnJnJu/HORwz996Z++ImvnK55XPFGINSSinfZ7M6gFJKKffQQldKKT+hha6UUn5CC10ppfyEFrpSSvmJAKtWHB0dbZKSkqxavVJK+aS1a9ceMcbElDfPskJPSkoiNTXVqtUrpZRPEpH9Fc3TQy5KKeUntNCVUspPaKErpZSf0EJXSik/oYWulFJ+otJCF5F3RCRLRH6qYL6IyKsisltENolIZ/fHVEopVZmq7KG/B/Q/x/wBQHPnnwnAvy4+llJKqfNVaaEbY5YCx86xyGDgA+OwEqgrIo3cFfBsGzOO89w329Fhf5VSvqaszPD011vZejDPI5/vjmPosUCGy+tM57TfEJEJIpIqIqnZ2dkXtLJNmcd54/s9bMrMvaD3K6WUVVakHeWtZXvZcdh7C73KjDHTjTEpxpiUmJhy71yt1OBOsdQKtDNrVbqb0ymllGfNWpVO3dBABrT1zEEMdxT6ASDe5XWcc5pHRIQEMrhjY+ZuPEju6WJPrUYppdwq60QB87ccYkjnOEIC7R5ZhzsKfS5wq/Nqlx5ArjHmZzd8boVGd0/kdHEp/13vsd8bSinlVp+kZlJSZhjZPcFj66jKZYuzgRVASxHJFJFxInKXiNzlXGQekAbsBt4C7vZYWqd2cXVoH1eHmav268lRpZTXKy0zzFqVzqVNo2gaE+6x9VQ62qIxZmQl8w1wj9sSVdGobglM+nwza/fnkJIUWd2rV0qpKlu6K5sDx08z+bpWHl2Pz94p+rsOjakdHMBMPTmqlPJyM1emEx0exDWtG3p0PT5b6GHBAdzYOZavN/9MTn6R1XGUUqpcB4+fZvH2wwxLiScowLOV67OFDjCqewJFJWV8ti7T6ihKKVWuj9ZkYICR3Tx3MvQMny70Vg0j6JJYj5mr0vXkqFLK65SUljFnTTq9m8cQHxnq8fX5dKEDjO6ewN4j+azYc9TqKEop9SuLtmdxOK+Q0R68VNGVzxf6de0aUTc0UE+OKqW8zsxV6TSMCOGqVvWrZX0+X+ghgXZu7hzH/C2HyD5RaHUcpZQCIP3oKZbtymZ413gC7NVTtT5f6OA4OVpSZvg4NaPyhZVSqhrMXpOOACO6xVe6rLv4RaE3jQmnZ3IUs1enU1amJ0eVUtYqKinjk9QM+l7SgEZ1alXbev2i0MGxl56Zc5qluy5sWF6llHKXb7ce4sjJIkZV08nQM/ym0K9t05CosCA9OaqUstzMlenE1atF7+YXNkz4hfKbQg8KsDGsazyLth3m59zTVsdRStVQu7NOsiLtKCO7JWC3SbWu228KHWBk1wQMMGe1nhxVSllj9up0AmzCsJTqOxl6hl8VekJUKL2ax/DRmgxKSsusjqOUqmEKikv5bF0m17ZpSEzt4Gpfv18VOjjuHD2UV8Di7VlWR1FK1TDzNv/M8VPF1XZn6Nn8rtD7tqpPg4hgPTmqlKp2M1elkxwdRs+mUZas3+8KPcBuY0TXBJbuyibj2Cmr4yilaojth/JYuz+HUd0TEKnek6Fn+F2hg+POLMFxckIpparDrFXpBAXYuLlznGUZ/LLQG9WpxVWtGvBxagZFJXpyVCnlWfmFJXy+7gAD2zWiXliQZTn8stABRvdI4MjJIr7desjqKEopP/flxoOcLCyx7GToGX5b6L2bxxBbtxaz9OSoUsrDZq1Op0WDcLok1rM0h98Wut0mjOqewI97jrIn+6TVcZRSfmpT5nE2ZeYyunuiZSdDz/DbQgcYmhJHgE2YrXvpSikPmbUqnVqBdm7sHGt1FP8u9Pq1Q7i2TUM+XZdJQXGp1XGUUn4mr6CYuRsPMqhDYyJCAq2O49+FDo5hdY+fKuabn362OopSys98sf4Ap4pKq32Y3Ir4faH3TI6iSXQYM1fqYRellPsYY5i5Kp22sRG0j6tjdRygBhS6zSaM6pZA6v4cdhw6YXUcpZSfWJeew/ZDJ7ziZOgZfl/oADd3iSPIbmPWqv1WR1FK+YmZq9IJDw5gUIfGVkf5RY0o9MiwIK5r15DP1x3gVFGJ1XGUUj7u+Kkivtr0Mzd0akxYcIDVcX5RIwodYHSPRE4UlvDlxoNWR1FK+bhP12ZSVFLGqG6JVkf5lRpT6CmJ9WjRIFyH1VVKXRRjDLNWp9M5oS6tG0dYHedXqlToItJfRHaIyG4RmVTO/AQR+U5E1ovIJhG5zv1RL46I4+TopsxcNmUetzqOUspHrdhzlLTsfEZ19669c6hCoYuIHXgdGAC0BkaKSOuzFnsc+NgY0wkYAfzT3UHd4aYucYQHBzBj2V6royilfNSM5XuJCgvi+vaNrI7yG1XZQ+8G7DbGpBljioA5wOCzljHAmX971AG88kB1REggI7vF8/Xmn8nM0YdfKKXOz67DJ1i8PYtbeyYREmi3Os5vVKXQY4EMl9eZzmmupgBjRCQTmAfcV94HicgEEUkVkdTs7OwLiHvx/u+yJgjw7g/7LFm/Usp3zVi2l+AAG7f09L7DLeC+k6IjgfeMMXHAdcCHIvKbzzbGTDfGpBhjUmJiYty06vPTuG4trm/fiDmr08k9XWxJBqWU78k6UcB/1h9gaEockRY+xOJcqlLoB4B4l9dxzmmuxgEfAxhjVgAhQLQ7AnrCHb2SyS8q1UfUKaWq7IMf91NcVsa4y5OtjlKhqhT6GqC5iDQRkSAcJz3nnrVMOtAXQEQuwVHo1hxTqYK2sXW4rFkU7/6wVx9Rp5Sq1KmiEj5cuZ9rWjegSXSY1XEqVGmhG2NKgHuB+cA2HFezbBGRp0RkkHOxB4HxIrIRmA2MNcYYT4V2h/G9kjmcV6g3GimlKvVJaia5p4uZ0Nt7984BqnTPqjFmHo6Tna7TnnD5eitwmXujedYVLWJo2aA2by1L46bOsV4zuI5SyruUlhlmLE+jc0JduiRGWh3nnGrMnaJnExHu6NWE7YdOsGzXEavjKKW81Pwth8g4dtrr986hBhc6wKCOjalfO5i3lqVZHUUp5YWMMby5NI3EqFCubt3Q6jiVqtGFHhxgZ+xlSSzbdYStB/OsjqOU8jKp+3PYmHGcOy5vgt3m/Ydla3ShA4zulkhokJ0Zy3UvXSn1a28tTaNeaCBDusRXvrAXqPGFXic0kOFd45m74SA/5562Oo5SykukZZ9kwbbD3NIjkVpB3nebf3lqfKED3H5ZE8qM4b0f91kdRSnlJd5evpdAu41beiZZHaXKtNCB+MhQrmvXiFkr0zlRoMMBKFXTHT1ZyKdrM7m5cywxtYOtjlNlWuhOE3onc6KwhI/WZFS+sFLKr324cj+FJd59m395tNCd2sfVpVuTSN79YR/FpTocgFI1VUFxKR+s2E/fVvVpVj/c6jjnRQvdxYReyRw4fpp5m3+2OopSyiKfrcvkWH4R433gRqKzaaG7uKpVfZJjwnhrWRpePhSNUsoDysoMM5btpX1cHbo38e7b/Mujhe7CZhPG90rmpwN5rEg7anUcpVQ1W7jtMHuP5DO+V7JPju+khX6WGzvFEh0exFtL9UYjpWqat5alEVu3FgPaev9t/uXRQj9LSKCdW3sm8d2ObHYePmF1HKVUNVmXnsOafTmMu7wJAXbfrEbfTO1hY3okEhJoY4YO2qVUjTFjWRoRIQEM6+obt/mXRwu9HJFhQQztEs9/1x8kK6/A6jhKKQ/bfzSf//10iNE9EgkPrtJjIrySFnoFxl3ehOKyMt5fsc/qKEopD3tn+V7sNmHspUlWR7koWugVSIoO49rWDfn3ynTyC0usjqOU8pCc/CI+Ts1kcMdYGkSEWB3nomihn8P43snkni7mk1QdDkApfzVz1X5OF5cyvpfv3Uh0Ni30c+iSWI8uifV4+4e9lOhwAEr5nYLiUt77cb/jGcMNa1sd56JpoVdifK9kMo6dZv6Ww1ZHUUq52RcbDnDkZKFPPC+0KrTQK3F16wYkRYUyfekeHQ5AKT9SVmZ4a9leWjeK4NKmUVbHcQst9ErYbcK4XslszMxlzb4cq+Mopdxkyc4sdmedZEJv37zNvzxa6FUwpHMc9UIDma7DASjlN6YvTaNRnRAGtm9kdRS30UKvglpBdm7pmcSi7YfZk33S6jhKqYu0OTOXlWnHuP2yJgT66G3+5fGfv4mH3dozkSC7jX9+t8fqKEqpi/T6d7upHRzA8G6+e5t/ebTQqyg6PJgxPRL5z/pM3UtXyof9dCCX/205xO2XNyEiJNDqOG6lhX4eJvZpSnCAnVcW7rI6ilLqAr20YCd1agUyrlcTq6O4nRb6eYgOD2bsZUl8uekgOw7p0LpK+Zp16Tks3p7FhN7Jfrd3Dlro521Cr2TCggJ4ecFOq6Mopc7TS9/uJCosyOcH4apIlQpdRPqLyA4R2S0ikypYZpiIbBWRLSIyy70xvUe9sCDGXd6E/205xE8Hcq2Oo5SqolVpR1m++wgT+zQlzIeHyD2XSgtdROzA68AAoDUwUkRan7VMc2AycJkxpg1wv/ujeo9xvZpQp1ag7qUr5SOMMfz9253Ur+24uMFfVWUPvRuw2xiTZowpAuYAg89aZjzwujEmB8AYk+XemN4lIiSQCb2TWbQ9i3XpeveoUt5u+e4jrN53jHuubEZIoN3qOB5TlUKPBVzHj810TnPVAmghIj+IyEoR6V/eB4nIBBFJFZHU7OzsC0vsJcZemkRkWJDupSvl5c7snTeuE8IIP7vu/GzuOikaADQH+gAjgbdEpO7ZCxljphtjUowxKTExMW5atTXCggOYeEVTlu06wqq0o1bHUUpVYPH2LDZkHOe+vs0JDvDfvXOoWqEfAFx/rcU5p7nKBOYaY4qNMXuBnTgK3q+N6ZFI/drB/H3BTh2JUSkvZIzhpQU7SYgMZUiXOKvjeFxVCn0N0FxEmohIEDACmHvWMv/FsXeOiETjOATj9yNZ1Qqyc8+VzVi99xg/7Na9dKW8zfwth9hyMI8/9G3uV2O2VKTSv6ExpgS4F5gPbAM+NsZsEZGnRGSQc7H5wFER2Qp8BzxsjKkRDTeiWzyN64Twt2936F66Ul6ktMyxd54cE8YNnc4+7eefqvQryxgzzxjTwhjT1BjztHPaE8aYuc6vjTHmAWNMa2NMO2PMHE+G9ibBAXbu69ucDRnH+W6HX1/co5RP+WrTQXYePskf+7XAbvOP8c4r4///BqkGQ7rEkRAZyt+/1WPpSnmDktIypi3cRauGtRnYzn/GO6+MFrobBNpt/L5vc7YczGP+lkNWx1Gqxls57V1umfMSU6LzsNWQvXPQQnebGzo2JjkmjJcW7KS0TPfSlbJK8bIf6PnwBP5v7Zd0nzAMVqywOlK10UJ3kwC7jfv7tWDn4ZN8temg1XGUqrH2v/gqdlOGAFJUBEuWWB2p2mihu9H17RrRskFtXlm4i5LSMqvjKFXjFOSfJnTZ9xjA2O0QFAR9+lgdq9poobuRzSb88eoWpB3J578bdC9dqeq2+ZGpND5+mH2PPIlMnQqLFkHPnlbHqjb+OYakha5t04C2sRG8smgngzs2rhE3MyjlDU7vy6D1jGmsbXcpnZ97EqTmnAw9Q9vGzUSEB69uScax03ySmml1HKVqjMzx9xJQWkLAq68iNbDMQQvdI/q0jKFTQl3+sXgXBcWlVsdRyu+dWrCY5gvnMm/ALXTo08XqOJbRQveAM3vpP+cWMGd1utVxlPJvxcWcvnMimRH1Sf77VKvTWEoL3UMuaxZF9yaRvL5kD6eLdC9dKU85Pe1Vovbu5L+3PUSHFo2tjmMpLXQPEREevKYl2ScK+XDlPqvjKOWfDh3CNmUKS5p04cpHx1udxnJa6B7UrUkkvZpH88b3aZwsLLE6jlJ+p/CBh5DCQhbf/RhtYutaHcdyWuge9sDVLTiWX8T7P+6zOopS/mX5coJnz+Stbjcy5pZ+VqfxClroHtYpoR59W9Xnze/3kHu62Oo4SvmHkhJKJt7NwYgY9k34Ay0a1LY6kVfQQq8Gf7y6BXkFJcxY5vcPcVKqevzrXwT8tJmn+97BxIHtrU7jNbTQq0Hb2Dpc374R05emkZlzyuo4Svm2w4cpffzPLE/qSO1Rw0mOCbc6kdfQQq8mk6+7BBF4+uttVkdRyqeZSZMw+fm8MPAeHurfyuo4XkULvZrE1q3FvVc245ufDrFsV7bVcZTyTStWIO+9x4yUwdw0si/R4cFWJ/IqWujV6I5eySRGhTJl7haKSnR4XaXOS2kpZffcS1ZENP+78Q7G9Ei0OpHX0UKvRiGBdp64vjV7svP1Mkalztf06djWr+OpPrczaWhXAnQk09/QLVLN+l7SgKta1Wfawp1k5RVYHUcp33DkCKV/eowVSR2wDR9Oj+QoqxN5JS10CzxxfWuKSw3PfbPd6ihK+YbJkzEn8nim/0T+NLC11Wm8lha6BZKiwxjfuwmfrz/Amn3HrI6jlHdbvRrz9tu803kQ1w3vR8M6IVYn8lpa6Ba558pmNKoTwpNfbKG0zFgdRynvVFpK2d33cLR2JHMHjeP2y5OsTuTVtNAtEhoUwGMDL2Hrz3nM0jHTlSrf229jW5vK1N5jeWhoV4ID7FYn8mpa6BYa2K4RPZOj+Nv8HRzLL7I6jlLe5ehRyiZNZk1CW04NGU6flvWtTuT1tNAtJCL8ZXAbThaW8Ldvd1gdRynv8vjjmNxc/nLtRJ74XRur0/gELXSLtWhQm9t6JjF7dTqbM3OtjqOUd1i7FvPmm7zfaSBXDe1HfGSo1Yl8gha6F7j/6uZEhQXxxNyfKNMTpKqmKyvD3H0POeH1mHP9OCZe0dTqRD5DC90LRIQE8mj/VqxPP87n6w9YHUcpa733HrJ6FX/tfRsPDOlKrSA9EVpVVSp0EekvIjtEZLeITDrHcjeLiBGRFPdFrBlu7hxHp4S6PPfNNvIK9EEYqobKyaHskUdZF9+a7BuGcW2bhlYn8imVFrqI2IHXgQFAa2CkiPzmVi0RqQ38AVjl7pA1gc0mPDWoLUfzi3hl4S6r4yhljT//GY4d44mrJ/LkoLaIiNWJfEpV9tC7AbuNMWnGmCJgDjC4nOWmAs8DOkDJBWoXV4cRXRN478d97Dx8wuo4SlWvDRsw//oXH3S6jktv7kuz+vrgivNVlUKPBTJcXmc6p/1CRDoD8caYr8/1QSIyQURSRSQ1O1vHBC/Pw9e2JDw4gCe/2IIxeoJU1RBlZZh77iE3NIIPBozjvquaWZ3IJ130SVERsQEvAQ9WtqwxZroxJsUYkxITE3Oxq/ZLkWFBPHRtS1akHWXe5kNWx1Gqevz738iPP/J0r9v4/c1dqR0SaHUin1SVQj8AxLu8jnNOO6M20BZYIiL7gB7AXD0xeuFGdUugdaMI/vr1Vk4VlVgdRynPys2l7OFH2BR3CfuvH8Lgjo2tTuSzqlLoa4DmItJERIKAEcDcMzONMbnGmGhjTJIxJglYCQwyxqR6JHENYLcJTw1uw8+5Bfzzuz1Wx1HKs558EsnO4vF+d/LkDe30ROhFqLTQjTElwL3AfGAb8LExZouIPCUigzwdsKZKSYrkxk6xTF+axr4j+VbHUcozNm/GvPYaszr2p8PgvrRpXMfqRD6tSsfQjTHzjDEtjDFNjTFPO6c9YYyZW86yfXTv3D0mD2hFoF146qutVkdRyv2MwdxzDydCwpl+ze08eE0LqxP5PL1T1IvVjwjh/n4tWLw9i0XbDlsdRyn3mjULWbaMZy6/hbtu6kbd0CCrE/k8LXQvd9ulSTSNCWPKl1vIL9QTpMpP5OVR9uBDbIltydYBQxmWEl/5e1SltNC9XFCAjWdvak9mzmn++rUeelH+wUyZAlmHebzfXTwzpAN2m54IdQctdB/QrUkkd/ZuyuzVGSzYqodelI/bsgXz6qvMaX8N/W67nraxeiLUXbTQfcQDV7egdaMIJn22iewThVbHUerCGEPBxLvJCwplwejfc5cOjetWWug+IijAxrQRHTlRWMKjn23SYQGUTyqb8xEhy5byypW38dS4Pnqoxc200H1Iiwa1mTygFYu3Z+mDpZXvOXmSU7+/n80NmtL6iQf1KUQeoIXuY27rmUSv5tH89attpGWftDqOUlV25NHHCT9ymHl3Ps6QbolWx/FLWug+xmYTXhzSgaAAG3/8eCPFpWVWR1KqUoWbt1D3jdf4onN/xj8ySm/v9xAtdB/UsE4Iz9zYjo0Zx3lt8W6r4yh1bsZw4Nbx5AeGEP2PvxEZpjcQeYoWuo8a2L4RN3WK5bXvdrMuPcfqOEpVaNtr75K8YQXLbv09l13axuo4fk2suloiJSXFpKbqkC8XI6+gmAHTlhFgF+b9vhdhwQFWR1LqV45n51DQvCUnwyKI3bWFWqHBVkfyeSKy1hhT7vDkuofuwyJCAnlpWAfSj53Su0iV1zHGsOqOB2mYm4289rqWeTXQQvdx3ZOj9C5S5ZXWPP8Gfb98n33detP0xmutjlMjaKH7Ab2LVHmb7G8W0uVP92A3ZSRuWg0rVlgdqUbQQvcDrneRTtK7SJXFysoMx+79I3ZjEECKi2HJEqtj1Qha6H6iRYPaTOrfikXbs5i9OsPqOKoG++GhqbRM+4kyux3sdggKgj59rI5VI2ih+5GxlyZxebNopn61Ve8iVZbYP+e/9Jz2FzZ1vBz57juYOhUWLYKePa2OViPoZYt+5lBuAddOW0pSdBif3tWTQLv+zlbVo3DTTxR378HBug2IXr+ayIZRVkfyS3rZYg3SsE4IT9/YVu8iVdXryBHyr72O0/ZAjsz6RMvcIlrofuj69o1/uYt0vd5FqjytsJDcAb8jLPsQ/3nydS69srPViWosLXQ/NWVwGxpGhPDHjzbos0iV5xhD0bjx1EldyQsjJ3HLH4ZbnahG00L3U2fuIt1/7BR//Xqb1XGUnzLPP0/QzA+Z1ms0Nzz7ILWC7FZHqtG00P1Y9+QoJvROZvbqdD5J1UsZlZt9/jkyeTJfXHIFAU8+Qbs4fTao1bTQ/dxD17Tk8mbR/Ok/m1mx56jVcZS/WLuW0tFjWNe4Jf974GkmXtnc6kQKLXS/F2i38froziRFhXHXv9eyR69PVxfrwAGKr/8dh4Nr8/LdL/D3W7vrs0G9hBZ6DVCnViDvjO1KgE24/b01HMsvsjqS8lX5+RQP/B1Fx47zyG1/5W/3XUNokA7b7C200GuI+MhQ3rothUO5BUz4IJWC4lKrIylfU1ZGyZhbsG3ayEM3TeKxScNpEBFidSrlQgu9BumcUI+XhnUkdX8Oj+ogXuo8lT32GAH//Q/PXDWOYU/exSWNIqyOpM6ihV7DDGzfiIevbckXGw7y8sJdVsdRvuL997E99xwzO/Yn8anJXNmqvtWJVDmqVOgi0l9EdojIbhGZVM78B0Rkq4hsEpFFIpLo/qjKXe7u05RhKXG8umgXn6/LtDqO8nbLllF6x3iWJ3Yg7YlnufXSJlYnUhWotNBFxA68DgwAWgMjRaT1WYutB1KMMe2BT4EX3B1UuY+I8Ncb2tEzOYpHP9vEqjS9nFFVYM8eigbdwP6I+nz8yN/50+AOVidS51CVPfRuwG5jTJoxpgiYAwx2XcAY850x5pTz5Uogzr0xlbsFBdh4Y0wXEiJDufPfa9l7JN/qSMrbHD9O4YCBnCos5tmJL/DsuCv08kQvV5VCjwVcbzPMdE6ryDjgm/JmiMgEEUkVkdTs7Oyqp1QeUSc0kHfHdsMmwv+9u5ocvZxRnVFSQuHNQ7Gl7WHy6ClMffgmwoL18kRv59aToiIyBkgBXixvvjFmujEmxRiTEhMT485VqwuUEBXKW7d24WBuAXd+uJbCEr2cUUHxfb8nePFCplx3H/f+ZRwN6+jlib6gKoV+AIh3eR3nnPYrItIPeAwYZIzRJxX7kC6Jkbw4pD2r9x1j0meb9XLGGq7sH/8g8I1/8Wb3m7jq+Udp01jHaPEVVfk31BqguYg0wVHkI4BRrguISCfgTaC/MSbL7SmVxw3uGMv+o6d4acFOkqLC+EM/HZujRvrf/+D++/m2eQ+Cnn+evpc0sDqROg+V7qEbY0qAe4H5wDbgY2PMFhF5SkQGORd7EQgHPhGRDSIy12OJlcfcd1Uzbuocy8sLd/LFht/8I0z5u59+ovjmoWyLTmL1U68wtndTqxOp81SlsxzGmHnAvLOmPeHydT8351IWEBGevakdmTmnefiTTTSuW4uuSZFWx1LVISuL0/0HkmcLYsYj03hxaAoiekWLr9E7RdWvBAfYeXNMF2Lr1WLCB6ns08sZ/V9BAaevH4QcPsTT459h6r0DCNCHi/sk/a6p36gXFsQ7Y7tiQEdn9HfGUDD2dmqtWcVfhjzCpCduJVwvT/RZWuiqXE2iw5h+SwqZx08z7M0VHDx+2upIygNyHptCyEezeaXPrYx68SEa161ldSR1EbTQVYW6NYnkg9u7cTi3gJv/9SO7Dp+wOpJyo/Q336fes0/xVfu+XPH+NH2EnB/QQlfn1CM5io/u7ElJmWHIGytYuz/H6kjKDTZ9voCYe+9kU2IbLvn6Izom1LM6knIDLXRVqdaNI/h84qXUCw1k9IyVLN5+2OpI6iKkvvIOySMGcapWGA0XzqNpXJTVkZSbaKGrKomPDOXTiZfSvH5txn+wls/W6rC7vmj5lGl0vn8cYcUFRBblUz9b7zfwJ1roqsqiw4OZPaEHPZIjefCTjUxfusfqSKqKjDF89edX6Tb1YQQcf0pKYMkSi5Mpd9JCV+clPDiAd8Z2ZWD7RjwzbzvPzNtGWZmO/eLNSotL+H7oBK7/6x842jgRQkLAboegIOjTx+p4yo30glN13oID7PxjRCeiw4KYvjSNIycLef7m9gTqzShepyD7KLuuHkyfjT+wsf8Q2v/nQ2T9eseeeZ8+0LOn1RGVG2mhqwtiswlTBrUhOjyYvy/YSU5+Ea+P7kxokP5IeYuT6zdyov/vaHXkID88OJXLXnwMRBwlrkXul3SXSl0wEeG+vs159qZ2fL8zm9EzVulDMrzE8dkfY+vZk8CTJ1j51sdc9rfHHWWu/JoWurpoI7sl8M/RXdhyMI+heleptcrKOP7wZOqOGs7uqHh2fbOEXrffZHUqVU200JVb9G/bUO8qtVpuLnn9B1L3b8/xRcerMUu+p2dvfahzTaKFrtzG9a7SoW/qXaXVascOTnVOIXTht7z0u/toO/8zOjRvaHUqVc200JVbnbmrtG4tx12lX248aHUkv2fmzqW4S1dOH8pm0sSXGfXhCzStX9vqWMoCWujK7c7cVdqqYQT3zV7PXR+uJSuvwOpY/qesjLzHnkAGD2Zb7Qb8+fF3+fOLE/WBzjWYFrryiOjwYD69qyeTBrTiux1Z9Hvpez5OzdAHULtJWW4e6VcOIOKZqXzRri8bZ87ltUdvoE5ooNXRlIW00JXHBNht3HVFU775Qy9aNYzgkU83ccvbq8k4dsrqaD4tfdUGDlzSkcbLFvLhiD/S+bsvuOWqS7DZ9LLEmk4LXXlcckw4cyb0YOoNbVmfnsM1Ly/l3R/2UqpDBpyXktIyvn7ubepecTnhx4+y/PWZjJn1d+KjwqyOpryEFrqqFjabcEuPRL594Aq6J0fyly+3MvSNH9mdpZc3VkXa3G/Z2LYHAybfwfEGsZSuXk2fiSP0Qc7qV7TQVbWKrVuLd8d25eXhHdh7JJ/rXlnOa4t3UVxaZnU0r1RQXMpXj75A4g396bJ9DdhsJLz3BtFtW1odTXkhLXRV7USEGzvFseCBK7imTQP+9u1OfveP5WzOzLU6mlfZvHw9S3sO4LoXHsXmPJlsE4GVKy1OpryVFrqyTHR4MK+N6sz0W7pwLL+IG/75A899s52C4lKro1nq1L4MVg8cRcsrunHFxu/Jvv5mRIe8VVWgha4sd02bhix44AqGdonjje/3MOCVZaxKO2p1rOqXk0P6Xfdja9Gczt98xKarb6R0x04afPkpLF4MU6fCokU6UqKqkFh1XXBKSopJTU21ZN3Ke/2w+wiTPt9ExrHT3NgpljE9EuicUM+vT/6ZkyfJeOoFol+fRuipEyzqdBX1X3qedn1SrI6mvJCIrDXGlPvDoYWuvM6pohKmLdzFzJX7yS8qpVn9cIanxHNT51iiwoOtjuc2WUfy2D71Jdq8+w+iThzj++bd2P/HPzHsjusJCbRbHU95KS105ZPyC0v4etPPzFmTzrr04wTahatbN2BYSjy9msdg98EbaUpKy1g/+yt49VXit62j4ckctjXrQPbkJ+g6ZjC1grTI1blpoSuft+vwCT5ak8Hn6w9wLL+IxnVCGJoSz9CUOOLqhVodr1L7j+bz7RfLiP7XPxic+g02DEaEo1OeJvrPk/ThE6rKtNCV3ygsKWXh1iw+Ss1g2a5sAC5vFs2Irgn0a12f4ADv2cMtKC5l2fxVHHnn37T/YT5tstIAMICA46qVqVNh8mQrYyofo4Wu/FJmzik+Sc3kk9QMDuYWEBkWxE2dYhneNZ7mDSwaPnbFCrL//RFp+7II27CWtgd3AnDokg6EjhlFRJuWMHIkFBU5LkHUq1bUebroQheR/sArgB2YYYx57qz5wcAHQBfgKDDcGLPvXJ+pha7cpbTMsHz3ET5ak86CrYcpLjW0alibpvXDia8XSkJkKPGRtYivF0rjurUICrj4q3WLl/3AifkL2d8mhQwJwbZmNbFLF9Bh7ZJfbgI61CiRgjG3kXDXWGzJTf7/m1esgCVLHNeTa5mr83RRhS4idmAncDWQCawBRhpjtrosczfQ3hhzl4iMAG40xgw/1+dqoStPOHqykM/XHWDprmwyjp0iM+c0JS6DgNkEGtWpRVy9WsRH/rrs4yNDiYmrj+3ECcqAsrAw8i69gpza9djY9Uryf84mbvUyivNPcdVPS7Ebx3AFZ45+F9kDCSwtRgBjtyN6OEV5wMUWek9gijHmWufryQDGmGddlpnvXGaFiAQAh4AYc44P10JX1aG0zHAor4CMY6dIP3aKzGOnyMg5/cvrrBOFvyy74aVh1Ckuf2jfs09ZnjkObkQouH4Qwc88jS33OFx9tR5OUR51rkIPqML7Y4EMl9eZQPeKljHGlIhILhAFHDkryARgAkBCQkKVwit1Mew2IbZuLWLr1qJHctRv5hcUl5KZc5qMnFNEPO8o8/KuNzE4xqDBGLDZEJsNjEGCgqg1+VFo28ax4KJFejhFWaYqhe42xpjpwHRw7KFX57qVKk9IoJ1m9cNpVj8cateGE+UP5yuBgY5LC0tLHXvf06bB0aO/Le6ePbXIlWWqUugHgHiX13HOaeUtk+k85FIHx8lRpXxHXh5ERPz/Ug8Ph379oGFDuPVWxzTd+1ZerCqFvgZoLiJNcBT3CGDUWcvMBW4DVgBDgMXnOn6ulNfKyzv3fC1y5cUqLXTnMfF7gfk4Llt8xxizRUSeAlKNMXOBt4EPRWQ3cAxH6SullKpGVTqGboyZB8w7a9oTLl8XAEPdG00ppdT50PHQlVLKT2ihK6WUn9BCV0opP6GFrpRSfsKy0RZFJBvYf4Fvj+asu1C9hOY6P5rr/HlrNs11fi4mV6IxJqa8GZYV+sUQkdSKxjKwkuY6P5rr/HlrNs11fjyVSw+5KKWUn9BCV0opP+GrhT7d6gAV0FznR3OdP2/NprnOj0dy+eQxdKWUUr/lq3voSimlzqKFrpRSfsJrC11EhorIFhEpE5EKL+8Rkf4iskNEdovIJJfpTURklXP6RyIS5KZckSKyQER2Of9br5xlrhSRDS5/CkTkBue890Rkr8u8jtWVy7lcqcu657pMt3J7dRSRFc7v9yYRGe4yz63bq6KfF5f5wc6//27n9khymTfZOX2HiFx7MTkuINcDIrLVuX0WiUiiy7xyv6fVlGusiGS7rP8Ol3m3Ob/vu0TktmrO9bJLpp0ictxlnie31zsikiUiP1UwX0TkVWfuTSLS2WXexW8vY4xX/gEuAVoCS4CUCpaxA3uAZCAI2Ai0ds77GBjh/PoNYKKbcr0ATHJ+PQl4vpLlI3EMKRzqfP0eMMQD26tKuYCTFUy3bHsBLYDmzq8bAz8Ddd29vc718+KyzN3AG86vRwAfOb9u7Vw+GGji/Bx7Nea60uVnaOKZXOf6nlZTrrHAa+W8NxJIc/63nvPretWV66zl78Mx7LdHt5fzs3sDnYGfKph/HfANjicd9gBWuXN7ee0eujFmmzFmRyWLdQN2G2PSjDFFwBxgsIgIcBXwqXO594Eb3BRtsPPzqvq5Q4BvjDHlP33Yfc431y+s3l7GmJ3GmF3Orw8CWUC5d8JdpHJ/Xs6R91Ogr3P7DAbmGGMKjTF7gd3Oz6uWXMaY71x+hlbieHKYp1Vle1XkWmCBMeaYMSYHWAD0tyjXSGC2m9Z9TsaYpTh24CoyGPjAOKwE6opII9y0vby20KuovAdYx+J4QPVxY0zJWdPdoYEx5mfn14eABpUsP4Lf/jA97fzn1ssiElzNuUJEJFVEVp45DIQXbS8R6YZjr2uPy2R3ba+Kfl7KXca5Pc488Lwq7/VkLlfjcOzlnVHe97Q6c93s/P58KiJnHlfpFdvLeWiqCbDYZbKntldVVJTdLdurWh8SfTYRWQg0LGfWY8aYL6o7zxnnyuX6whhjRKTC6z6dv3nb4Xja0xmTcRRbEI5rUR8FnqrGXInGmAMikgwsFpHNOErrgrl5e30I3GaMKXNOvuDt5Y9EZAyQAlzhMvk331NjzJ7yP8HtvgRmG2MKReROHP+6uaqa1l0VI4BPjTGlLtOs3F4eZWmhG2P6XeRHVPQA66M4/ikT4NzLKu/B1heUS0QOi0gjY8zPzgLKOsdHDQP+Y4wpdvnsM3urhSLyLvBQdeYyxhxw/jdNRJYAnYDPsHh7iUgE8DWOX+YrXT77grdXOS7mgedVea8ncyEi/XD8krzCGFN4ZnoF31N3FFSluYwxrg+Dn4HjnMmZ9/Y5671L3JCpSrlcjADucZ3gwe1VFRVld8v28vVDLr88wFocV2WMAOYax1mG73AcvwbHA6zdtcd/5oHYVfnc3xy7c5bamePWNwDlng33RC4RqXfmkIWIRAOXAVut3l7O791/cBxb/PSsee7cXuX+vJwjr+sDz+cCI8RxFUwToDmw+iKynFcuEekEvAkMMsZkuUwv93tajbkaubwcBGxzfj0fuMaZrx5wDb/+l6pHczmztcJxgnGFyzRPbq+qmAvc6rzapQeQ69xpcc/28tTZ3ov9A9yI4zhSIXAYmO+c3hiY57LcdcBOHL9hH3OZnozjf7jdwCdAsJtyRQGLgF3AQiDSOT0FmOGyXBKO37q2s96/GNiMo5j+DYRXVy7gUue6Nzr/O84bthcwBigGNrj86eiJ7VXezwuOQziDnF+HOP/+u53bI9nlvY8537cDGODmn/fKci10/n9wZvvMrex7Wk25ngW2ONf/HdDK5b23O7fjbuD/qjOX8/UU4Lmz3ufp7TUbx1VaxTj6axxwF3CXc74Arztzb8blCj53bC+99V8ppfyErx9yUUop5aSFrpRSfkILXSml/IQWulJK+QktdKWU8hNa6Eop5Se00JVSyk/8P2O4xqLPu9cnAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
    " + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "from numpy import asarray\n", + "from numpy import arange\n", + "from numpy.random import rand\n", + "from numpy.random import seed\n", + "from matplotlib import pyplot\n", + " \n", + "# objective function\n", + "def objective(x):\n", + "\treturn x**2.0\n", + " \n", + "# derivative of objective function\n", + "def derivative(x):\n", + "\treturn x * 2.0\n", + " \n", + "# gradient descent algorithm\n", + "def gradient_descent(objective, derivative, bounds, n_iter, step_size, momentum):\n", + "\t# track all solutions\n", + "\tsolutions, scores = list(), list()\n", + "\t# generate an initial point\n", + "\tsolution = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0])\n", + "\t# keep track of the change\n", + "\tchange = 0.0\n", + "\t# run the gradient descent\n", + "\tfor i in range(n_iter):\n", + "\t\t# calculate gradient\n", + "\t\tgradient = derivative(solution)\n", + "\t\t# calculate update\n", + "\t\tnew_change = step_size * gradient + momentum * change\n", + "\t\t# take a step\n", + "\t\tsolution = solution - new_change\n", + "\t\t# save the change\n", + "\t\tchange = new_change\n", + "\t\t# evaluate candidate point\n", + "\t\tsolution_eval = objective(solution)\n", + "\t\t# store solution\n", + "\t\tsolutions.append(solution)\n", + "\t\tscores.append(solution_eval)\n", + "\t\t# report progress\n", + "\t\tprint('>%d f(%s) = %.5f' % (i, solution, solution_eval))\n", + "\treturn [solutions, scores]\n", + " \n", + "# seed the pseudo random number generator\n", + "seed(4)\n", + "# define range for input\n", + "bounds = asarray([[-1.0, 1.0]])\n", + "# define the total iterations\n", + "n_iter = 30\n", + "# define the step size\n", + "step_size = 0.1\n", + "# define momentum\n", + "momentum = 0.3\n", + "# perform the gradient descent search with momentum\n", + "solutions, scores = gradient_descent(objective, derivative, bounds, n_iter, step_size, momentum)\n", + "# sample input range uniformly at 0.1 increments\n", + "inputs = arange(bounds[0,0], bounds[0,1]+0.1, 0.1)\n", + "# compute targets\n", + "results = objective(inputs)\n", + "# create a line plot of input vs result\n", + "pyplot.plot(inputs, results)\n", + "# plot the solutions found\n", + "pyplot.plot(solutions, scores, '.-', color='red')\n", + "# show the plot\n", + "pyplot.show()" + ] + }, + { + "cell_type": "markdown", + "id": "3f3919ce", + "metadata": {}, + "source": [ + "## Using Autograd with OLS\n", + "\n", + "We conclude the part on optmization by showing how we can make codes\n", + "for linear regression and logistic regression using **autograd**. The\n", + "first example shows results with ordinary leats squares." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "558d6dc4", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Own inversion\n", + "[[4.08805934]\n", + " [2.86802426]]\n", + "Eigenvalues of Hessian Matrix:[0.26504701 4.42519896]\n", + "theta from own gd\n", + "[[4.08805934]\n", + " [2.86802426]]\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAq0klEQVR4nO3deXxcdb3/8dcn3dNSuu9NSgWKLMoSQUFs3S6LKC5XhVuxKFqpP67oRblqVUAuyHW5uFwVyyLVRkBxF7yCkKSF0tYUS6Eb3Qule6F0o2mSz++Pc1ImYWZyZj0zyfv5eMwjM2f9nDOT7+d8v9+zmLsjIiLSmYq4AxARkfKghCEiIpEoYYiISCRKGCIiEokShoiIRKKEISIikShhSJdlZteb2Zy448gnM5tiZs/HHYd0T0oYUlRmtsHMDprZPjPbamZ3m9mAuOMSkc4pYUgc3uvuA4BTgdOAr8QbTvdkZj3jjkHKixKGxMbdtwJ/I0gcAJjZl81srZntNbPlZvaBhHGXm9ljZvZdM3vRzNab2QUJ448xs4Zw3oeBYYnrM7P3mdkyM3vJzOrN7PUJ4zaY2ZfMbKmZ7TezO81spJn9NVze381scLLtaGsmMrNrzGy7mW0xs08kjK83s0913I6Ez25mnzWz1eG6bjSz15nZfDN72cx+bWa9O6zzq2a2M4x7asLwPuH+2WRm28zsNjPr1yHO/zSzrcDPzWyYmf0l3Ce7zWyemalckKT0w5DYmNk44AJgTcLgtcC5wNHADcAcMxudMP4sYBVBMvg2cKeZWTjuV8DicNyNwLSEdR0P3AN8HhgOPAj8uUNB/CHg3cDxwHuBvwJfDaevAD6XZnNGhTGPBa4AfpwqwaRwHnAG8GbgWmAW8DFgPHAycGmHdQ0L1zUNmGVmk8Jxt4TxnwocG07zjQ7zDgGqgenANcDz4TaODLdX9wuSpJQwJA5/MLO9wHPAduC6thHu/ht3f8HdW939PmA1cGbCvBvd/XZ3bwFmA6OBkWZWBbwJ+Lq7H3L3ucCfE+b7KPCAuz/s7oeB7wL9gLMTpvmRu29z983APGChu//T3V8Bfk/QfJbKYeCb7n7Y3R8E9gGT0kzf0bfd/WV3XwY8Azzk7uvcfQ9B4uq47rbtbAAeAD4SJs7pwBfcfbe77wVuBi5JmK8VuC6c92AY92igOox9nusGc5KCEobE4f3ufhQwBTiBhKYjM/u4mS0Jm0heIji6Tmxa2tr2xt0PhG8HAGOAF919f8K0GxPej0n87O6tBAlrbMI02xLeH0zyOV3n/C53b074fKCT6TvKZN3JtnMMQS2hElicsP/+LxzeZkeYANt8h6CG95CZrTOzL2cQs3QzShgSm/Do+G6Co33MrBq4HbgKGOrugwiOti3FIhJtAQabWf+EYVUJ718gaIYhXJcRNPdszn4LIttPUJC3GZXj8pJt5wvAToLkcpK7DwpfR4cnGLRpV3tw973ufo27TwTeB/yHmb0zx/iki1LCkLh9H3i3mb0R6E9QoO0ACDuOT46yEHffCDQCN5hZbzN7K0E/RJtfA+8xs3eaWS+CtvtDwPx8bUgaS4APmlmlmR1L0MeRq7btPBe4CPhNWGu6HbjVzEYAmNlYMzsv1ULM7CIzOzZMoHuAFoJmK5HXUMKQWLn7DuAXwDfcfTnwPeAJgiaZU4DHM1jcvxF0iu8m6Bf5RcJ6VhF0Iv+I4Ej8vQSn9zblYTM6cyvQRLBNs4HaHJe3FXiRoFZRC1zp7ivDcf9J0MS0wMxeBv5O+r6U48Jp9hHs95+4e12O8UkXZerfEhGRKFTDEBGRSJQwREQkEiUMERGJRAlDREQiKaubjw0bNswnTJgQdxgiImVl8eLFO919eOdTpldWCWPChAk0NjbGHYaISFkxs42dT9U5NUmJiEgkShgiIhKJEoaIiESihCEiIpEoYYiISCRKGCIiEokShoiIRKKEISIikShhiIhIJEoYIiISScEThpndZWbbzeyZJOOuMTM3s2GFjkNERHJTjBrG3cD5HQea2XjgX4BNRYhBRERyVPCE4e5zCZ6x3NGtwLWAnhErIlIGYunDMLOLgc3u/lSEaaebWaOZNe7YsaMI0YmISDJFTxhmVgl8FfhGlOndfZa717h7zfDhOd/OXUREshRHDeN1wDHAU2a2ARgHPGlmo2KIRUREIir6A5Tc/WlgRNvnMGnUuPvOYsciIiLRFeO02nuAJ4BJZva8mV1R6HWKiEj+FbyG4e6XdjJ+QqFjEBGR3OlKbxERiUQJQ0REIlHCEBGRSJQwREQkEiUMERGJRAlDREQiUcIQEZFIlDBERCQSJQwREYlECUNERCJRwhARkUiUMEREJBIlDBERiUQJQ0REIlHCEBGRSJQwREQkEiUMERGJRAlDREQiUcIQEZFIlDBERCSSgicMM7vLzLab2TMJw75jZivNbKmZ/d7MBhU6DhERyU0xahh3A+d3GPYwcLK7vwF4FvhKEeIQEZEcFDxhuPtcYHeHYQ+5e3P4cQEwrtBxiIhIbkqhD+OTwF9TjTSz6WbWaGaNO3bsKGJYIiKSKNaEYWYzgWagNtU07j7L3WvcvWb48OHFC05ERNrpGdeKzexy4CLgne7uccUhIiLRxJIwzOx84FpgsrsfiCMGERHJTDFOq70HeAKYZGbPm9kVwP8CRwEPm9kSM7ut0HGIiEhuCl7DcPdLkwy+s9DrFRGR/CqFs6RERKQMKGGIiEgkShgiIhKJEoaIiESihCEiIpEoYYiIlKvaWpgwASoqgr+1KW+akRexXektIiI5qK2F6dPhQHjt88aNwWeAqVMLskrVMEREytHMma8mizYHDgTDC0QJQ0TKV5GbZErKpk2ZDc8DJQwRKU9tTTIbN4L7q00y3SVpVFVlNjwPlDBEpDzF0CRTUm66CSor2w+rrAyGF4gShoiUpxiaZErK1KkwaxZUV4NZ8HfWrIJ1eIPOkhKRclVVFTRDJRveXUydWtAE0ZFqGCJSnmJokunulDBEpDzF0CTT3alJSkTKV5GbZLo71TBERCQSJQwREYlECUNEpBSV4FXs6sMQESk1MdxYMIqC1zDM7C4z225mzyQMG2JmD5vZ6vDv4ELHISJSNkr0KvZiNEndDZzfYdiXgUfc/TjgkfCziIhAyV7FXvCE4e5zgd0dBl8MzA7fzwbeX+g4RKSMlWB7fkHFcGPBKOLq9B7p7lvC91uBkakmNLPpZtZoZo07duwoTnQiUjrivittsZJV4nr27YNevdqPL4Gr2GM/S8rdHfA042e5e4271wwfPryIkYlISShWe36yxFCsZNVxPbt2BVevDx2a01XsO1bs5LdfeiJvYcZ1ltQ2Mxvt7lvMbDSwPaY4RKTUFaM9P9VZSf36pU5WuZ6tVFsbLGfTpiBJtbS0H9/UBAMGwM6dkRe5fdkOGu5cQ8Pfm6h/dgzLDh0HDMstzgRxJYw/AdOAW8K/f4wpDhEpdcW4K22qWkzHYW1yTVYdE1THZBFxPdue2UHDnaupf7iZ+tVjWdH0OmA4/dnHOUNXMnVyPZM/OJRzrswt3DbFOK32HuAJYJKZPW9mVxAkineb2WrgXeFnEUmmu3X4dlSMu9JmmgByTVbJElSE9WxZso17PzefK0+cy+v7rGXUKcP56PfP5pVlq6lrfiutVPDKkDG8dNfv+NvOGr7ytymc/ZlTcos1kbuXzeuMM85wkW5lzhz3ykr3oGU7eFVWBsO7kzlz3Kur3c2Cv/ne/urq9vu47TV0aGH2v1ny9XVYz65v/q/XfvYxn35Cgx/fa92RUUexxy8cvsj/+4I6X/3pW7y1kxiBRs9DGRx7EsjkpYQh3U6qgqy6Ou7IupZ0iTmXZJVq3hTfa2tFD2/FfEfPkf75iluPjBrIS37RiIX+nffU+T9mL/PDBw+/uo4IvxElDJHuINWRqFnckZWHTAr7fNdiOktCHcbtp59fyhwH90H2or935AL/3vvqvPGXy735UHPq9UT4jShhiHQHqmFkL9/NeZkmlBTf3d4BI/0Tx831q+0Hvp5qb8F8I+P8uwOv9/+5uM6f/NWK9gmis/WqhqGEIeLu6sPIRT6TbRbfQ2uKI/8WzIfYLn//6Cf8+x+o9yX3rfSWwy3ZrzfCNEoYIt1FoTt8u6p8NudFSD6tLa2+5pENfuflc/2yifP8OcYmnefQsDGpE0QW63X3Tn8j+UoYFiyrPNTU1HhjY2PcYYgItL/wrKoqOM21lB6XOmFC8us3qqthw4bMllVRERTVHbgZd17WQMM8o37TRM5tqeNmZlLFJvbZQCrZR09PuMaisjKzK7ZTrBczaG2NHL6ZLXb3msgzpAon1wWISDcU9/2dosjn9RsprrvY6FV8+hfn8tCGSXzp6FnM7nEFE9hIBc5A30PPXj1yu71Hid2EUAlDRDJXos9raGfq1KCArq7OqsD2Vmflg+v42dS5/Gjv5RygX7vxh6wPa95yGSseWMfW5mF87qi76dVyqP1C2m7v0doa1GoyrYEV46LFDKhJSkQyl6emklLSliDqazdT/3hPGjYfy7bWEQCMrtjKVwf/lGkHf8aAg9th/Hjs5pvbJ4BC7ZM8NP3lq0lKCUMkLqXeB5BOPvsHYuKtzvI/r6W+djMNT/SmYfOxbPfgjthjK7YwpWodk89tYcplVRz7zmqswtIvsIT3Sb4SRuxnPmXy0llS0mWU++myZRh/y+EWX3r/Kv/Rv9b7h8bO9+G2/Ujo43ps9o8dM8/vmDbX1zyywVtbWtvPHOVMtainwMZwxhs6rVakjHWFC/JK/HTflsMt/tRvVvkPPljvHxj9hA+1nUd2c1WP5/wnw7/ue/uP9FYzb62qSh1/Jskx3T6JMcnmK2GoSUokDl2wDyBurc2tLP3taurv3UrDwj7M3Xo8u30IABN6PseUCRuYPNmYcvkEJmxsaH97cUh9ymuqpqahQ5M/qyJVU2OMTVbqwxApZyXc3l0uWppaeOo3z9Lwm23UL+zH3G2TeMkHATCx50YmH7OJKVOMydMmUH3OuPYzZ7L/UyV3gDlz2ieYjs+5gFcT0WWXxXaQoD4MkXJWin0AuTYxFbiJqvlQszf+crl/96I6v2jEQj+al47sumN7rfcrjp/rv7zyMd+0YHPnC8vkKvBUzYfJmhDTNTXG2AyJ+jBEykiywrSU+gByTWAFSICHDx72RXcv829fWOfvGbHQByYkiON6rfNPn9Dgc2Y85s8teiHzhWdSeM+ZkzphdEww6RJRF+jDiD0JZPJSwpCyFGdtImpSyvXoNw9Hz037m3zBHU/7f19Q5xcMX+RHsefIYib1XuvTT2jwX131uG9evCXyMlPK9DsZOjTa9nW2H3SWlBKGSFpxNUVk8lCgqEfQqWRxo7+m/U0+/2dL/Vvn1fl5Q//hA3j5yGyv773Grzyxwe/93OP+wj+35mV3vEamz8qIkmBKsanRlTBEykehH4KU4VPdkj52NFWMeaxhHNp7yB+/banf9O46/1b/b/pGxnsL5uup9i/2+J5/9uR6v+/zj/vWp7fnZ79ko7PTYqMkmFJqagwpYYiUi0LWMNId0UZ5bnS6pJFjH0Zrv0pfMe1m/6931fm7hjR6Jfsc3C9ljh+gb/braltfvgvlEq0dpJTBPugSCQP4ArAMeAa4B+ibbnolDClLhSyIsjkrJ90rh0K46fa7/eDg0d6K+ZaKUT6Nu44s9g19V/q/v6Hef/ulJ7x5zLjME2hi4Th0qHuvXvnfn+V0MWWGv6myTxjAWGA90C/8/Gvg8nTzKGFI2SpUM0U2Z+VE7cDtxMEXD3r99//p10+p87cPetL7ciBYNS3+xr4r/epT6/33X17gO5/dFT3mZJJtRyEK9nJ6fnqGyS1fCaNnzhdy5KYn0M/MDgOVwAsxxyNSGFOnFubGglVVyS9Aq6p6dX0zZwbT9OgRXFBmKW6id+GFaVd1cPdBFty9kvo/7qHnP+bz8YO3cS7PU814DvW6ilNP38Pk8/px7hXHM+R1k4BJmcecTLJbqSezaVPn06STaVxxSrWtue6DzuQj62T7Aq4G9gE7gNoU00wHGoHGqqqqbPOxSNeU7TOfIxyd7t+x3//+7cX+9XPr/NyBS7w3rzi4/xu/yK0PItMmuqh9MbnWMMqpDyOmGkacyWIw8CgwHOgF/AH4WLp51CQlOSnBs1c61VnMc+a0b2IaOvS100Tsy2g184f/e7HPPKfO3zpwiffikIN7Bc1eU7nMv1hT53+5bpG3jB2fe4GdyXcRJf58Fezl8hvphn0YHwbuTPj8ceAn6eZRwpCsldPRY5vOYo66TRGP0DdQ5eDeg8N+Zv+n/doz6/yB6xf5nuf2RFteodr6k21n795Bciz1gr2QutNZUsBZBGdIVQIGzAb+Pd08ShiStWzOgIn7aLOzmKNuU4Qj9AP09XteN9P/euM//OXNL+cWV6J87cO4v4syV/YJI9gGbgBWEpxW+0ugT7rplTAka/k4MyebGkkuBV1nMaerObj7nuf2+IM3LPJfTZz5mj6HQ/T0fT0HeivmLePGZ34NRBlf9dwddYmEkelLCUOylmkNIx/n5OdaYGZZw2gF/1rvW7yCZgf3Xhzy6/rc7Lv7jMwuQaTats4SYTld19DFFS1hAA8Db8zHynJ9KWFI1vJ1Zk4m7fS5FpidxLzv+7d7a4oaxhYb7V97a53//duLff+O/dFjzqdyuq6hiytmwjgdqAN+DozOx0qzfSlhSE7ycWZOJkfH+SgwE2JuHjPOG993vX/h9Do/rd9yN1q8JVWTVKHvUxWFahglo+hNUsCHgKeA6wivzi72SwlDiiYf7e85Fpi71uz23395gV99ar2f2m+FGy0O7n046FMGPenXT6nzg4NHF65QLsFnZEh2ipowwrOYTgauBHYCzwOX5SOATF5KGFJU+XgCXQYF5p7v3Ob7Bo7yVsw3M8YvZY6De18O+DsGL/Yb3l7nDT9c4gdfPJj1OjKSr34cnd0Uu2I2ST1OcMuOh4EbgYuAY4EfAbPyEUTUlxJGN5aq4Cn1AilNfNuX7/D7vzjfrzql3q/t8V3fR/uC/1BFH18x7WZ/Zc8rWa8jJ+qD6DLylTAsWFZqZnYSsNyTTGhmK9z99WkXkEc1NTXe2NhYrNVJqaithenT299PqLISpk2D2bNfO3zWrMLctylH25ftoOHONdQ/3ETD6jEsO3QcAJXsZ13FsYxs3framaqrYcOG4gbaZsKE5PdWijMmyYqZLXb3mlyXU9HZBO6+LFmyCL0n1wAkQW1t8E9aURH8ra2NO6LSkOzmcwcOwE9/mnz4zJmFi6Wz7yhhfPOY8Txx/vV89uS5nNhnLSNPHs5Hbn0Ls585g3ED9nDzv9Qz/2dP8+LeXoz0bcnXV+ibyaVz001BAk5UWRkMl+4pH9WUYr26dJNUd+sgzKQZJZsHARUq5jTf0e6bfuyHe/ZpN34flf4J7vALhi/yW86v8wV3PO1N+5teu+xSPaOo1Jv8JBJ04V4XU6oFRiFkmhwzfRBQofZZijh29hzpx/da5+tJPr51fIS7LOfjjKRyKtjLLd4yp4TR1XSnDsZMk2PU23N3VsjmWEi1pviOWjC/aMRCbyWD7zBZLNnGl2z/mLnPmJHR9hVNd6tNlwAljK6mO9UwskmOiYVpjx7J5+/RI32y6PhYz1690hZSG+c/77M/Pc8/edxcn9hzQ+oaRNtzWlJ9h0OHtk8EM2bkt8BMtd62p+513H9xH9F3p996iVDC6Gq601FXqgKjR49oBVo2+yrVY0mHDj0yyfp5z/ndn5rnlx8714/pufHIJINtt188aoE/eMZMb+6d5sFBqW7D3TFRpUqY2RaY6fp42pJrx2ni/G11p9p0iVDC6IqKeRQY5xFnlCamzgq0TONPsZ5W8Gmvm+fVPZ47MniI7fIPjH7Cf/DBel9y30pvOdySfr2Jw4YObf+chlSJKp8FZqZ9PHEf0auGUXRKGJK9UqjNRGliymMBkuomfa3gw2yHf2jsfP/hh+p96f2r2ieIzsyYkf7oPZMzvLLd3jlzMj+TLM4j+lL4/XUzShiSvVI7witAE0VrS6s/+9B6v/3jc33qhMd8J4OTruPwUYMzSxCJ0hXUnT3kKN9NRMkSV6nWMNxLq0+lG1DC6Cri+McptTbkPCSw1pZWX/V/6/xnUxv80urHfEzFC0cWM8K2+w+GXOfNFR36Enr3zm1/p2sKatuXqY6mZ8zI//cepdamI/puSQmjK4iral5qNYwsTgttbWn1FQ+s9Z9e2uCXVD3uoyq2Hpl1VMVWv6Tqcf/ppQ2+4oG13trS+up68llIpzuiT9yXcRwUpNqnbbEpWXQrShhdQb4K7kwLpFJsQ+6kL6C1pdWX/XG1/+SSBv/I+Md9ZMW2I5ONqXjBL61+zH82tcFfuPbW4DTXYhTOUU5njZOafSSkhNEV5OsBO1EL/3Rn88RdmKQofPdUjvR/HTvfR9j2I4PHVrzgUyc85rd/fK4/+9D69jWIYibCcrtgTrotJYyuIB81jKjLiKtWEfEoN91V1ON7PO+XTZznd0yb62se2fBqgugojqY2HcVLGegSCQMYBNwPrARWAG9JN32XSxj5KMSj1lLSddAWqqBLs30th1t8yX0r/QcfrPcPjH7CNzEuaWxNI8emThAdlVpnvkiJ6CoJYzbwqfB9b2BQuum7XMJwz/0INepRdWenXBaitpEitm0Vo3yw7T4y6JieG/22EV/3wz365BZTqXXmi5SIsk8YwNHAegge4hTl1SUTRq6i1lKiXA2cx4K1+VBzypvxtWD+yePm+uxPz/MNjz3XfluyTZ5z5iS/qrptX3SXq+hFkugKCeNUYBFwN/BP4A6gf5LppgONQGNVVYTbRHdHUQqoKLfjyKHp5vDBw/6P2cv8uxfV+UUjFvrRvJTyZn15S0xt290We8f1DB36arIoVv9NKZ6BJt1eV0gYNUAzcFb4+QfAjenmUQ0jR4kFbI4F+eGDh33hXc/4ty+s8wuHL/KBvHRkMcf3WuefPqHBH3vXN7ylb7/CFJ5REmBnV1sXoqlKzWJSgrpCwhgFbEj4fC7wQLp5lDDyJIuj4Kb9Tb7gjqf9lvPr/ILhi/wo9hyZdVLvtf6Z1zf4r6563Dcv3vLadRWieSZKE1tbjamYneHp4klGzVdSBGWfMIJtYB4wKXx/PfCddNPnPWF053/WTra9aX+Tz//ZUv/WeXV+3tB/eH/2Hin3Xt97jc84qcHv/dzjvuWpbbGEH+m+SXHUMNI9q6MjNV9JkXSVhHFq2D+xFPgDMDjd9HlNGOX4z1rABHdo7yF/7CdP+U3vrvN3D2lslyBO6vOsf/bkev/1F+b71qe3522dOemshtHZcyoK9V1nUsNQ85UUSZdIGJm+8powyu2fNc+F3it7XvF5P37Kb3xnnb9rSKP3Y/+RxZ7c51m/6pR6v/+L83378h153pA8yfReScWqTWbyu9J1I1IkShi5iuOfNZdCK8cE98qeV7zhh0v8hrfX+TsGL/a+HHBwv5Q5vpkx3or5voGjfM93bstiwwqks/1Vik2KmST2cjtokbKlhJGrdP+shSiIcq0hZJjgDr540Otu/adfP6XOpwx68kiCMFr81H4r/OpT633Re6/31kKdxZSrYt4WPN+i/n7KsVlUypISRq7SFUiF+CfO9Wiyk/kP7Drgj37vSf/G2+p88tH/9D4cPJIgTuu33L9wep3/8asLfPe6F/MXUyEV68FDcSvFWpJ0OUoY+ZDsn7VQhWhnz07I4pbkzb36+q+Pn+nnDlzivXnFwb2CZj+jcplfc0ad/+lrC/3FDS9lHlMptKEX49GmIt2EEkZUmR7BFaoQzeSsniTxt5r5ocqjfW+Pgd6C+Qaq/FLmeAXN/qb+z/gXa+r8L9ct8pc2pkkQUWMqhQI4ynUWpZTgREqYEkYU2bQRF6oQzeTKZHffu2WvP/XhG/1QRfsb8u2j0r/W+xa/9sw6f/CGRb7nuT35jalUmnjSnQVViglOpIQpYUSRTeFfyEK0k1tztGL+n2fV+ZsHLPWeNBX+XkyJMZViG3rH2ArVvyTSxeUrYViwrPJQU1PjjY2N0WeoqAiKlY7MoLU19Xy1tTBzJmzaBFVVcNNNMHVq5gGnMmECbNz4msEbqOY4VnPmUSuYfPJubnriHRhZxN+VFfq7EemCzGyxu9fkvJwunTBSFMxUV8OGDfkKK5I9m/bw2F3PUv/Afvo8tYivHL6B/hw4Mr6pog+rPvx1Jv7w8/Qf0T8YWELxi0j5ylfCqMhHMEVVWxsUpBUVwd/a2tTT3nQTVFa2H1ZZGQwvsJc27uHPX1/ENTX11PRfzpDqAVx0w5v4YeNbmNvvPB6c9AUODhqFm0F1Nb1/cSen3Dvz1WSRKn4zuPDCgscvIvIa+WjXKtbrjGOOybwNO5OLqHJoy9+1Zrf/4SsL/Aun1/lp/Za70eLg3oeDPvnof/o33lbnj37vST+w60BGy/UZM7retQciUlR0yz6MPn28sanptSNybaKprYXp0+HAq01EVFbCrFkp28d3rd7NvLtWU/9/r1C/chRLXzkOp4I+vMLZg1Yw+dQ9THn/IM6adgJ9B/XNPjY1S4lIjrpnH4aZJ+3ByLUTOEKhvHPVLubeuZr6vx2i4dlRLH1lEgB9OcjZg1cw5bSXmfz+wZx52aTcEkRH2XbclyN1aIsURL4SRs98BFM0vXtDshpGVVVuy920Kelg37iJz72hgfpnx/DMoeOAofTjAOcMWcGN59Qz5YNDeNPHJtFn4Om5rT+dqqrkySzXbS41HWt5GzcGn0FJQ6RElFcNY+JEb9y2LaOmo0jSnOZ6Est469AVTD59H1M+NJSaqZPoPaB39uvKVBbNZWVJTW8iBdM9z5IaMiQoKKurgyaZ6uqsC86tS7dz39XzmXHSXL70/NXsp/3ZSIcr+tD0ySt5aX9v/razhq8+NIWzP3NKYZJFujO/pk7N2zaXtBS1vJTDRaToyquGkel1GAleeHIrDT9fR8MjzdSvHceqpokADGAv5w5fyYzRf+Rfnr+L3i9uxYrZft5dahCdUQ1DpGC6Z6d3Bgljc+MW6u9aR0NdC/Vrx7P68DEADGQP545YxeSaA0z5yAhO++jx9Oybx66cTDtuVVAGlDhFCqZ7dnqn8dzCF2i4ez31j7bSsH48aw5PAEZzdJggPnPmRiZ/eASnfuR4evY9szBBZNNxq6aYQNv+0VlSIiWrbGsYm57YTP3P11Nf7zSsr2JdczUAg+wl3jZyFVPOOsjkD4/kjR8+nh69exQnwGxqC6phiEiBdZkmKTPrATQCm939onTTHjP8BJ8y6HbqN0xgQ/N4AAbbi0wetYrJZx1iyiWjOOUDxxYvQXSUzTUTaooRkQLrSmdJXQ2siDLhhp0D+PPaEzlt+GZ+8MEGnvrNs+xsOprfv/BmPv/7yZz60UnJk0Um95/KRaprI9JdM5HqLCgoTswiIlHl4/4i2b6AccAjwDuAv3Q2/YkTT/KWwy2Z3USlmA8Jyte6SvnBRiJSdsjTvaTirmF8H7gWSHmPCzObbmaNZta4r3kvFT0zDHnmzPbNPRB8njkz42A7la9rJooZs4hIRLH1YZjZRcCF7v5ZM5sCfNE76cPI6jqMcrwXUznGLCIlqyv0YZwDvM/MNgD3Au8wszl5X0s2/QpxK8eYRaTLiy1huPtX3H2cu08ALgEedfePpZ1p9+7MO4JjfIhS1soxZhHp8uLuw8jMxo3By/3Vi+I6SxrleC+mcoxZRLq82K/DyETS52HoAjcRkbS6Qh9GfsRxC41iXdchIlJCyv9eUsXuCNaDfkSkmyqvGkZFh3Dj6AjWNRIi0k2VV8Koro6/I1h3lxWRbqq8mqSGDIEsH6CUN6mesT1kSPFjEREpovKqYZSCm26CXr1eO3zv3vx1fqtTXURKkBJGpqZOhYEDXzu8qSk//RhtneqZXm8iIlJgShjZ2L07+fB89GOoU11ESpQSRjYKea8ndaqLSIlSwshGIe/1pBsPikiJUsLIRiHv9aQbD4pIiSqv02pLydSphbkGpG2ZM2cGzVBVVUGy0FXkIhIzJYxSVKhkJCKSAzVJiYhIJEoYIiISiRKGiIhEooQhIiKRKGGIiEgkShgiIhKJEoaIiEQSW8Iws/FmVmdmy81smZldHVcsIiLSuTgv3GsGrnH3J83sKGCxmT3s7stjjElERFKIrYbh7lvc/cnw/V5gBTA2rnhERCS9kujDMLMJwGnAwiTjpptZo5k17tixo+ixiYhIIPaEYWYDgN8Cn3f3lzuOd/dZ7l7j7jXDhw8vfoAiIgLEnDDMrBdBsqh199/FGYuIiKQX51lSBtwJrHD3/4krDhERiSbOGsY5wGXAO8xsSfi6MMZ4REQkjdhOq3X3xwCLa/0iIpKZ2Du9RUSkPChhiIhIJEoYIiISiRKGiIhEooQhIiKRKGGIiEgkShgiIhKJEoaIiESihCEiIpEoYYiISCRKGCIiEokShoiIRKKEISIikShhiIhIJEoYIiISiRKGiIhEooQhIiKRKGGIiEgkShgiIhKJEoaIiEQSa8Iws/PNbJWZrTGzL8cZi4iIpBdbwjCzHsCPgQuAE4FLzezEuOIREZH04qxhnAmscfd17t4E3AtcHGM8IiKSRs8Y1z0WeC7h8/PAWR0nMrPpwPTw4yEze6YIseVqGLAz7iAiUJz5Uw4xguLMt3KJc1I+FhJnwojE3WcBswDMrNHda2IOqVOKM7/KIc5yiBEUZ76VU5z5WE6cTVKbgfEJn8eFw0REpATFmTD+ARxnZseYWW/gEuBPMcYjIiJpxNYk5e7NZnYV8DegB3CXuy/rZLZZhY8sLxRnfpVDnOUQIyjOfOtWcZq752M5IiLSxelKbxERiUQJQ0REIimZhNHZbULMrI+Z3ReOX2hmExLGfSUcvsrMzosxxv8ws+VmttTMHjGz6oRxLWa2JHwVtHM/QpyXm9mOhHg+lTBumpmtDl/TYo7z1oQYnzWzlxLGFWV/mtldZrY91fU/FvhhuA1Lzez0hHHF3JedxTk1jO9pM5tvZm9MGLchHL4kX6df5hDnFDPbk/DdfiNhXNFuJRQhzi8lxPhM+HscEo4ryv40s/FmVheWOcvM7Ook0+T39+nusb8IOr3XAhOB3sBTwIkdpvkscFv4/hLgvvD9ieH0fYBjwuX0iCnGtwOV4fsZbTGGn/eV0L68HPjfJPMOAdaFfweH7wfHFWeH6f+d4MSIYu/PtwGnA8+kGH8h8FfAgDcDC4u9LyPGeXbb+glux7MwYdwGYFiJ7M8pwF9y/b0UOs4O074XeLTY+xMYDZwevj8KeDbJ/3pef5+lUsOIcpuQi4HZ4fv7gXeamYXD73X3Q+6+HlgTLq/oMbp7nbsfCD8uILi2pNhyueXKecDD7r7b3V8EHgbOL5E4LwXuKVAsKbn7XGB3mkkuBn7hgQXAIDMbTXH3Zadxuvv8MA6I77cZZX+mUtRbCWUYZ1y/zS3u/mT4fi+wguAOGony+vsslYSR7DYhHTf8yDTu3gzsAYZGnLdYMSa6giCzt+lrZo1mtsDM3l+A+NpEjfNDYRX1fjNru4CyWPsyo3WFTXvHAI8mDC7W/uxMqu0o5r7MVMffpgMPmdliC27FE7e3mNlTZvZXMzspHFaS+9PMKgkK2t8mDC76/rSgif40YGGHUXn9fZb8rUHKkZl9DKgBJicMrnb3zWY2EXjUzJ5297XxRMifgXvc/ZCZfYag5vaOmGKJ4hLgfndvSRhWSvuzbJjZ2wkSxlsTBr813JcjgIfNbGV4hB2HJwm+231mdiHwB+C4mGKJ4r3A4+6eWBsp6v40swEECevz7v5yodYDpVPDiHKbkCPTmFlP4GhgV8R5ixUjZvYuYCbwPnc/1Dbc3TeHf9cB9QRHA4XQaZzuvishtjuAM6LOW8w4E1xChyp/EfdnZ1JtR8nd+sbM3kDwfV/s7rvahifsy+3A7ylMk24k7v6yu+8L3z8I9DKzYZTg/gyl+20WfH+aWS+CZFHr7r9LMkl+f5+F7piJ2HnTk6DT5Rhe7dA6qcM0/4/2nd6/Dt+fRPtO73UUptM7SoynEXTMHddh+GCgT/h+GLCaAnXYRYxzdML7DwAL/NWOsPVhvIPD90PiijOc7gSCTkSLY3+G65hA6k7a99C+U3FRsfdlxDirCPr3zu4wvD9wVML7+cD5McY5qu27JihoN4X7NtLvpVhxhuOPJujn6B/H/gz3yy+A76eZJq+/z4Lt7Cw2/kKCXv61wMxw2DcJjtQB+gK/CX/0i4CJCfPODOdbBVwQY4x/B7YBS8LXn8LhZwNPhz/yp4ErYt6X3wKWhfHUASckzPvJcB+vAT4RZ5zh5+uBWzrMV7T9SXD0uAU4TNDOewVwJXBlON4IHgS2NoylJqZ92VmcdwAvJvw2G8PhE8P9+FT4m5gZc5xXJfw2F5CQ4JL9XuKKM5zmcoITbhLnK9r+JGhWdGBpwvd6YSF/n7o1iIiIRFIqfRgiIlLilDBERCQSJQwREYlECUNERCJRwhARkUiUMEREJBIlDBERiUQJQyQH4fMI3h2+/y8z+1HcMYkUim4+KJKb64BvhjeaOw14X8zxiBSMrvQWyZGZNQADgCkePJdApEtSk5RIDszsFIInnzUpWUhXp4QhkqXwyWW1BE8122dmBXuinkgpUMIQyUL4pLXfAde4+wrgRoL+DJEuS30YIiISiWoYIiISiRKGiIhEooQhIiKRKGGIiEgkShgiIhKJEoaIiESihCEiIpH8f4iocrHeUNFwAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
    " + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Using Autograd to calculate gradients for OLS\n", + "from random import random, seed\n", + "import numpy as np\n", + "import autograd.numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from autograd import grad\n", + "\n", + "def CostOLS(beta):\n", + " return (1.0/n)*np.sum((y-X @ beta)**2)\n", + "\n", + "n = 100\n", + "x = 2*np.random.rand(n,1)\n", + "y = 4+3*x+np.random.randn(n,1)\n", + "\n", + "X = np.c_[np.ones((n,1)), x]\n", + "XT_X = X.T @ X\n", + "theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)\n", + "print(\"Own inversion\")\n", + "print(theta_linreg)\n", + "# Hessian matrix\n", + "H = (2.0/n)* XT_X\n", + "EigValues, EigVectors = np.linalg.eig(H)\n", + "print(f\"Eigenvalues of Hessian Matrix:{EigValues}\")\n", + "\n", + "theta = np.random.randn(2,1)\n", + "eta = 1.0/np.max(EigValues)\n", + "Niterations = 1000\n", + "# define the gradient\n", + "training_gradient = grad(CostOLS)\n", + "\n", + "for iter in range(Niterations):\n", + " gradients = training_gradient(theta)\n", + " theta -= eta*gradients\n", + "print(\"theta from own gd\")\n", + "print(theta)\n", + "\n", + "xnew = np.array([[0],[2]])\n", + "Xnew = np.c_[np.ones((2,1)), xnew]\n", + "ypredict = Xnew.dot(theta)\n", + "ypredict2 = Xnew.dot(theta_linreg)\n", + "\n", + "plt.plot(xnew, ypredict, \"r-\")\n", + "plt.plot(xnew, ypredict2, \"b-\")\n", + "plt.plot(x, y ,'ro')\n", + "plt.axis([0,2.0,0, 15.0])\n", + "plt.xlabel(r'$x$')\n", + "plt.ylabel(r'$y$')\n", + "plt.title(r'Random numbers ')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "8e6f5c03", + "metadata": {}, + "source": [ + "## Including Stochastic Gradient Descent with Autograd\n", + "In this code we include the stochastic gradient descent approach discussed above. Note here that we specify which argument we are taking the derivative with respect to when using **autograd**." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "2a9f764d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Own inversion\n", + "[[3.85362611]\n", + " [3.15660109]]\n", + "Eigenvalues of Hessian Matrix:[0.37174525 3.67841217]\n", + "theta from own gd\n", + "[[3.86040722]\n", + " [3.14992207]]\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsNklEQVR4nO3deXxU9bnH8c+TsAZc2EUkCW4o7pJW63WhVXHXrvdqse5S7W691Vba2kXtYm/tbi+1uBSqtYu9tWoNAopCEYKAiCCgAqJhV1nCluS5f5xJHOLM5CQ5M2dm8n2/XvPKzNnmN2cmv+f81mPujoiISGtK4k6AiIgUBgUMEREJRQFDRERCUcAQEZFQFDBERCQUBQwREQlFAUOKlpl9x8wmxp2OKJnZKDNbHXc6pHNSwJCcMrMVZrbdzLaa2Rozu8/MesedLhFpnQKGxOECd+8NHAscB3wj3uR0TmbWJe40SGFRwJDYuPsa4EmCwAGAmX3dzF41sy1m9rKZfSxp3RVm9pyZ/cTM3jaz183snKT1w8zsmcS+k4H+ye9nZhea2SIze8fMnjazw5PWrTCzr5nZi2a2zcx+b2aDzOyJxPGeMrM+qT5HUzWRmd1oZuvMrNbMrkxa/7SZXdPycyS9djP7nJktS7zX983sIDObaWabzexhM+vW4j1vMbMNiXSPSVrePXF+VpnZWjP7rZn1bJHOm81sDXCvmfU3s38mzskmM3vWzJQvSEr6YUhszOwA4BxgedLiV4FTgH2A7wITzWxw0voTgFcIgsGPgd+bmSXW/RGYm1j3feDypPc6FHgQ+AowAHgceLRFRvwJ4EzgUOAC4AnglsT2JcCXMnyc/RJpHgJcDfw6XYBJ4yxgJHAicBMwHrgUGAocCVzS4r36J97rcmC8mQ1PrPthIv3HAgcntvl2i337AhXAWOBGYHXiMw5KfF7NFyQpKWBIHP5uZluAN4B1wK1NK9z9z+7+lrs3uvufgGXAB5P2Xenuv3P3BuB+YDAwyMzKgQ8A33L3ne4+HXg0ab//Ah5z98nuvhv4CdATOClpm1+6+1p3fxN4Fnje3ee5+w7gEYLqs3R2A99z993u/jiwFRieYfuWfuzum919EfASUO3ur7n7uwSBq+V7N33OZ4DHgP9MBM6xwA3uvsndtwB3ABcn7dcI3JrYd3si3YOBikTan3VNMCdpKGBIHD7q7nsBo4DDSKo6MrPLzGx+oorkHYKr6+SqpTVNT9y9LvG0N7A/8La7b0vadmXS8/2TX7t7I0HAGpK0zdqk59tTvM7UOL/R3euTXte1sn1LbXnvVJ9zf4JSQhkwN+n8/SuxvMn6RABscidBCa/azF4zs6+3Ic3SyShgSGwSV8f3EVztY2YVwO+ALwD93H1fgqttS3OIZLVAHzPrlbSsPOn5WwTVMCTeywiqe95s/ycIbRtBRt5kvw4eL9XnfAvYQBBcjnD3fROPfRIdDJrsUXpw9y3ufqO7HwhcCHzVzE7vYPqkSClgSNx+BpxpZscAvQgytPUAiYbjI8McxN1XAjXAd82sm5mdTNAO0eRh4DwzO93MuhLU3e8EZkb1QTKYD3zczMrM7GCCNo6OavqcpwDnA39OlJp+B9xlZgMBzGyImZ2V7iBmdr6ZHZwIoO8CDQTVViLvo4AhsXL39cADwLfd/WXgf4B/E1TJHAXMaMPhPk3QKL6JoF3kgaT3eYWgEfmXBFfiFxB0790VwcdozV3ALoLPdD8wqYPHWwO8TVCqmARc5+5LEutuJqhimmVmm4GnyNyWckhim60E5/037j6tg+mTImVq3xIRkTBUwhARkVAUMEREJBQFDBERCUUBQ0REQimoycf69+/vlZWVcSdDRKSgzJ07d4O7D2h9y8wKKmBUVlZSU1MTdzJERAqKma1sfavWqUpKRERCUcAQEZFQFDBERCQUBQwREQlFAUNEREJRwBARkVAUMEREJBQFDBERCUUBQ0REQlHAEBGRULIeMMxsgpmtM7OXUqy70czczPpnOx0iItIxuShh3Aec3XKhmQ0FRgOrcpAGERHpoKwHDHefTnCP5ZbuAm4CdI9YEZECEEsbhpldBLzp7gtCbDvWzGrMrGb9+vU5SJ2IiKSS84BhZmXALcC3w2zv7uPdvcrdqwYM6PB07iIi0k5xlDAOAoYBC8xsBXAA8IKZ7RdDWkREJKSc30DJ3RcCA5teJ4JGlbtvyHVaREQkvFx0q30Q+Dcw3MxWm9nV2X5PERGJXtZLGO5+SSvrK7OdBhER6TiN9BYRkVAUMEREJBQFDBERCUUBQ0REQlHAEBGRUBQwREQkFAUMEREJRQFDRERCUcAQEZFQFDBERCQUBQwREQlFAUNEREJRwBARkVAUMEREJBQFDBERCUUBQ0REQlHAEBGRUBQwREQkFAUMEZFcmTQJKiuhpCT4O2lS3Clqk6zf01tERAiCw9ixUFcXvF65MngNMGZMfOlqg6yXMMxsgpmtM7OXkpbdaWZLzOxFM3vEzPbNdjpERGI1btx7waJJXV2wvEDkokrqPuDsFssmA0e6+9HAUuAbOUiHiEh8Vq1q2/I8lPWA4e7TgU0tllW7e33i5SzggGynQ0QkVuXlbVueh/Kh0fsq4Il0K81srJnVmFnN+vXrc5gsEZEI3X47lJXtuaysLFheIGINGGY2DqgH0nYVcPfx7l7l7lUDBgzIXeJERKI0ZgyMHw8VFWAW/B0/vmAavCHGXlJmdgVwPnC6u3tc6RARyZkxYwoqQLQUSwnDzM4GbgIudPe61rYXEZGEGMdyZL2EYWYPAqOA/ma2GriVoFdUd2CymQHMcvfrsp0WEZGCFvNYDiuk2qCqqiqvqamJOxkiIvGorAyCREsVFbBiRdrdzGyuu1d19O3zoZeUiIiEEfNYDgUMEZFCEfNYDgUMEZFCEfNYDgUMEZFCEfNYDs1WKyJSSGIcy6EShohIXArs/hgqYYiIxKEA74+hEoaISBwK8P4YChgiInHI5piKFlVd/aFvxw+qgCEiEo9sjaloqupauRLcYeVKyqGiYwcNKGCIiMQhW2MqUlR1WUR5vQKGiEgcsjWmIovThKiXlIhIXLIxpqK8PPUEhRFQCUNEpJice27WDq2AISJSTB5/PGuHVsAQESkmWWzDUMAQESkmWZzqXAFDRKSYpOiu69AYxaEVMESkcyuwCQBblaK77iqIpNuUAoaIdF4pRkUzdmxxBI0VK9j6Tj2P3fAUG2BTFIfNesAwswlmts7MXkpa1tfMJpvZssTfPtlOh4jI+xTgBICZNDbC3CfW8YNPzuXDgxbRd596zv/KwZEdPxcljPuAs1ss+zowxd0PAaYkXouI5FY2JwDMkdVL67j3hhe5ZPgLDOq2iapzB3LLX0fyzqZGbjiimsk3TY7svbI+0tvdp5tZZYvFFwGjEs/vB54Gbs52WkRE9pBuVHQWexp11LYtjUy/7zWq//Q21fMH8PK2SuBo9qOWcwfPZfSoXZxx7TAGjToS7Khgpx9H895xTQ0yyN1rE8/XAIPSbWhmY4GxAOV5/CWKSAG6/fY9b2IE0UwAGKHGRlgweR3V96yi+tkePLf2EHZxMD3Yzqm953HVaQsZ/en+HHnpsVjZmVlNS+xzSbm7m5lnWD8eGA9QVVWVdjsRkTZrmsdp3LigGqq8PAgWMd/x7q3ldUy+eznVj+1m8vJK1jcMBAZydJdFfGnEFEaf342TrzuSnsNOymm64goYa81ssLvXmtlgYF1M6RCRzi4bEwC2Ud3WRp6979Wgmmlef17adiBwNANZy1mD5zH6tF2ccU0lgz8yAuyI2NIZV8D4B3A58MPE3/+LKR0i+W/SpLy7ApaOaWyEhVPWBtVMz/Tg2bWHsJND6M4OTuk9j8tOfZnRl/TjqM8cS0mvM+JObrOsBwwze5Cggbu/ma0GbiUIFA+b2dUEA0r+M9vpEClITeMEmurYm8YJgIJGPgkR1Ne8Vsfk3yxrrmZaWz8IGMSRXRbz+cOnMfq8rpxy3RGUHfSheD5DCOZeOM0CVVVVXlNTE3cyRHKnsjJ1L56KClixItepkSbJAaJvX9iyBXbtem99WRm7fvm/PLPzRKof2kj1C/15cetBAAxgHWfu9xKjT9vJmddUsP/phwcjsrPIzOa6e1VHjxN7o7eIZFAE4wSKTstS38aN79+mro7aq29hNKvoxlBO7j2fH56yhNGX9OOYy46hpNdHspu+lqWdiGhqEJF8lq4reTa6mBfbnErZkmp0eArlrObxG6ewafnbTNlyAjdPP4/jrj+Rkl49s5e2NFOd9Ie+URxeAUMkn6WYeTQr4wSKdU6liO2oa8RXhivdWUU55/zkdHodtF+WU5UkzVQn+8OQKA6vgCGSz1LMPMr48dE3eBfZnEpRcYeXpqzlp5fM5uzBC+jbawcrCVG6i2vwX5qqyq7QLYrDq9FbRIJqqFR5gVnQB7QTWff6Np66exnVj+2iemkltfUDATi8dCmjD3mNyw+czrGP/4C0zdQVFfF1fU7TSeIY2LXAvXtHD69GbxEpyDmVorJzeyMzHniV6gc3Uj23H/O2HgIcS182cuaghYw+9QXOvKaCoWceBnYocDb0H5+6sTvu3mtppjp5q67uzSgOryopEcldW0kecIeXp63lZ5+ezbn7z6dv2XZOv+4Q/ueZkezt7/DIoTexo89gNtgAHupxBVddtJGho1t0ff35z/PzfKWpwozqfhi4e8E8Ro4c6SKSJRMnuldUuJsFfydOjDtFkVm/Yqs/ePM8v/KIWT6kS60HYcN9eOlS/+LwJ/3RG6b45qW1wWcuK/PmDSB4nepcFND5Amo8gjxYbRgiUthSjDvY9YlLmPmHV6n+4waq5/blhS2H4JTQh02cMWgho0/ZwZlXl1Nx1mF7lhyKdKBkVAP3VCUlItmVzfEdKboD7/jMNXy25318eOwh3Pl0FWWNW/nef1Tz/M9nsX5LTx5ecxrX/PksKs5OMcI6FwMlC3i8ixq9RSR7sjgX1sY36ujxua/Rq0V34B6+g590+QYf+8IwRl1/OHsfOjL8QbPd+F/gc4OphCEi2fPlL0c2vqP+vons6DcEtxJqS4bwpfK/0XPzmpTb9mtYz4V3fZi9D23joLlsN/4X+HgXBQwRyY5Jk1J3PYVQVTzusPTZtfzq0ln8pM9t7LpyLD02vYXhDPa3uLfkGhp67Z165/aWCMIMlOxIlVKBzw2mKikRyY5MV81pMvRNq+uYevcrVP9jB9WvlLNy9xBgEKv5BGVs32Pbbo07oUdv8LJob7Ga6YZKHa1SKvDxLiphiEjbhbnKznTVnMjQd+9s5LkJS/n26TM4cZ+XGTC0O5+64zj+9NIIju/zOnd/fDLL/7mEIVab+jibNuVm6pQmHa1SKvTxLlH0zc3VQ+MwRPJA2LEKFRV7bpN4NGLeiPna0v38Su5xcC+h3j9UNs9v/dCTPuOns3z35rpQx/KKihx96ASz1OkwC3+MGMZvENE4jNiDQFseChgieSBs5p0isDS22GeHdffnR4/zrT/+deZMtC0D6rIpXwJXG0UVMFQlJSJt00rDbf2uRmbcu5Rb763k1sZbWUk5jRj1lL5vwr7uvpMPzv0tvb7ztcxTq+dq1t7WFHqVUgcpYIhI26RpoN1SNoiPDZ1Dvx5bOfmqQ7ltyolU21lMOPF3/Punsyi1NLPebtwYrl1gzJhgtHVjY/C3o8GiPb2d8iVwxURTg4hI20yahF97Lbb9vV5L2yjjWsYzs+QUzjpwOaPPKeUj1w+nz+FJ4yDSTbuRTjanVm/Z2wmCkkKRZv5RTQ0Sa8AwsxuAawAHFgJXuvuOdNsrYIjEo35XI3P+uIzqSeupnrMvw96dx218i3JWsaFkIPOPu5Jh376Mg88bjpWmqbhIl0n37Jn7qcKLdM6odKIKGLGNwzCzIcCXgBHuvt3MHgYuBu6LK00i8p7XZ62l+revUT21lCmrh/OuD8c4hA/0XMSwE/fjjU8+zJBrj2Lg3j0ZHeaATVfuLSYKBFIHkmy2CxT4ALq4xD1wrwvQ08x2A2XAWzGnRzqLFDOcFmNVRFtsXlPHtLuXUP1/dVQvHsryXRXAIIaWrOZTB73A6LNL+Mh1w+l3xFHAUe17k0yD4nL5fRT4ALrYRNHVqr0P4MvAVmA9MCnNNmOBGqCmvLw8ok5m0qllq4tmy/71118fvr99e/rmd7A/f/2uBp91/xL/3pnT/eR9XvRSdju492KLn99/pv/iosm+5O+LvbG+oU3HLQj50k03Ryj0cRhAH2AqMADoCvwduDTTPhqHIZHIRl/6VBlQy0emG/G0NfNqZ4b3+qw1Pv6KGf7JobN8X3s7GHNGg1f1XOi3nDDZn75ztu98py7jMYpGAd0AqaOKIWB8Cvh90uvLgN9k2kcBQyIRxWjdltIFoTBBqT0BLN0+JSV7ZHyb12zzf9xa4184drof2u215s0OKFntVx30tD/0uWd8/cLa9n/ufNKJAkBbRRUw4hyHsQo40czKzMyA04HFMaZHikmmPvbp6qk7Un8dtrE0Vb15expg061rbKTxiit55OhvcVqfBfTdrysXfnckE+Yfz8F7rePx429h54AhrPKh/L7+cv7rpDfof2QbpwDPRylupPS+wX/ScVFEnfY+gO8CS4CXgD8A3TNtrxKGhNJadU026q/DljBKS8Pv254SRuLxOhV+fI9F/vUPPuVTfzTbd7xdV9z19gU6ZUeuUOhVUu15KGBIKGEyj6irL8K0YTQ9wuybISPfsnabz/349983L1PyozFV9VoxZ6rZqGYsIjkLGMBk4Jgo3qyjDwUMCSWuzCM5CJWWti1zzhDAGnY3+JyJi/32s57x0/ad713Z6eC+nn7pg1Kq9ynmTLWYg2EEogoYYdowbgZ+Zmb3mtngbFSLiUQqG20UYSTPdXT//a1PUpfczjJuXLAuMU/SG4eezoSrZ3Bxxb8Z2O1tPnDpYYx78lQ27+jGVz/wHFN+VMPe/3sndO36/nR065Z60Ftc5yUXOvmkgDkTNrIAnwAWALcCPaOIVm19qIQhoeRLXX2maq8UadzdtYffX/FNP7z78ubFg0tq/fIDp/uk66b72gUpejNNnOjeL6mk0a9f+s+ZL+clW9RLKi1y2YYBGHAkcB2wAVgNfCaKBLTloYAhoeV55tFYXv7+6hPwFQz1s/rN9v85f6ov/MuS6AfN5fl5keyIKmC0Ovmgmc0AhgGLgFnA8wQ9m75M0KtpbPTlntQ0+WAWaIqMnHlz7hom372c6qdKmLjyZEp4//+em2HZmqE1LP0mik5Ukw+GKV0cQWJW2xTrFkcRtcI+VMKIWLFXUbjHekW9bf02f+L7c/yG46f5Ed2XNp/iQSVrfX3poJQljNgbaTvDb6ITIh+61QIHRpGIsA8FjIjF0f00l3Kc+TXsbvB5Dy3xH50zzU/vM9e7scPBvTvb/cy+NX7nedN8wcOJaqZ8zZjV26goRRUwdAOlzqykJMgOWmq6cU2h32QmB/c8qJ23hsl3L6N6sjF55XDW+QAAjuqxlNFHvMXoT+zFKWMPp2e/svfvnI9VP639JqQg5axKKp8eRV3CiONKvrWrybivNjt6TrIw7qBuwzZ/8rbZfuPIaX5U91eaDznQ1vmYyuf8/muf9TfnFvDcTHF/55IV5EOVVK4fRRsw4qqeaO194xzoFcU5iSDza6xv8AV/Wux3njvVz+w7x7uz3cG9Gzv89L5z/UfnPu3zHlriDbvzYArwKC468rWqTDpEAaOQtPaPHOdVXaa0xZmutrx3us/Qzsyvdl6t/2HsdP9M5XQfVLK2edcjui/zG0Y+7U/cVuPb1m+L+AN3UJQZfSG3W0lKChiFIsw/cr5O2RDn1WbYcxJmosFWMr+6Ddu8+o45/t8jp/oxPRY3H6a/rfdLKmb4vdc856tr8ryaSVVJkoECRqEI84+cz//sUV5ttuVYYc9JO85dY32Dv/jwYv/JeVN9dN853oM6B/eu7PQP93nBf3DO0z73wVfyo5oprHy96JC8EFXAiPue3sUvzL0Obr89dW+kYpoHp2WPq6b7FUDqnkFhz0nIe0msXbCGp36zlOrJMHnlodQ2HgYcxojuy/nscXMY/bFenHbd4fQacFz7Pl/cdI9qyYUook6uHkVbwnDPz3rjKKuk2lOKCnNO0hy34YCh/tQPZvtNVVP82B4vN6/qZxv84oqZPuHq5/yNOXlezdQWaqyWDFCVVIEohH/kdBlzlFVl2aoySXF+t9PDL2dCczXTqD7z/I6zn/GaSXnSmylb8vGiQ/KCAkYhyed/5EwBrbVMPhttEm2w7sVa/+N1z/hvB3zT32CIN2D+OhV+Y+lP/UvHPeP//G6Nb1mbZ72ZRGKggCFtlyqDz5SRZ1rX1pJTBCWtHZu2+dQfzfavf+ApP77HoubD9LWN/p/l//Z7rprhK2ev6fBpEik2nTtg5PMVe75Kl2GnCghNpYhMmXy6YFJamv57aeP31ljf4C//9WX/2QVP+bn9ZnkZWx3cu7DLT913gd921nSfPXGJ1+/Kg2qmttyXoml7/YYlRzpvwCiENoF8lCmDz1RVlC5jS1dd1cHvZf3CWn/oc8/4VQdN8wNKVjcf6quld/nGroO8EfOGA4a2/fvOZgY9caJ7t27v//xdu6Z+H/2GJceKImAA+wJ/Ibi/xmLgQ5m2HzlyZH6PWchnmTL4VJnX9de3b3R6G7+XnW9v86d//Lzf8sHJXtVzoRsNDu772tv+yfJZPv7Kmb7uO7/uWAab7Qw607lI9fn1G5YcK5aAcT9wTeJ5N2DfTNuPHDlSA5Taq7X2iOTgcP31rWewqTLhEN9LY32DL3nkZf/FhZP9/P4zvRdbgoIOu/3kfV/0741+1mc98Mqe1UwdzWCznUFnCsapfpf6DUuOFXzAAPYBXifNzZlSPVTC6IC2XGW3Z+xIhqqtjYtq/eHPP+3XHDTVy0tWNa86uNsK/9wxz/nfb53r79Zm6M3U0Qw22xm0ShiS54ohYBwLzAbuA+YB9wC9Umw3FqgBasrLy1X/2xFh6/Hbk8Gm+F52lnT3b3a9o7maaR97xz8+9Hn/7RUz/dWZbejNlM8ljJaN3ckPtWFIniiGgFEF1AMnJF7/HPh+pn3USypH2pjBNtY3+NK/L/Inj7/Z15YMah4PMYYH/KR9XvLvnPGsz7zvFd+9s529mTqawWYrg85ULadeUpJHiiFg7AesSHp9CvBYpn00DiNHQmSwm16u9T9/YZpfe/AUryxZ2bzZgV1X+XVHz/C/fesFf/vNdg6aS5WZdjSDzUYGraolKRAFHzCCz8CzwPDE8+8Ad2baXgEjh1pksLt/N8Gf/cks/9YJT/oJPed7CfUO7nvbu/7RA+b4by77ty+fEcGguThvJtXWgKLGaykQxRIwjk20T7wI/B3ok2n7ThswYqi+aKxv8GX/t8h//dEn/aIBz/levOvgXkK9f2ifl/zWM57z5+5d6rt2RDxoLtvtDRHeaEklDCkURREw2vqIJGAUWt1xDq+4315c63/94jT/7MFP+bCSFc1vV9n1Df/s0TP8r998wTetzvLcTDmcpLDVUeutZfxqvJYCoYDRHrn+B48iOIXJzNr5Prvf3eYzfjrLbz3xX/6hnvOaq5n2ss1+0QFz/NeXzfJlz63xxsa2J7vdsnXVnum4HQlS+XwBks9pk5xSwGiP9sx/1F5RBacwM8aGfZ+GBn/10UV+98ee9I8NeNb35p3maqYT9n7Zv3X6DH92QhaqmdoiW0E903ksxqollX4kiQJGe2Rp/qOUosqEWjtOK+vfWVLrj3xpql9/yGQ/qPS15tXlXd/0a4+a6X8eN883vpFnU4DnukdTMWauxRgEpd0UMNojovmPQomqLj5dX/+mfv5p3qcR8/8om+ul7HZw721b/IIhc/2Xl87yV6bnuJopH7QWFIqt+kY9uCSJAkZ7tHP+o7THas/kfO0JRulGE5eVpR1l/Drl/oG9F/u4j8zwZ+5Z6ju358EU4NkUJsMvtqCQiUoYkkQBo71Czn/U6jEyXfWn28YsmNivPdJkABvZ17ey5/vs7tLDN/98Quufv1gyzWKsUuoonRNJooARhaj737fc//rr31810I5/2vrN27yR1FUMDZg/cdJ3fNegA7yxtSBQrJmIrqZTK8aLA2kXBYyoRDnCt2VG1d6MrKHBVzy+yH/3iSf8UwOf9j5s9Ndp57GSxZ2xZisDU329SEYKGHFqrfG8KaNqQ0a2eWmtP3rDU/7F4f/y4aVLmzcd0mWNX3nE8z7jgju8sWfPPY+TqnSQKVOOM2PNZukm7kAokucUMHJR3E4eBdzU3pHuBkNtLGHUb97mc34x028/+TE/tdcc78IuB/eebPNzBs/zuy553hdNW7tnb6bWPnNrmXKcGWu2p/woxqo2kYh07oCRiwwiU4+qpluYpuu5lKHhe1dpd//53uO8H+ubFx+311K/edS/fcr/LvMddR3ozdRaphxnxprt0o3q60XS6twBIxdXyq1VOyVnwmkyqi3Lav2Fs7/hG7sMbL5HxCVM9MFd1vrlI2b7H2+e7+/cNSG6jC5MphxXxqpqI5HYdO6AkanROSqtNWynuDJu2FrnNb+a6Xec/JiP6vW8d2VnczXT2YPn+08vme0LpyRVM4XpntuaKLoJd0TYAKRqI5HYdO6Ake5q1Sy6DChMCaOhwd94cpFP+NTjfvGgqd6fdc2rj+29zG86bZY/9dtlvn1bmmqmsN1z0wkzEDHbkyu2JQio2kgkFp07YGSYEiOyq+kMmfHuLt39gcFf8xGlS5oX79dlnV82YrZP/Np8X/PatveOkSmDDNs9N51cTqbYlvdXNZNIXokqYHShEI0ZA5demnrdqlUdO/akSTBuXHCcvn3xHj1g0yacEqCRNyjnG/V38Ejtxzl1v2VcddocRl9TzpGnD8JswJ7HGTsW6uqC1ytXBq+b0g9QXh4sT6e1z5JufUMDmIX6uB2S7v07+h2ISF4qiTsB7VZRkXp5eXn7j9mUya9cGVwrb9zI9k3bGcNESmnguN7L+dWpf+bKu09g09buPFl7NDdesJSjrjkBKy2BysrgGBAEnaZg0aSuLlje5PbboawsfXpa+yyZ1ru/F6Sa0hS1dO/fke9ARPJXFMWUXD3eNw4jokbUbcvf8n/992Tf2GVAyiqWLXsP9reWp5gCPFMawnYjzTSxYBRtGNmsIlJDtkhBoFO3YTRp2UZw/fWhGlUbttb5/N/M8B+f+qif0Xumd2e7g3tDG3pEuXvmOvy21u+3t0E4eb+2pj8KasgWyXtRBQwLjlUYqqqqvKamJvXKlm0GEFT3jB8Pn/40a6a+zOTfraD6me5MXnMka9kPgCN7vc7o4zdw2YHPcfQDNwZT/LVUUQErVrx/eUlJkCW3ZAZ/+EP69DS1YUStsjJ1m0i69ItIp2Bmc929qsMHiiLqdOQBlALzgH+2tm3GuaTSXNFv6jLAj+7yUvOiAaUb/NOH1fh9X13gby6va3X/jF11w4yszuXVt6qIRCQFiqiX1JeBxcDeHTmIr1pFqn5B+9RvoP8g54enzmX01UM55syBlJT023OjSZPS91ZyT18iuP321KWI228Pno8Zk73SRCpN79XUy6u8PEhLLtMgIkUr1l5SZnYAcB5wT5t3dmft1EVMvOQxLt+/mtW+f+rthg5lypojufnhkRx31kBKWn7ipqqsdPr1S79uzJigiqmiIqiGqqiIvspp0qSgqqmkRS+sTGlasQIaG4O/ChYiEpUoiintfQB/AUYCo0hTJQWMBWqAmor9D/DJN1X710b804/psrC51qV/6Ub/5X63+e4uPdpeHdPaiO5u3eKr0lEVU/ap0V46AQq90dvMzgfOdffPmdko4L/d/fxM+5TY8e68QFd2cfKgZYw+ZQf/NWw2lX/6EfZGMNAOgE2bWq+OaRqgl2ngXJO4Go3ViJ1dmTpKqGQmRSSqRu84A8YPgM8A9UAPgjaMv7l7miHcMKj3CL/3B49y6uXD6L13Sfv/4VPtlzmxQRVPrmXqhRVHeoqNArJ0EgUfMPZIRMgSRpWZ11RUvFdyaO8/fLr90lEJozgpIEsnEVXAKLypQZKnu0g3Z9HKlZkbhzPNddRyDqbkXk+5lmrqkDjTU2w0tYlIm+RFwHD3p1srXeyhaU6mTP/YmeZQSrdfRUUw4C6bvZ7aIhe9sDozBWSRNsmLKqmwqsy8eZx3utHUydJV3aixU5okz06scStSpIqqDSOsPQJGUzCYNCn9VOeZ6qKVUYhIJ9G5A0bL0oAah0VE0uq8jd6p6vFVFy0iknX5MJdUeCNHQqrZajWHkohI1hVeCSOdsHMotXVuJhERAQqthNFRYe6zLSIiKRVPCSOMMPfZFhGRlDpXwEg3wjvTyG8REQE6W8AIMxWE2jhERFLqXAGjte63TW0cK1cGk9Ilz1slItLJda6A0drcTGrjEBFJq7BGeldVeU2qcRhR0XTXIlKEOu9I72zSdNciImkpYCTTFCMiImkpYCTT/SdERNJSwGgp7BQjrVH3XBEpMp1rapBc0RQkIlKEVMLIBnXPFZEipICRDZqCRESKUGwBw8yGmtk0M3vZzBaZ2ZfjSkvk1D1XRIpQnCWMeuBGdx8BnAh83sxGxJie6Kh7rogUodgChrvXuvsLiedbgMXAkLjSEyl1zxWRIpQXU4OYWSUwHTjS3Te3WDcWGAtQXl4+cuXKlblPoIhIASuaqUHMrDfwV+ArLYMFgLuPd/cqd68aMGBA7hMoIiJAzAHDzLoSBItJ7v63ONMiIiKZxdlLyoDfA4vd/adxpUNERMKJs4TxH8BngI+Y2fzE49wY0yMiIhnENjWIuz8HWFzvLyIibRN7o7eIiBQGBQwREQlFAUNEREJRwBARkVAUMEREJBQFDBERCUUBQ0REQlHAEBGRUBQwREQkFAUMEREJRQFDRERCUcAQEZFQFDBERCQUBQwREQlFAUNEREJRwBARkVAUMEREJBQFDBERCUUBQ0REQlHAEBGRUGINGGZ2tpm9YmbLzezrcaZFREQyiy1gmFkp8GvgHGAEcImZjYgrPSIiklmcJYwPAsvd/TV33wU8BFwUY3pERCSDLjG+9xDgjaTXq4ETWm5kZmOBsYmXO83spRykraP6AxviTkQISmd0CiGNoHRGrVDSOTyKg8QZMEJx9/HAeAAzq3H3qpiT1CqlM1qFkM5CSCMonVErpHRGcZw4q6TeBIYmvT4gsUxERPJQnAFjDnCImQ0zs27AxcA/YkyPiIhkEFuVlLvXm9kXgCeBUmCCuy9qZbfx2U9ZJJTOaBVCOgshjaB0Rq1TpdPcPYrjiIhIkdNIbxERCUUBQ0REQsmbgNHaNCFm1t3M/pRY/7yZVSat+0Zi+StmdlaMafyqmb1sZi+a2RQzq0ha12Bm8xOPrDbuh0jnFWa2Pik91yStu9zMliUel8eczruS0rjUzN5JWpeT82lmE8xsXbrxPxb4ReIzvGhmxyety+W5bC2dYxLpW2hmM83smKR1KxLL50fV/bID6RxlZu8mfbffTlqXs6mEQqTza0lpfCnxe+ybWJeT82lmQ81sWiLPWWRmX06xTbS/T3eP/UHQ6P0qcCDQDVgAjGixzeeA3yaeXwz8KfF8RGL77sCwxHFKY0rjh4GyxPPrm9KYeL01j87lFcCvUuzbF3gt8bdP4nmfuNLZYvsvEnSMyPX5PBU4HngpzfpzgScAA04Ens/1uQyZzpOa3p9gOp7nk9atAPrnyfkcBfyzo7+XbKezxbYXAFNzfT6BwcDxied7AUtT/K9H+vvMlxJGmGlCLgLuTzz/C3C6mVli+UPuvtPdXweWJ46X8zS6+zR3r0u8nEUwtiTXOjLlylnAZHff5O5vA5OBs/MknZcAD2YpLWm5+3RgU4ZNLgIe8MAsYF8zG0xuz2Wr6XT3mYl0QHy/zTDnM52cTiXUxnTG9dusdfcXEs+3AIsJZtBIFunvM18CRqppQlp+8OZt3L0eeBfoF3LfXKUx2dUEkb1JDzOrMbNZZvbRLKSvSdh0fiJRRP2LmTUNoMzVuWzTeyWq9oYBU5MW5+p8tibd58jluWyrlr9NB6rNbK4FU/HE7UNmtsDMnjCzIxLL8vJ8mlkZQUb716TFOT+fFlTRHwc832JVpL/PvJ8apBCZ2aVAFXBa0uIKd3/TzA4EpprZQnd/NZ4U8ijwoLvvNLPPEpTcPhJTWsK4GPiLuzckLcun81kwzOzDBAHj5KTFJyfO5UBgspktSVxhx+EFgu92q5mdC/wdOCSmtIRxATDD3ZNLIzk9n2bWmyBgfcXdN2frfSB/Shhhpglp3sbMugD7ABtD7purNGJmZwDjgAvdfWfTcnd/M/H3NeBpgquBbGg1ne6+MSlt9wAjw+6by3QmuZgWRf4cns/WpPsceTf1jZkdTfB9X+TuG5uWJ53LdcAjZKdKNxR33+zuWxPPHwe6mll/8vB8JmT6bWb9fJpZV4JgMcnd/5Zik2h/n9lumAnZeNOFoNFlGO81aB3RYpvPs2ej98OJ50ewZ6P3a2Sn0TtMGo8jaJg7pMXyPkD3xPP+wDKy1GAXMp2Dk55/DJjl7zWEvZ5Ib5/E875xpTOx3WEEjYgWx/lMvEcl6Rtpz2PPRsXZuT6XIdNZTtC+d1KL5b2AvZKezwTOjjGd+zV91wQZ7arEuQ31e8lVOhPr9yFo5+gVx/lMnJcHgJ9l2CbS32fWTnY7Pvy5BK38rwLjEsu+R3ClDtAD+HPiRz8bODBp33GJ/V4BzokxjU8Ba4H5icc/EstPAhYmfuQLgatjPpc/ABYl0jMNOCxp36sS53g5cGWc6Uy8/g7wwxb75ex8Elw91gK7Cep5rwauA65LrDeCG4G9mkhLVUznsrV03gO8nfTbrEksPzBxHhckfhPjYk7nF5J+m7NICnCpfi9xpTOxzRUEHW6S98vZ+SSoVnTgxaTv9dxs/j41NYiIiISSL20YIiKS5xQwREQkFAUMEREJRQFDRERCUcAQEZFQFDBERCQUBQwREQlFAUOkAxL3Izgz8fw2M/tl3GkSyRZNPijSMbcC30tMNHcccGHM6RHJGo30FukgM3sG6A2M8uC+BCJFSVVSIh1gZkcR3Plsl4KFFDsFDJF2Sty5bBLBXc22mlnW7qgnkg8UMETaIXGntb8BN7r7YuD7BO0ZIkVLbRgiIhKKShgiIhKKAoaIiISigCEiIqEoYIiISCgKGCIiEooChoiIhKKAISIiofw/bJ0I69c7pOcAAAAASUVORK5CYII=\n", + "text/plain": [ + "
    " + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "theta from own sdg\n", + "[[3.87985928]\n", + " [3.17479642]]\n" + ] + } + ], + "source": [ + "# Using Autograd to calculate gradients using SGD\n", + "# OLS example\n", + "from random import random, seed\n", + "import numpy as np\n", + "import autograd.numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from autograd import grad\n", + "\n", + "# Note change from previous example\n", + "def CostOLS(y,X,theta):\n", + " return np.sum((y-X @ theta)**2)\n", + "\n", + "n = 100\n", + "x = 2*np.random.rand(n,1)\n", + "y = 4+3*x+np.random.randn(n,1)\n", + "\n", + "X = np.c_[np.ones((n,1)), x]\n", + "XT_X = X.T @ X\n", + "theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y)\n", + "print(\"Own inversion\")\n", + "print(theta_linreg)\n", + "# Hessian matrix\n", + "H = (2.0/n)* XT_X\n", + "EigValues, EigVectors = np.linalg.eig(H)\n", + "print(f\"Eigenvalues of Hessian Matrix:{EigValues}\")\n", + "\n", + "theta = np.random.randn(2,1)\n", + "eta = 1.0/np.max(EigValues)\n", + "Niterations = 30\n", + "\n", + "# Note that we request the derivative wrt third argument (theta, 2 here)\n", + "training_gradient = grad(CostOLS,2)\n", + "\n", + "for iter in range(Niterations):\n", + " gradients = (1.0/n)*training_gradient(y, X, theta)\n", + " theta -= eta*gradients\n", + "print(\"theta from own gd\")\n", + "print(theta)\n", + "\n", + "xnew = np.array([[0],[2]])\n", + "Xnew = np.c_[np.ones((2,1)), xnew]\n", + "ypredict = Xnew.dot(theta)\n", + "ypredict2 = Xnew.dot(theta_linreg)\n", + "\n", + "plt.plot(xnew, ypredict, \"r-\")\n", + "plt.plot(xnew, ypredict2, \"b-\")\n", + "plt.plot(x, y ,'ro')\n", + "plt.axis([0,2.0,0, 15.0])\n", + "plt.xlabel(r'$x$')\n", + "plt.ylabel(r'$y$')\n", + "plt.title(r'Random numbers ')\n", + "plt.show()\n", + "\n", + "n_epochs = 50\n", + "M = 5 #size of each minibatch\n", + "m = int(n/M) #number of minibatches\n", + "t0, t1 = 5, 50\n", + "def learning_schedule(t):\n", + " return t0/(t+t1)\n", + "\n", + "theta = np.random.randn(2,1)\n", + "\n", + "for epoch in range(n_epochs):\n", + "# Can you figure out a better way of setting up the contributions to each batch?\n", + " for i in range(m):\n", + " random_index = M*np.random.randint(m)\n", + " xi = X[random_index:random_index+M]\n", + " yi = y[random_index:random_index+M]\n", + " gradients = (1.0/M)*training_gradient(yi, xi, theta)\n", + " eta = learning_schedule(epoch*m+i)\n", + " theta = theta - eta*gradients\n", + "print(\"theta from own sdg\")\n", + "print(theta)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "e2389c3f", + "metadata": {}, + "outputs": [], + "source": [ + "import jax.numpy as jnp\n", + "from jax import grad, jit, vmap\n", + "\n", + "def sum_logistic(x):\n", + " return jnp.sum(1.0 / (1.0 + jnp.exp(-x)))\n", + "\n", + "x_small = jnp.arange(3.)\n", + "derivative_fn = grad(sum_logistic)\n", + "print(derivative_fn(x_small))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/doc/src/week39/sgd.py b/doc/src/week39/sgd.py new file mode 100644 index 000000000..34e56c1bf --- /dev/null +++ b/doc/src/week39/sgd.py @@ -0,0 +1,74 @@ +# Using Autograd to calculate gradients using SGD +# OLS example +from random import random, seed +import numpy as np +import autograd.numpy as np +import matplotlib.pyplot as plt +from autograd import grad + +# Note change from previous example +def CostOLS(y,X,theta): + return np.sum((y-X @ theta)**2) + +n = 100 +x = 2*np.random.rand(n,1) +y = 4+3*x+np.random.randn(n,1) + +X = np.c_[np.ones((n,1)), x] +XT_X = X.T @ X +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y) +print("Own inversion") +print(theta_linreg) +# Hessian matrix +H = (2.0/n)* XT_X +EigValues, EigVectors = np.linalg.eig(H) +print(f"Eigenvalues of Hessian Matrix:{EigValues}") + +theta = np.random.randn(2,1) +eta = 1.0/np.max(EigValues) +Niterations = 100 + +# Note that we request the derivative wrt third argument (theta, 2 here) +training_gradient = grad(CostOLS,2) + +for iter in range(Niterations): + gradients = (1.0/n)*training_gradient(y, X, theta) + theta -= eta*gradients +print("theta from own gd") +print(theta) + + +n_epochs = 50 +M = 5 #size of each minibatch +m = int(n/M) #number of minibatches +t0, t1 = 5, 50 +def learning_schedule(t): + return t0/(t+t1) + +theta = np.random.randn(2,1) + +change = 0.0 +delta_momentum = 0.3 + +for epoch in range(n_epochs): +# Can you figure out a better way of setting up the contributions to each batch? + for i in range(m): + random_index = M*np.random.randint(m) + xi = X[random_index:random_index+M] + yi = y[random_index:random_index+M] + gradients = (1.0/M)*training_gradient(yi, xi, theta) + eta = learning_schedule(epoch*m+i) + # calculate update + new_change = eta*gradients+delta_momentum*change + # take a step + theta -= new_change + # save the change + change = new_change +print("theta from own sdg") +print(theta) + + + + + + diff --git a/doc/src/week39/test.do.txt b/doc/src/week39/test.do.txt new file mode 100644 index 000000000..52a49347a --- /dev/null +++ b/doc/src/week39/test.do.txt @@ -0,0 +1,293 @@ +TITLE: Codes +AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo & Department of Physics and Astronomy and Facility for Rare Isotope Beams, Michigan State University +DATE: today +!bc pycod +from numpy import asarray +from numpy import arange +from numpy.random import rand +from numpy.random import seed +from matplotlib import pyplot + +# objective function +def objective(x): + return x**2.0 + +# derivative of objective function +def derivative(x): + return x * 2.0 + +# gradient descent algorithm +def gradient_descent(objective, derivative, bounds, n_iter, step_size): + # track all solutions + solutions, scores = list(), list() + # generate an initial point + solution = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0]) + # run the gradient descent + for i in range(n_iter): + # calculate gradient + gradient = derivative(solution) + # take a step + solution = solution - step_size * gradient + # evaluate candidate point + solution_eval = objective(solution) + # store solution + solutions.append(solution) + scores.append(solution_eval) + # report progress + print('>%d f(%s) = %.5f' % (i, solution, solution_eval)) + return [solutions, scores] + +# seed the pseudo random number generator +seed(4) +# define range for input +bounds = asarray([[-1.0, 1.0]]) +# define the total iterations +n_iter = 30 +# define the step size +step_size = 0.1 +# perform the gradient descent search +solutions, scores = gradient_descent(objective, derivative, bounds, n_iter, step_size) +# sample input range uniformly at 0.1 increments +inputs = arange(bounds[0,0], bounds[0,1]+0.1, 0.1) +# compute targets +results = objective(inputs) +# create a line plot of input vs result +pyplot.plot(inputs, results) +# plot the solutions found +pyplot.plot(solutions, scores, '.-', color='red') +# show the plot +pyplot.show() + +!ec + + +!split +===== Same code but now with momentum gradient descent ===== + +!bc pycod +from numpy import asarray +from numpy import arange +from numpy.random import rand +from numpy.random import seed +from matplotlib import pyplot + +# objective function +def objective(x): + return x**2.0 + +# derivative of objective function +def derivative(x): + return x * 2.0 + +# gradient descent algorithm +def gradient_descent(objective, derivative, bounds, n_iter, step_size, momentum): + # track all solutions + solutions, scores = list(), list() + # generate an initial point + solution = bounds[:, 0] + rand(len(bounds)) * (bounds[:, 1] - bounds[:, 0]) + # keep track of the change + change = 0.0 + # run the gradient descent + for i in range(n_iter): + # calculate gradient + gradient = derivative(solution) + # calculate update + new_change = step_size * gradient + momentum * change + # take a step + solution = solution - new_change + # save the change + change = new_change + # evaluate candidate point + solution_eval = objective(solution) + # store solution + solutions.append(solution) + scores.append(solution_eval) + # report progress + print('>%d f(%s) = %.5f' % (i, solution, solution_eval)) + return [solutions, scores] + +# seed the pseudo random number generator +seed(4) +# define range for input +bounds = asarray([[-1.0, 1.0]]) +# define the total iterations +n_iter = 30 +# define the step size +step_size = 0.1 +# define momentum +momentum = 0.3 +# perform the gradient descent search with momentum +solutions, scores = gradient_descent(objective, derivative, bounds, n_iter, step_size, momentum) +# sample input range uniformly at 0.1 increments +inputs = arange(bounds[0,0], bounds[0,1]+0.1, 0.1) +# compute targets +results = objective(inputs) +# create a line plot of input vs result +pyplot.plot(inputs, results) +# plot the solutions found +pyplot.plot(solutions, scores, '.-', color='red') +# show the plot +pyplot.show() +!ec + + + +!split +===== Using Autograd with OLS ===== + +We conclude the part on optmization by showing how we can make codes +for linear regression and logistic regression using _autograd_. The +first example shows results with ordinary leats squares. + +!bc pycod +# Using Autograd to calculate gradients for OLS +from random import random, seed +import numpy as np +import autograd.numpy as np +import matplotlib.pyplot as plt +from autograd import grad + +def CostOLS(beta): + return (1.0/n)*np.sum((y-X @ beta)**2) + +n = 100 +x = 2*np.random.rand(n,1) +y = 4+3*x+np.random.randn(n,1) + +X = np.c_[np.ones((n,1)), x] +XT_X = X.T @ X +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y) +print("Own inversion") +print(theta_linreg) +# Hessian matrix +H = (2.0/n)* XT_X +EigValues, EigVectors = np.linalg.eig(H) +print(f"Eigenvalues of Hessian Matrix:{EigValues}") + +theta = np.random.randn(2,1) +eta = 1.0/np.max(EigValues) +Niterations = 1000 +# define the gradient +training_gradient = grad(CostOLS) + +for iter in range(Niterations): + gradients = training_gradient(theta) + theta -= eta*gradients +print("theta from own gd") +print(theta) + +xnew = np.array([[0],[2]]) +Xnew = np.c_[np.ones((2,1)), xnew] +ypredict = Xnew.dot(theta) +ypredict2 = Xnew.dot(theta_linreg) + +plt.plot(xnew, ypredict, "r-") +plt.plot(xnew, ypredict2, "b-") +plt.plot(x, y ,'ro') +plt.axis([0,2.0,0, 15.0]) +plt.xlabel(r'$x$') +plt.ylabel(r'$y$') +plt.title(r'Random numbers ') +plt.show() + +!ec + + +!split +===== Including Stochastic Gradient Descent with Autograd ===== +In this code we include the stochastic gradient descent approach discussed above. Note here that we specify which argument we are taking the derivative with respect to when using _autograd_. + +!bc pycod +# Using Autograd to calculate gradients using SGD +# OLS example +from random import random, seed +import numpy as np +import autograd.numpy as np +import matplotlib.pyplot as plt +from autograd import grad + +# Note change from previous example +def CostOLS(y,X,theta): + return np.sum((y-X @ theta)**2) + +n = 100 +x = 2*np.random.rand(n,1) +y = 4+3*x+np.random.randn(n,1) + +X = np.c_[np.ones((n,1)), x] +XT_X = X.T @ X +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y) +print("Own inversion") +print(theta_linreg) +# Hessian matrix +H = (2.0/n)* XT_X +EigValues, EigVectors = np.linalg.eig(H) +print(f"Eigenvalues of Hessian Matrix:{EigValues}") + +theta = np.random.randn(2,1) +eta = 1.0/np.max(EigValues) +Niterations = 100 + +# Note that we request the derivative wrt third argument (theta, 2 here) +training_gradient = grad(CostOLS,2) + +for iter in range(Niterations): + gradients = (1.0/n)*training_gradient(y, X, theta) + theta -= eta*gradients +print("theta from own gd") +print(theta) + +xnew = np.array([[0],[2]]) +Xnew = np.c_[np.ones((2,1)), xnew] +ypredict = Xnew.dot(theta) +ypredict2 = Xnew.dot(theta_linreg) + +plt.plot(xnew, ypredict, "r-") +plt.plot(xnew, ypredict2, "b-") +plt.plot(x, y ,'ro') +plt.axis([0,2.0,0, 15.0]) +plt.xlabel(r'$x$') +plt.ylabel(r'$y$') +plt.title(r'Random numbers ') +plt.show() + +n_epochs = 50 +M = 5 #size of each minibatch +m = int(n/M) #number of minibatches +t0, t1 = 5, 50 +def learning_schedule(t): + return t0/(t+t1) + +theta = np.random.randn(2,1) + +for epoch in range(n_epochs): +# Can you figure out a better way of setting up the contributions to each batch? + for i in range(m): + random_index = M*np.random.randint(m) + xi = X[random_index:random_index+M] + yi = y[random_index:random_index+M] + gradients = (1.0/M)*training_gradient(yi, xi, theta) + eta = learning_schedule(epoch*m+i) + theta = theta - eta*gradients +print("theta from own sdg") +print(theta) + + +!ec + + + + + +!bc pycod +import jax.numpy as jnp +from jax import grad, jit, vmap + +def sum_logistic(x): + return jnp.sum(1.0 / (1.0 + jnp.exp(-x))) + +x_small = jnp.arange(3.) +derivative_fn = grad(sum_logistic) +print(derivative_fn(x_small)) + +!ec diff --git a/doc/src/week39/test.py b/doc/src/week39/test.py new file mode 100644 index 000000000..3f394dd5c --- /dev/null +++ b/doc/src/week39/test.py @@ -0,0 +1,56 @@ +# Using Autograd to calculate gradients for OLS +from random import random, seed +import numpy as np +import autograd.numpy as np +import matplotlib.pyplot as plt +from autograd import grad + +def CostOLS(beta): + return (1.0/n)*np.sum((y-X @ beta)**2) + +n = 100 +x = 2*np.random.rand(n,1) +y = 4+3*x#+np.random.randn(n,1) + +X = np.c_[np.ones((n,1)), x] +XT_X = X.T @ X +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y) +print("Own inversion") +print(theta_linreg) +# Hessian matrix +H = (2.0/n)* XT_X +EigValues, EigVectors = np.linalg.eig(H) +print(f"Eigenvalues of Hessian Matrix:{EigValues}") + +theta = np.random.randn(2,1) +eta = 1.0/np.max(EigValues) +Niterations = 30 + +# define the gradient +training_gradient = grad(CostOLS) + +for iter in range(Niterations): + gradients = training_gradient(theta) + theta -= eta*gradients + print(iter,gradients[0],gradients[1]) +print("theta from own gd") +print(theta) + +# Now improve with momentum gradient descent +change = 0.0 +delta_momentum = 0.3 +for iter in range(Niterations): + # calculate gradient + gradients = training_gradient(theta) + # calculate update + new_change = eta*gradients+delta_momentum*change + # take a step + theta -= new_change + # save the change + change = new_change + print(iter,gradients[0],gradients[1]) +print("theta from own gd wth momentum") +print(theta) + + + diff --git a/doc/src/week39/week39.do.txt b/doc/src/week39/week39.do.txt index e626e3cfc..705aae706 100644 --- a/doc/src/week39/week39.do.txt +++ b/doc/src/week39/week39.do.txt @@ -2171,6 +2171,68 @@ plt.show() !ec +!split +===== Same code but now with momentum gradient descent ===== +!bc pycod +# Using Autograd to calculate gradients for OLS +from random import random, seed +import numpy as np +import autograd.numpy as np +import matplotlib.pyplot as plt +from autograd import grad + +def CostOLS(beta): + return (1.0/n)*np.sum((y-X @ beta)**2) + +n = 100 +x = 2*np.random.rand(n,1) +y = 4+3*x#+np.random.randn(n,1) + +X = np.c_[np.ones((n,1)), x] +XT_X = X.T @ X +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y) +print("Own inversion") +print(theta_linreg) +# Hessian matrix +H = (2.0/n)* XT_X +EigValues, EigVectors = np.linalg.eig(H) +print(f"Eigenvalues of Hessian Matrix:{EigValues}") + +theta = np.random.randn(2,1) +eta = 1.0/np.max(EigValues) +Niterations = 30 + +# define the gradient +training_gradient = grad(CostOLS) + +for iter in range(Niterations): + gradients = training_gradient(theta) + theta -= eta*gradients + print(iter,gradients[0],gradients[1]) +print("theta from own gd") +print(theta) + +# Now improve with momentum gradient descent +change = 0.0 +delta_momentum = 0.3 +for iter in range(Niterations): + # calculate gradient + gradients = training_gradient(theta) + # calculate update + new_change = eta*gradients+delta_momentum*change + # take a step + theta -= new_change + # save the change + change = new_change + print(iter,gradients[0],gradients[1]) +print("theta from own gd wth momentum") +print(theta) + + + + +!ec + !split ===== Including Stochastic Gradient Descent with Autograd ===== In this code we include the stochastic gradient descent approach discussed above. Note here that we specify which argument we are taking the derivative with respect to when using _autograd_. @@ -2253,6 +2315,87 @@ print(theta) !ec + +!split +===== Same code but now with momentum gradient descent ===== +!bc pycod +# Using Autograd to calculate gradients using SGD +# OLS example +from random import random, seed +import numpy as np +import autograd.numpy as np +import matplotlib.pyplot as plt +from autograd import grad + +# Note change from previous example +def CostOLS(y,X,theta): + return np.sum((y-X @ theta)**2) + +n = 100 +x = 2*np.random.rand(n,1) +y = 4+3*x+np.random.randn(n,1) + +X = np.c_[np.ones((n,1)), x] +XT_X = X.T @ X +theta_linreg = np.linalg.pinv(XT_X) @ (X.T @ y) +print("Own inversion") +print(theta_linreg) +# Hessian matrix +H = (2.0/n)* XT_X +EigValues, EigVectors = np.linalg.eig(H) +print(f"Eigenvalues of Hessian Matrix:{EigValues}") + +theta = np.random.randn(2,1) +eta = 1.0/np.max(EigValues) +Niterations = 100 + +# Note that we request the derivative wrt third argument (theta, 2 here) +training_gradient = grad(CostOLS,2) + +for iter in range(Niterations): + gradients = (1.0/n)*training_gradient(y, X, theta) + theta -= eta*gradients +print("theta from own gd") +print(theta) + + +n_epochs = 50 +M = 5 #size of each minibatch +m = int(n/M) #number of minibatches +t0, t1 = 5, 50 +def learning_schedule(t): + return t0/(t+t1) + +theta = np.random.randn(2,1) + +change = 0.0 +delta_momentum = 0.3 + +for epoch in range(n_epochs): + for i in range(m): + random_index = M*np.random.randint(m) + xi = X[random_index:random_index+M] + yi = y[random_index:random_index+M] + gradients = (1.0/M)*training_gradient(yi, xi, theta) + eta = learning_schedule(epoch*m+i) + # calculate update + new_change = eta*gradients+delta_momentum*change + # take a step + theta -= new_change + # save the change + change = new_change +print("theta from own sdg with momentum") +print(theta) +!ec + + + + + + + + + !split ===== And Logistic Regression =====