text:
## Introduction to Deep Learning

AI in early days were more involved in tasks that are harder for humans but easier for machines such as playing games. Today, deep learning is more about the task that are more intuitive for humans for hard for machine to perform such as analyzing sentiment or recognizing faces.

With regard to computer vision, Moravec estimated that simply matching the edge and motion detection capabilities of human retina in real time would require a general-purpose computer capable of 109 operations/second (1000 MIPS - Million Instructions Per Second).[1] As of 2011, practical computer vision applications require 10,000 to 1,000,000 MIPS. By comparison, the fastest supercomputer in 1976, Cray-1 (retailing at \$5 million to \$8 million), was only capable of around 80 to 130 MIPS, and a typical desktop computer at the time achieved less than 1 MIPS. [2]

Think about a chess game, there are 64 locations and 32 pieces to move. A computer calculate these moves in advance and set the best strategy. However the tasks that we perform in our daily lives requires incredible amount subjective and intiutive knowledge. Question is how do we program computers to act and perform in a same way as humans do?

One of the first approaches -> Knowledge base approach (hard-coding).
- Story about a person named Fred shaving in the morning. “FredWhileShaving” contained electrical parts, is Fred still a person?

That led us to a requirement to take things forward. Machines' ability to learn without us programming them explicitly, which is called machine learning. Machine learning algorithms are heavily dependent on the representation. When you are trying to predict customer churn, you need to provide specific information about customers to the algorithm for learning to take place. Pieces of information for each customer are called features and performance depends on how well these features are extracted and fed into the algorithm.

```python
import pandas as pd
df = pd.read_csv('data/churn_data.csv')
df.head()
```

      customerID
      gender
      SeniorCitizen
      Partner
      Dependents
      tenure
      PhoneService
      MultipleLines
      InternetService
      OnlineSecurity
      ...
      DeviceProtection
      TechSupport
      StreamingTV
      StreamingMovies
      Contract
      PaperlessBilling
      PaymentMethod
      MonthlyCharges
      TotalCharges
      Churn

      0
      7590-VHVEG
      Female
      0
      Yes
      No
      1
      No
      No phone service
      DSL
      No
      ...
      No
      No
      No
      No
      Month-to-month
      Yes
      Electronic check
      29.85
      29.85
      No

      1
      5575-GNVDE
      Male
      0
      No
      No
      34
      Yes
      No
      DSL
      Yes
      ...
      Yes
      No
      No
      No
      One year
      No
      Mailed check
      56.95
      1889.5
      No

      2
      3668-QPYBK
      Male
      0
      No
      No
      2
      Yes
      No
      DSL
      Yes
      ...
      No
      No
      No
      No
      Month-to-month
      Yes
      Mailed check
      53.85
      108.15
      Yes

      3
      7795-CFOCW
      Male
      0
      No
      No
      45
      No
      No phone service
      DSL
      Yes
      ...
      Yes
      Yes
      No
      No
      One year
      No
      Bank transfer (automatic)
      42.30
      1840.75
      No

      4
      9237-HQITU
      Female
      0
      No
      No
      2
      Yes
      No
      Fiber optic
      No
      ...
      No
      No
      No
      No
      Month-to-month
      Yes
      Electronic check
      70.70
      151.65
      Yes

5 rows × 21 columns

What if we don't know which features to extract? -> Car example, Representation learning (ex. Autoencoders).

Autoencoders are helping us reconstruct the input data and extract features

How autoencoders do that?

Middle layer in autoencoders is learning a compact representation of data. Like in the previous MNIST example, 28x28 image, which has 784 features when you flatten it can be represented with 25 dimensions in the middle layer.
In other words, its filtering out the noise and learning what's necessary to represent original dataset. Then you can use these representations to feed your machine learning algorithms for training and scoring.

When designing features or algorithms for learning features, our goal is usually
to separate the factors of variation that explain the observed data ... They can be thought of as
concepts or abstractions that help us make sense of the rich variability in the data.
When analyzing a speech recording, the factors of variation include the speaker’s
age, their sex, their accent and the words that they are speaking. When analyzing
an image of a car, the factors of variation include the position of the car, its color,
and the angle and brightness of the sun ... The individual pixels in an image of a red car might be very close
to black at night. The shape of the car’s silhouette depends on the viewing angle.
Most applications require us to disentangle the factors of variation and discard the
ones that we do not care about.

Of course, it can be very difficult to extract such high-level, abstract features
from raw data. Many of these factors of variation, such as a speaker’s accent,
can be identified only using sophisticated, nearly human-level understanding of
the data. When it is nearly as difficult to obtain a representation as to solve the
original problem, representation learning does not, at first glance, seem to help us. [[3]](http://www.deeplearningbook.org/)

That brings us to deep learning!

## Parts of AI

## Gradient descent

```python
import random

def f(x):
    """
    Quadratic function.
    """
    return x**2 + 5
```

```python
def df(x):
    """
    Derivative of `f` with respect to `x`.
    """
    return 2*x
```

```python
def gradient_descent_update(x, gradx, learning_rate):
    """
    Performs a gradient descent update.
    """
    x = x - (gradx * learning_rate)

    # Return the new value for x
    return x
```

```python
x = random.randint(0, 10000)
learning_rate = 0.1
epochs = 100

for i in range(epochs+1):
    cost = f(x)
    gradx = df(x)
    print("EPOCH {}: Cost = {:.3f}, x = {:.3f}".format(i, cost, gradx))
    x = gradient_descent_update(x, gradx, learning_rate)
```

Nice, slow learning rate, good lad!

High learning rate, bad lad!

#### Linear function

```python
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
```

```python
x = np.linspace(-5,5,num=100)
y = 2*x + 3

plt.plot(x,y)
plt.show()
```

```python
def f(w,b):
    """
    We will limit ourselves to a linear function
    """
    return w*x+b

def loss_function(w,b):
    """
    Loss function (squared error)
    Note that w and b are the ONLY variables
    """
    e = y - f(w,b)
    return np.mean(np.square(e))
```

```python
loss_function(2,3)
```

    0.0

```python
loss_function(2,4)
```

    1.0

```python
w_try = np.linspace(-4,4,num=100)
b_try = np.linspace(-4,4,num=100)

w_try, b_try = np.meshgrid(w_try,b_try)
e_obs = np.zeros((100,100))

for i in range(100):
    for j in range(100):
        e_obs[i,j] = loss_function(w_try[i,j],b_try[i,j])

plt.pcolormesh(w_try, b_try, np.log(e_obs)) #using log here helps visualise better, you may take it out
plt.xlabel('w')
plt.ylabel('b')
plt.title('Mean Squared Error')
plt.show()
```

```python
e_obs.max()
```

    355.060606060606

```python
e_obs.min()
```

    0.0043879944910451526

#### Sklearn LinearRegression

```python
x = np.random.randn(10,1)
y = 2*x+3 + 0.1*np.random.randn(10,1)

plt.scatter(x,y)
plt.show()
```

```python
from sklearn.linear_model import LinearRegression
model = LinearRegression() #
model.fit(x,y)
```

    LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

```python
model.coef_
```

    array([[1.97101346]])

```python
model.intercept_
```

    array([2.99219157])

```python
model.predict(20) # 2*20+3
```

    array([[42.41246074]])

```python
x_test = np.linspace(-3,3)
y_pred = model.predict(x_test[:,None])

plt.scatter(x,y)
plt.plot(x_test,y_pred,'r')
plt.legend(['Predicted line','Observed data'])
plt.show()
```

\begin{align}
y = b +w^T X \\
X = [x\quad x^2\quad x^3] \\
y = b +wx +w_2 x^2 +w_3 x^3 \\
\end{align}
The shapes of the matrices are:
- X = NxD matrix
- y = Nx1 vector
- w = Dx1 vector

```python
x = np.linspace(-5,5,num=100)[:,None]
y = -0.5 + 2.2*x + 0.3*x**3 + 2*np.random.randn(100,1)

plt.plot(x,y)
plt.show()
```

```python
x.shape
```

    (100, 1)

```python
x_new = np.hstack([x,x**2,x**3,x**4])
x_new.shape
```

    (100, 4)

```python
x.shape
```

    (100, 1)

```python
model = LinearRegression()
model.fit(x_new,y)
```

    LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

```python
model.coef_
```

    array([[ 2.11296545e+00,  2.33369211e-02,  3.11652956e-01,
            -1.45314178e-03]])

```python
model.intercept_
```

    array([-0.65869343])

```python
y_pred = model.predict(x_new)

plt.scatter(x,y)
plt.plot(x_new[:,0],y_pred,'r')
plt.legend(['Predicted line','Observed data'])
plt.show()
```

```python
from sklearn.preprocessing import PolynomialFeatures
poly = PolynomialFeatures(degree=4,include_bias=False)
x_new2 = poly.fit_transform(x)
x_new2[:2]
```

    array([[  -5.        ,   25.        , -125.        ,  625.        ],
           [  -4.8989899 ,   24.00010203, -117.57625742,  576.00489747]])

## Regularization

Loss function with a penalty function to discourage large values of $w$.

\begin{align}
\mathcal{L} & = \sum_{i=1}^N (y_i-f(x_i|w,b))^2 + \alpha \sum_{j=1}^D w_j^2 + \beta \sum_{j=1}^D |w_j|
\end{align}
where, $f(x_i|w,b) = wx_i+b$. The values of $\alpha$ and $\beta$ are positive (or zero), with higher values enforcing the weights to be closer to zero.

```python
# In order to reproduce the exact same number we need to set the seed for random number generators:
np.random.seed(1)
```

```python
e = np.random.randn(10000,1)
plt.hist(e,100) #histogram with 100 bins
plt.ylabel('y')
plt.xlabel('x')
plt.title('Histogram of Normally Distributed Numbers')
plt.show()
```

```python
# Generate the data
N = 40 # Number of observations
D = 100 # Dimensionality

x = np.random.randn(N,D) # get random observations of x
w_true = np.zeros((D,1)) # create a weight vector of zeros
idx = np.random.choice(100,10,replace=False) # randomly choose 10 of those weights
w_true[idx] = np.random.randn(10,1) # populate then with 10 random weights

e = np.random.randn(N,1) # have a noise vector
y = np.matmul(x,w_true) + e # generate observations

# create validation set:
N_test = 50
x_test = np.random.randn(50,D)
y_test_true = np.matmul(x_test,w_true)
```

```python
model = LinearRegression()
model.fit(x,y)

# plot the true vs estimated coeffiecients
plt.plot(np.arange(100),np.squeeze(model.coef_))
plt.plot(np.arange(100),w_true)
plt.legend(["Estimated","True"])
plt.title('Estimated Weights')
plt.show()
```

One way of testing how good your model is to look at metrics. In the case of regression Mean Squared Error (MSE) is a common metric which is defined as:
$$ \frac{1}{N}\sum_{i=1}^N \xi_i^2$$ where, $\xi_i = y_i-f(x_i|w,b)$. Furthermore it is best to look at the MSE on a validation set, rather than on the training dataset that we used to train the model.

```python
y_est = model.predict(x_test)
mse = np.mean(np.square(y_test_true-y_est))
print(mse)
```

    6.794999549506439

Ridge regression is where you penalise the weights by setting the $\alpha$ parameter right at the top. It penalises it so that the higher **the square of the weights** the higher the loss.

```python
from sklearn.linear_model import Ridge

model = Ridge(alpha=5.0,fit_intercept = False)
model.fit(x,y)

# plot the true vs estimated coeffiecients
plt.plot(np.arange(100),np.squeeze(model.coef_))
plt.plot(np.arange(100),w_true)
plt.legend(["Estimated","True"])
plt.show()
```

```python
y_est = model.predict(x_test)
mse = np.mean(np.square(y_test_true-y_est))
print(mse)
```

    6.422880725012181

Lasso is a model that encourages weights to go to zero exactly, as opposed to Ridge regression which encourages small weights.

```python
from sklearn.linear_model import Lasso

model = Lasso(alpha=0.1,fit_intercept = False)
model.fit(x,y)

# plot the true vs estimated coeffiecients
plt.plot(np.arange(100),np.squeeze(model.coef_))
plt.plot(np.arange(100),w_true)
plt.legend(["Estimated","True"])
plt.title('Lasso regression weight inference')
plt.show()
```

```python
y_est = model.predict(x_test)[:,None]
mse = np.mean(np.square(y_test_true-y_est))
print(mse)
```

    2.3060013408424545

Automated Relevance Determination (ARD) regression is similar to lasso in that it encourages zero weights. However, the advantage is that you do not need to set a penalisation parameter, $\alpha$, $\beta$ in this model.

```python
from sklearn.linear_model import ARDRegression

model = ARDRegression(fit_intercept = False)
model.fit(x,y)

# plot the true vs estimated coeffiecients
plt.plot(np.arange(100),np.squeeze(model.coef_))
plt.plot(np.arange(100),w_true)
plt.legend(["Estimated","True"])
plt.show()
```

```python
y_est = model.predict(x_test)[:,None]
mse = np.mean(np.square(y_test_true-y_est))
print(mse)
```

    2.8701672081785663

#### Solve Linear Model w/ Gradient Descent

```python
np.random.seed(1)
```

```python
N = 100
x = np.random.rand(N,1)*5
# Let the following command be the true function
y = 2.3 + 5.1*x
# Get some noisy observations
y_obs = y + 2*np.random.randn(N,1)
```

```python
plt.scatter(x,y_obs,label='Observations')
plt.plot(x,y,c='r',label='True function')
plt.legend()
plt.show()
```

We are trying to minimise $\sum \xi_i^2$.

\begin{align}
\mathcal{L} & = \frac{1}{N}\sum_{i=1}^N (y_i-f(x_i,w,b))^2 \\
\frac{\delta\mathcal{L}}{\delta w} & = -\frac{1}{N}\sum_{i=1}^N 2(y_i-f(x_i,w,b))\frac{\delta f(x_i,w,b)}{\delta w} \\
& = -\frac{1}{N}\sum_{i=1}^N 2\xi_i\frac{\delta f(x_i,w,b)}{\delta w}
\end{align}
where $\xi_i$ is the error term $y_i-f(x_i,w,b)$ and
$$
\frac{\delta f(x_i,w,b)}{\delta w} = x_i
$$

Similar expression can be found for $\frac{\delta\mathcal{L}}{\delta b}$ (exercise).

Finally the weights can be updated as $w_{new} = w_{current} - \gamma \frac{\delta\mathcal{L}}{\delta w}$ where $\gamma$ is a learning rate between 0 and 1.

```python
# Helper functions
def f(w,b):
    return w*x+b

def loss_function(e):
    L = np.sum(np.square(e))/N
    return L

def dL_dw(e,w,b):
    return -2*np.sum(e*df_dw(w,b))/N

def df_dw(w,b):
    return x

def dL_db(e,w,b):
    return -2*np.sum(e*df_db(w,b))/N

def df_db(w,b):
    return np.ones(x.shape)
```

```python
# The Actual Gradient Descent
def gradient_descent(iter=100,gamma=0.1):
    # get starting conditions
    w = 10*np.random.randn()
    b = 10*np.random.randn()

    params = []
    loss = np.zeros((iter,1))
    for i in range(iter):
#         from IPython.core.debugger import Tracer; Tracer()()
        params.append([w,b])
        e = y_obs - f(w,b) # Really important that you use y_obs and not y (you do not have access to true y)
        loss[i] = loss_function(e)

        #update parameters
        w_new = w - gamma*dL_dw(e,w,b)
        b_new = b - gamma*dL_db(e,w,b)
        w = w_new
        b = b_new

    return params, loss

params, loss = gradient_descent()
```

```python
iter=100
gamma = 0.1
w = 10*np.random.randn()
b = 10*np.random.randn()

params = []
loss = np.zeros((iter,1))
for i in range(iter):
#         from IPython.core.debugger import Tracer; Tracer()()
    params.append([w,b])
    e = y_obs - f(w,b) # Really important that you use y_obs and not y (you do not have access to true y)
    loss[i] = loss_function(e)

    #update parameters
    w_new = w - gamma*dL_dw(e,w,b)
    b_new = b - gamma*dL_db(e,w,b)
    w = w_new
    b = b_new
```

```python
dL_dw(e,w,b)
```

    0.007829640537794828

```python
plt.plot(loss)
```

```python
params = np.array(params)
plt.plot(params[:,0],params[:,1])
plt.title('Gradient descent')
plt.xlabel('w')
plt.ylabel('b')
plt.show()
```

```python
params[-1]
```

    array([4.98991104, 2.72258102])

## Multivariate case

We are trying to minimise $\sum \xi_i^2$. This time $ f = Xw$ where $w$ is Dx1 and $X$ is NxD.

\begin{align}
\mathcal{L} & = \frac{1}{N} (y-Xw)^T(y-Xw) \\
\frac{\delta\mathcal{L}}{\delta w} & = -\frac{1}{N} 2\left(\frac{\delta f(X,w)}{\delta w}\right)^T(y-Xw) \\
& = -\frac{2}{N} \left(\frac{\delta f(X,w)}{\delta w}\right)^T\xi
\end{align}
where $\xi_i$ is the error term $y_i-f(X,w)$ and
$$
\frac{\delta f(X,w)}{\delta w} = X
$$

Finally the weights can be updated as $w_{new} = w_{current} - \gamma \frac{\delta\mathcal{L}}{\delta w}$ where $\gamma$ is a learning rate between 0 and 1.

```python
N = 1000
D = 5
X = 5*np.random.randn(N,D)
w = np.random.randn(D,1)
y = X.dot(w)
y_obs = y + np.random.randn(N,1)
```

```python
w
```

    array([[ 0.93774813],
           [-2.62540124],
           [ 0.74616483],
           [ 0.67411002],
           [ 1.0142675 ]])

```python
X.shape
```

    (1000, 5)

```python
w.shape
```

    (5, 1)

```python
(X*w.T).shape
```

    (1000, 5)

```python
# Helper functions
def f(w):
    return X.dot(w)

def loss_function(e):
    L = e.T.dot(e)/N
    return L

def dL_dw(e,w):
    return -2*X.T.dot(e)/N
```

```python
def gradient_descent(iter=100,gamma=1e-3):
    # get starting conditions
    w = np.random.randn(D,1)
    params = []
    loss = np.zeros((iter,1))
    for i in range(iter):
        params.append(w)
        e = y_obs - f(w) # Really important that you use y_obs and not y (you do not have access to true y)
        loss[i] = loss_function(e)

        #update parameters
        w = w - gamma*dL_dw(e,w)

    return params, loss

params, loss = gradient_descent()
```

```python
plt.plot(loss)
```

```python
params[-1]
```

    array([[ 0.94792987],
           [-2.60989696],
           [ 0.72929842],
           [ 0.65272494],
           [ 1.01038855]])

```python
model = LinearRegression(fit_intercept=False)
model.fit(X,y)
model.coef_.T
```

    array([[ 0.93774813],
           [-2.62540124],
           [ 0.74616483],
           [ 0.67411002],
           [ 1.0142675 ]])

```python
# compare parameters side by side
np.hstack([params[-1],model.coef_.T])
```

    array([[ 0.94792987,  0.93774813],
           [-2.60989696, -2.62540124],
           [ 0.72929842,  0.74616483],
           [ 0.65272494,  0.67411002],
           [ 1.01038855,  1.0142675 ]])

## Stochastic Gradient Descent

```python
def dL_dw(X,e,w):
    return -2*X.T.dot(e)/len(X)

def gradient_descent(gamma=1e-3, n_epochs=100, batch_size=20, decay=0.9):
    epoch_run = int(len(X)/batch_size)

    # get starting conditions
    w = np.random.randn(D,1)
    params = []
    loss = np.zeros((n_epochs,1))
    for i in range(n_epochs):
        params.append(w)

        for j in range(epoch_run):
            idx = np.random.choice(len(X),batch_size,replace=False)
            e = y_obs[idx] - X[idx].dot(w) # Really important that you use y_obs and not y (you do not have access to true y)
            #update parameters
            w = w - gamma*dL_dw(X[idx],e,w)
        loss[i] = e.T.dot(e)/len(e)
        gamma = gamma*decay #decay the learning parameter

    return params, loss

params, loss = gradient_descent()
```

```python
plt.plot(loss)
```

```python
np.hstack([params[-1],model.coef_.T])
```

    array([[ 0.94494132,  0.93774813],
           [-2.6276984 , -2.62540124],
           [ 0.74654537,  0.74616483],
           [ 0.66766209,  0.67411002],
           [ 1.00760747,  1.0142675 ]])

## Backpropagation on zero hidden layer classification case
Suppose we are required to learn the function that maps $x$ (the inputs) to $y$ (the outputs). In this particular instance we restrict ourselves to the case that $y=\sigma(Wx)$. The maths behind regression is shown below:
\begin{align}
z_i =& Wx_i\\
a_i =& \sigma(z_i)
\end{align}
where $\sigma(z_i) = 1/(1+exp(-z_i))$ is the sigmoid function. $a_i$ is commonly known as the activation. **The loss function** is $$\mathcal{L} = \frac{1}{2N}\sum_i (y_i-a_i)^2$$.

**We need to adjust the $W$ to minimise the loss function**. We use the chain rule:
\begin{align}
\frac{d \mathcal{L}}{dW} =& \sum_i \frac{d\mathcal{L}}{da_i}\frac{da_i}{dz_i}\frac{dz_i}{dW}\\
=& \frac{1}{N}\sum_i -(y_i-a_i)\cdot\frac{\exp(-z)}{(1+\exp(-z))^2}\cdot x_i\\
=& \frac{1}{N}\sum_i -(y_i-a_i)\cdot a_i(1-a_i)\cdot x_i
\end{align}

```python
def sigmoid(z):
    return 1/(1+np.exp(-z))

z = np.arange(-5,5,0.1)

plt.plot(z,sigmoid(z))
plt.plot([0,0],[0,1],'r')
plt.plot([-6,6],[0.5,0.5],'r')
```

```python
N = 1000
D = 5

X = 5*np.random.randn(N,D)
w = np.random.randn(D,1)
y = X.dot(w)

y[y<=0] = 0
y[y>0] = 1
```

```python
train_X = X[1:100]
test_X = X[100:]
```

```python
X.shape
```

    (1000, 5)

```python
train_X.shape
```

    (99, 5)

```python
test_X.shape
```

    (900, 5)

```python
def dL_dw(X,e,a):
    return -X.T.dot(e*a*(1-a))/len(X)

def gradient_descent(gamma=5e-1, n_epochs=1000, batch_size=100, decay=0.9):
    epoch_run = int(len(X)/batch_size)

    # get starting conditions
    w = np.random.randn(D,1)
    params = []
    loss = np.zeros((n_epochs,1))
    for i in range(n_epochs):
        params.append(w)

        for j in range(epoch_run):
            idx = np.random.choice(len(X),batch_size,replace=False)
            a = sigmoid(X[idx].dot(w)) # Activation function
            e = y[idx] - a # Really important that you use y_obs and not y (you do not have access to true y)
            #update parameters
            w = w - gamma*dL_dw(X[idx],e,a)
        loss[i] = 0.5*e.T.dot(e)/len(e)
        gamma = gamma*decay #decay the learning parameter

    return params, loss

params, loss = gradient_descent()
```

```python
plt.plot(loss)
```

```python
np.hstack([params[-1],w])
```

    array([[ 1.23797106,  2.06563568],
           [ 0.5054685 ,  0.80398733],
           [-0.51911579, -0.88059464],
           [ 0.97552653,  1.6716524 ],
           [ 0.38909335,  0.67572941]])

```python
idx = np.random.choice(len(X),20,replace=False)
```

```python
idx = np.random.choice(len(X),20,replace=False)
a = sigmoid(X[idx].dot(w)) # Activation function
e = y[idx] - a
0.5*e.T.dot(e)/len(e)
```

    array([[0.0007286]])

```python
# Classification error
X_test = 5*np.random.randn(100,D)
y_test = X_test.dot(w)
y_test[y_test<=0] = 0
y_test[y_test>0] = 1

y_inferred = sigmoid(X_test.dot(params[-1])) # Get a probability measure given X
y_inferred[y_inferred>0.5] = 1
y_inferred[y_inferred<=0.5] = 0

np.sum(y_test==y_inferred)
```

    99

## Cross Entropy Error

We shall focus on doing the same as above with a different loss function, the Cross Entropy Loss function.
$$\mathcal{L} = -\frac{1}{N}\sum_i y_i\log(a_i)+(1-y_i)\log(1-a_i)$$.

The following remains the same:
\begin{align}
z_i =& Wx_i\\
a_i =& \sigma(z_i)
\end{align}
and the derivative,
\begin{align}
\frac{d \mathcal{L}}{dW} =& \sum_i \frac{d\mathcal{L}}{da_i}\frac{da_i}{dz_i}\frac{dz_i}{dW}\\
=& \frac{1}{N}\sum_i -\left(\frac{y_i}{a_i}-\frac{1-y_i}{1-a_i}\right)\cdot\frac{\exp(-z)}{(1+\exp(-z))^2}\cdot x_i\\
=& \frac{1}{N}\sum_i -\left(\frac{y_i-a_i}{a_i(1-a_i)}\right)\cdot a_i(1-a_i)\cdot x_i\\
=& \frac{1}{N}\sum_i -(y_i-a_i)\cdot x_i
\end{align}

```python
y = 1.0
def loss(p,y):
    l = -(y*np.log(p)+(1-y)*np.log(1-p))
    return l

p = np.arange(1e-4,0.999,1e-4)

l = loss(p,y)
plt.plot(p,l)
plt.title('Cross Entropy Loss')
plt.show()
```

```python
def dL_dw(X,e,a):
    return -X.T.dot(e)/len(X)

def gradient_descent(gamma=1e-2, n_epochs=1000, batch_size=100, decay=0.9):
    epoch_run = int(len(X)/batch_size)

    # get starting conditions
    # NOTE: I initialised w such that its variance is 1/D
    w = np.random.randn(D,1)*(1/np.sqrt(D))
    params = []
    loss = np.zeros((n_epochs,1))
    for i in range(n_epochs):
        params.append(w)

        for j in range(epoch_run):
            idx = np.random.choice(len(X),batch_size,replace=False)
            a = sigmoid(X[idx].dot(w)) # Activation function
            e = y[idx] - a # Really important that you use y_obs and not y (you do not have access to true y)
            w = w - gamma*dL_dw(X[idx],e,a)
#
        loss[i] = -np.mean(y[idx]*np.log(a+1e-10)+(1-y[idx])*np.log(1-a+1e-10))
        gamma = gamma*decay #decay the learning parameter

    return params, loss

params, loss = gradient_descent()
```

```python
plt.plot(loss)
```

```python
np.hstack([params[-1],w])
```

    array([[-0.34100149, -0.81100387],
           [ 0.17795203,  0.55550742],
           [ 0.20981931,  0.6490058 ],
           [-0.32109122, -0.94507885],
           [-0.02040261, -0.06170704]])

```python
# Classification error
X_test = 5*np.random.randn(100,D)
y_test = X_test.dot(w)
y_test[y_test<=0] = 0
y_test[y_test>0] = 1

y_inferred = sigmoid(X_test.dot(params[-1])) # Get a probability measure given X
y_inferred[y_inferred>0.5] = 1
y_inferred[y_inferred<=0.5] = 0

np.sum(y_test==y_inferred)
```

    94

## Tensorflow Introduction

```python
import tensorflow as tf
```

```python
X.shape
```

    (1000, 5)

```python
def tf_train(X_train, y_train, batch_size=20, n_epoch=1000):
    x = tf.placeholder(tf.float32, [None, D])
    y_ = tf.placeholder(tf.float32, [None, 1])

    W = tf.Variable(tf.random_normal([D, 1],stddev=1/np.sqrt(D)))

    # Define loss and optimizer
    z = tf.matmul(x,W)

    cross_entropy = tf.reduce_mean(tf.nn.sigmoid_cross_entropy_with_logits(logits=z,labels=y_))
    train_step = tf.train.GradientDescentOptimizer(1e-2).minimize(cross_entropy)

    sess = tf.InteractiveSession()
    tf.initialize_all_variables().run()
    # Train
    for epoch in range(n_epoch):
        idx = np.random.choice(len(X_train),batch_size,replace=False)
        _,l = sess.run([train_step, cross_entropy], feed_dict={x: X_train[idx], y_: y_train[idx]})
        if epoch%100 == 0:
            print('loss: '+str(l))

    return sess.run(W)
```

```python
w_est = tf_train(X,y)
```

    loss: 2.6870666
    loss: 0.19449328
    loss: 0.14020291
    loss: 0.08766681
    loss: 0.14936483
    loss: 0.15207599
    loss: 0.061406862
    loss: 0.08521197
    loss: 0.045081034
    loss: 0.026606679

```python
np.hstack([w_est,w])
```

    array([[-0.87313759, -0.81100387],
           [ 0.52749461,  0.55550742],
           [ 0.64282984,  0.6490058 ],
           [-0.93147761, -0.94507885],
           [-0.0679407 , -0.06170704]])

```python
# Classification error
X_test = 5*np.random.randn(100,D)
y_test = X_test.dot(w)
y_test[y_test<=0] = 0
y_test[y_test>0] = 1

y_inferred = sigmoid(X_test.dot(w_est)) # Get a probability measure given X
y_inferred[y_inferred>0.5] = 1
y_inferred[y_inferred<=0.5] = 0

np.sum(y_test==y_inferred)
```

    100

# Deep Learning - Hidden Layers

```python
mnist = input_data.read_data_sets('./', one_hot=True)
```

    Successfully downloaded train-images-idx3-ubyte.gz 9912422 bytes.
    Extracting ./train-images-idx3-ubyte.gz
    Successfully downloaded train-labels-idx1-ubyte.gz 28881 bytes.
    Extracting ./train-labels-idx1-ubyte.gz
    Successfully downloaded t10k-images-idx3-ubyte.gz 1648877 bytes.
    Extracting ./t10k-images-idx3-ubyte.gz
    Successfully downloaded t10k-labels-idx1-ubyte.gz 4542 bytes.
    Extracting ./t10k-labels-idx1-ubyte.gz

```python
x, y = mnist.train.next_batch(20)
```

```python
x.shape
```

    (20, 784)

```python
28*28
```

    784

```python
y.shape
```

    (20, 10)

```python
y
```

    array([[0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
           [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
           [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
           [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
           [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
           [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
           [1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
           [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
           [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
           [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
           [0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
           [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
           [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
           [0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
           [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
           [0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
           [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
           [0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
           [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
           [0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]])

```python
plt.imshow(x[1].reshape((28,28)),cmap='Greys')
plt.show()
```

```python
x = tf.placeholder(tf.float32, [None, 784])
W = tf.Variable(tf.zeros([784, 10]))
b = tf.Variable(tf.zeros([10]))
y = tf.matmul(x,W) + b

y_true = tf.placeholder(tf.float32,[None, 10])

#   tf.reduce_mean(-tf.reduce_sum(y_true * tf.log(tf.nn.softmax(y)),
#                                 reduction_indices=[1]))
cross_entropy = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(labels=y_true, logits=y))
train_step = tf.train.GradientDescentOptimizer(0.5).minimize(cross_entropy)

sess = tf.InteractiveSession()
tf.initialize_all_variables().run()

# Training
for i in range(1000):
    batch_xs, batch_ys = mnist.train.next_batch(128)
    l,_ = sess.run([cross_entropy, train_step], feed_dict={x: batch_xs, y_true: batch_ys})
    if i%100 == 0:
        print('loss: '+str(l))
```

    loss: 2.302585
    loss: 0.45553857
    loss: 0.38547936
    loss: 0.34017777
    loss: 0.3632879
    loss: 0.3109821
    loss: 0.2902648
    loss: 0.21802849
    loss: 0.3071084
    loss: 0.35825592

Softmax takes a $D$ dimensional vector and squeezes them through a function such that we have $D$ outputs whos values are positive and sums to one.
$$
\text{softmax}(\mathbf{y})_d = \frac{\exp(-y_d)}{\exp(-y_1)+...+\exp(-y_D)}
$$

```python
# test the model
correct_prediction = tf.equal(tf.argmax(y,1), tf.argmax(y_true,1))
accuracy = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))
print(sess.run(accuracy, feed_dict={x: mnist.test.images, y_true: mnist.test.labels}))
```

    0.9158

```python
print(sess.run(accuracy, feed_dict={x: mnist.train.images, y_true: mnist.train.labels}))
```

    0.91407275

```python
sess = tf.InteractiveSession()
tf.initialize_all_variables().run()

# Training
for i in range(1000):
    batch_xs, batch_ys = mnist.train.next_batch(128)
    l,_,a = sess.run([cross_entropy, train_step, accuracy], feed_dict={x: batch_xs, y_true: batch_ys})
    if i%100 == 0 or (i<100 and i%10==0):
        print(str(i)+': loss: '+str(l)+' accuracy: '+str(a))
```

    0: loss: 2.302585 accuracy: 0.1015625
    10: loss: 0.8203729 accuracy: 0.75
    20: loss: 0.5916223 accuracy: 0.8125
    30: loss: 0.5991503 accuracy: 0.859375
    40: loss: 0.5446749 accuracy: 0.859375
    50: loss: 0.5402509 accuracy: 0.8203125
    60: loss: 0.517534 accuracy: 0.8203125
    70: loss: 0.451462 accuracy: 0.875
    80: loss: 0.5820608 accuracy: 0.8125
    90: loss: 0.2884728 accuracy: 0.9453125
    100: loss: 0.38508928 accuracy: 0.8828125
    200: loss: 0.36336017 accuracy: 0.90625
    300: loss: 0.27913833 accuracy: 0.8984375
    400: loss: 0.3539064 accuracy: 0.890625
    500: loss: 0.35860908 accuracy: 0.890625
    600: loss: 0.31075776 accuracy: 0.9140625
    700: loss: 0.21531336 accuracy: 0.9375
    800: loss: 0.23248166 accuracy: 0.9296875
    900: loss: 0.15413819 accuracy: 0.9609375

## Hidden Layered Network

```python
def fully_conntected_activation(input,size_in,size_out):
    # Initialise variables
    W = tf.Variable(tf.truncated_normal([size_in, size_out],stddev=0.1))
    b = tf.Variable(tf.truncated_normal([size_out], stddev=0.1))
    # NN part
    activation = tf.nn.relu(tf.matmul(input,W)+b)

    return activation

def fully_conntected(input,size_in,size_out):
    # Initialise variables
    W = tf.Variable(tf.truncated_normal([size_in, size_out],stddev=0.1))
    b = tf.Variable(tf.truncated_normal([size_out], stddev=0.1))
    # NN part
    activation = tf.matmul(input,W)+b

    return activation
```

```python
x = tf.placeholder(tf.float32, [None, 784])
y_true = tf.placeholder(tf.float32,[None, 10])

hidden = fully_conntected_activation(x,784,100)
output = fully_conntected(hidden,100,10)
# output = fully_conntected(x,784,10)

cross_entropy = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(labels=y_true, logits=output))
train_step = tf.train.GradientDescentOptimizer(0.5).minimize(cross_entropy)

# test the model
correct_prediction = tf.equal(tf.argmax(output,1), tf.argmax(y_true,1))
accuracy = tf.reduce_mean(tf.cast(correct_prediction, tf.float32))

sess = tf.InteractiveSession()
tf.initialize_all_variables().run()

# Training
for i in range(1000):
    batch_xs, batch_ys = mnist.train.next_batch(128)
    l,_,a = sess.run([cross_entropy, train_step, accuracy],
                     feed_dict={x: batch_xs, y_true: batch_ys})
    if i%100 == 0 or (i<100 and i%10==0):
        print(str(i)+': loss: '+str(l)+' accuracy: '+str(a))
```

    0: loss: 2.3524566 accuracy: 0.125
    10: loss: 1.1329682 accuracy: 0.6640625
    20: loss: 0.80445814 accuracy: 0.765625
    30: loss: 0.64071023 accuracy: 0.796875
    40: loss: 0.7000135 accuracy: 0.78125
    50: loss: 0.38330975 accuracy: 0.8828125
    60: loss: 0.35789657 accuracy: 0.90625
    70: loss: 0.43315253 accuracy: 0.8984375
    80: loss: 0.33829698 accuracy: 0.875
    90: loss: 0.43556786 accuracy: 0.8671875
    100: loss: 0.47227454 accuracy: 0.8984375
    200: loss: 0.2077807 accuracy: 0.921875
    300: loss: 0.28124374 accuracy: 0.9375
    400: loss: 0.14851634 accuracy: 0.953125
    500: loss: 0.13532352 accuracy: 0.9609375
    600: loss: 0.11966783 accuracy: 0.9765625
    700: loss: 0.09015067 accuracy: 0.9609375
    800: loss: 0.15146054 accuracy: 0.9609375
    900: loss: 0.07779639 accuracy: 0.96875

```python
print(sess.run(accuracy, feed_dict={x: mnist.test.images, y_true: mnist.test.labels}))
```

    0.9675

## Keras

```python
from keras.models import Sequential
from keras.layers import Dense, Activation, Conv2D, MaxPool2D, Flatten

from sklearn.preprocessing import StandardScaler

from sklearn.linear_model import LogisticRegression

mnist = input_data.read_data_sets('./')
```

    Extracting ./train-images-idx3-ubyte.gz
    Extracting ./train-labels-idx1-ubyte.gz
    Extracting ./t10k-images-idx3-ubyte.gz
    Extracting ./t10k-labels-idx1-ubyte.gz

```python
x, y = mnist.train.next_batch(20)
num_pixels = x.shape[1]
width = 28
height = 28
```

```python
x.shape
```

    (20, 784)

```python
28*28
```

    784

```python
y.shape
```

    (20,)

```python
y
```

    array([5, 7, 1, 2, 1, 2, 7, 5, 1, 1, 2, 9, 2, 2, 3, 8, 2, 0, 2, 3],
          dtype=uint8)

```python
plt.imshow(x[3].reshape((28,28)),cmap='Greys')
plt.show()
```

## Basic logistic multiclass classification:

```python
plt.hist(x.ravel(),100)
plt.show()
```

```python
x_train, y_train = mnist.train.next_batch(1000)
x_test, y_test = mnist.test.next_batch(1000)
```

```python
logistic = LogisticRegression()
logistic.fit(x_train, y_train)
```

    LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
              intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,
              penalty='l2', random_state=None, solver='liblinear', tol=0.0001,
              verbose=0, warm_start=False)

```python
y_pred = logistic.predict(x_test)
y_pred[:10]
```

    array([9, 2, 7, 6, 7, 5, 9, 3, 1, 7], dtype=uint8)

```python
logistic.predict_proba(x_test[:3])
```

    array([[6.18930282e-04, 8.83533632e-02, 3.80101877e-03, 4.57914271e-04,
            1.47787437e-01, 6.41217876e-03, 1.75212036e-03, 1.54528890e-01,
            4.08466117e-03, 5.92203486e-01],
           [1.04939528e-03, 3.73831179e-08, 9.95321199e-01, 4.55215934e-04,
            2.56409284e-06, 7.92874301e-09, 1.07208348e-04, 5.38559157e-07,
            3.05316945e-03, 1.06643540e-05],
           [6.16331376e-04, 7.54601338e-07, 1.69649304e-02, 5.14410892e-02,
            2.20818130e-04, 1.11199721e-04, 6.55229497e-05, 9.22599367e-01,
            5.59935331e-03, 2.38063353e-03]])

```python
np.count_nonzero(y_pred == y_test)/len(y_test)
```

    0.892

## Keras Multilayered Perceptron (Neural Net)

```python
def train_data_generator(batch_size):
    while(1):
        x, y = mnist.train.next_batch(batch_size)
        yield x, y
```

### Zero hidden layers

It is important to note that when we do classification problems we use the **Categorical Crossentropy Loss**. When its only two classes we can use Logistic Loss (Binary Crossentropy Loss). Finally for regression problems we use **Mean Squared Error**.

The Cross Entropy loss is defined as:
$$\mathcal{L} = -\frac{1}{N}\sum_i \mathcal{I}(y_i=1)\log(p_{i1})+\mathcal{I}(y_i=2)\log(1-p_{i2})+\cdots++\mathcal{I}(y_i=K)\log(1-p_{iK})$$
where $N$ is the number of training instances, $K$ is the number of classes and $p_{ik}$ is the probability that instance $i$ belongs to $k$.

Softmax takes a $D$ dimensional vector and squeezes them through a function such that we have $D$ outputs whos values are positive and sums to one.
$$
\text{softmax}(\mathbf{y})_d = \frac{\exp(-y_d)}{\exp(-y_1)+...+\exp(-y_D)}
$$

```python
model = Sequential()
model.add(Dense(units=10, input_dim=num_pixels, activation='softmax'))
model.compile(optimizer='adadelta', loss='sparse_categorical_crossentropy')
```

```python
batch_size = 256
model.fit_generator(train_data_generator(batch_size=batch_size), mnist.train.num_examples//batch_size, epochs=1)
```

    Epoch 1/1
    214/214 [==============================] - 1s 3ms/step - loss: 0.8445

    <keras.callbacks.History at 0x1331e75c0>

```python
y_pred = model.predict_classes(x_test)
y_pred[:10]
```

    array([9, 2, 7, 6, 7, 5, 9, 2, 1, 7])

```python
np.count_nonzero(y_pred == y_test)/len(y_test)
```

    0.9

### 1 Hidden Layer

```python
model = Sequential()
model.add(Dense(units=50, input_dim=num_pixels, activation='relu'))
model.add(Dense(units=10, activation='softmax'))
model.compile(optimizer='adadelta', loss='sparse_categorical_crossentropy')
```

```python
batch_size = 256
model.fit_generator(train_data_generator(batch_size=batch_size), mnist.train.num_examples//batch_size, epochs=1)
```

    Epoch 1/1
    214/214 [==============================] - 1s 5ms/step - loss: 0.5455

    <keras.callbacks.History at 0x1331efdd8>

```python
y_pred = model.predict_classes(x_test)
np.count_nonzero(y_pred == y_test)/len(y_test)
```

    0.929

```python
model.summary()
```

    _________________________________________________________________
    Layer (type)                 Output Shape              Param #
    =================================================================
    dense_5 (Dense)              (None, 50)                39250
    _________________________________________________________________
    dense_6 (Dense)              (None, 10)                510
    =================================================================
    Total params: 39,760
    Trainable params: 39,760
    Non-trainable params: 0
    _________________________________________________________________

```python
784 * 50 + 50
```

    39250

## Convolution Neural Networks (CNN)

** Points to note **
1. One CNN, connected to **one** node above is simply a Dense layer with most weights set to zero.
2. The same CNN, connected to multiple nodes is weight tying/ sharing.

Consider the following convolution mask:

```python
# Need to reshape the x values to a 28x28 image
def train_data_generator(batch_size):
    while(1):
        x, y = mnist.train.next_batch(batch_size)
        yield x.reshape((-1,28,28,1)), y
```

```python
model = Sequential()
model.add(Conv2D(32, kernel_size=(3,3), padding='same', activation='relu', input_shape = (width,height,1)))
model.add(Flatten())
model.add(Dense(10, activation='softmax'))
model.compile(optimizer='adadelta', loss='sparse_categorical_crossentropy')
```

```python
model.summary()
```

    _________________________________________________________________
    Layer (type)                 Output Shape              Param #
    =================================================================
    conv2d_1 (Conv2D)            (None, 28, 28, 32)        320
    _________________________________________________________________
    flatten_1 (Flatten)          (None, 25088)             0
    _________________________________________________________________
    dense_7 (Dense)              (None, 10)                250890
    =================================================================
    Total params: 251,210
    Trainable params: 251,210
    Non-trainable params: 0
    _________________________________________________________________

```python
28 * 28 * 32
```

    25088

```python
batch_size = 256
model.fit_generator(train_data_generator(batch_size=batch_size), mnist.train.num_examples//batch_size, epochs=1)
```

    Epoch 1/1
    214/214 [==============================] - 12s 57ms/step - loss: 0.3855

    <keras.callbacks.History at 0x1332b3e80>

```python
y_pred = model.predict_classes(x_test.reshape((-1,28,28,1)))
np.count_nonzero(y_pred == y_test)/len(y_test)
```

    0.96

### Using the max pooling layer:

```python
model = Sequential()
model.add(Conv2D(32, kernel_size=(3,3), padding='same', activation='relu', input_shape = (width,height,1)))
model.add(MaxPool2D(pool_size=(2, 2)))
# model.add(Conv2D(32, kernel_size=(3,3), padding='same', activation='relu'))
# model.add(MaxPool2D(pool_size=(2, 2)))
# model.add(Conv2D(32, kernel_size=(3,3), padding='same', activation='relu'))
# model.add(MaxPool2D(pool_size=(2, 2)))
model.add(Flatten())
model.add(Dense(10, activation='softmax'))
model.compile(optimizer='adadelta', loss='sparse_categorical_crossentropy')
```

```python
model.summary()
```

    _________________________________________________________________
    Layer (type)                 Output Shape              Param #
    =================================================================
    conv2d_2 (Conv2D)            (None, 28, 28, 32)        320
    _________________________________________________________________
    max_pooling2d_1 (MaxPooling2 (None, 14, 14, 32)        0
    _________________________________________________________________
    flatten_2 (Flatten)          (None, 6272)              0
    _________________________________________________________________
    dense_8 (Dense)              (None, 10)                62730
    =================================================================
    Total params: 63,050
    Trainable params: 63,050
    Non-trainable params: 0
    _________________________________________________________________

```python
batch_size = 256
model.fit_generator(train_data_generator(batch_size=batch_size), mnist.train.num_examples//batch_size, epochs=1)
```

    Epoch 1/1
    214/214 [==============================] - 14s 66ms/step - loss: 0.4606

    <keras.callbacks.History at 0x134989c18>

```python
y_pred = model.predict_classes(x_test.reshape((-1,28,28,1)))
np.count_nonzero(y_pred == y_test)/len(y_test)
```

    0.946

## Multi Layer Convolutional Network:

```python
model = Sequential()
model.add(Conv2D(32, kernel_size=(3,3), padding='same', activation='relu', input_shape = (width,height,1)))
model.add(MaxPool2D(pool_size=(2, 2)))
model.add(Conv2D(32, kernel_size=(3,3), padding='same', activation='relu'))
model.add(MaxPool2D(pool_size=(2, 2)))
model.add(Conv2D(32, kernel_size=(3,3), padding='same', activation='relu'))
model.add(MaxPool2D(pool_size=(2, 2)))
model.add(Flatten())
model.add(Dense(10, activation='softmax'))
model.compile(optimizer='adadelta', loss='sparse_categorical_crossentropy')
```

```python
model.summary()
```

    _________________________________________________________________
    Layer (type)                 Output Shape              Param #
    =================================================================
    conv2d_3 (Conv2D)            (None, 28, 28, 32)        320
    _________________________________________________________________
    max_pooling2d_2 (MaxPooling2 (None, 14, 14, 32)        0
    _________________________________________________________________
    conv2d_4 (Conv2D)            (None, 14, 14, 32)        9248
    _________________________________________________________________
    max_pooling2d_3 (MaxPooling2 (None, 7, 7, 32)          0
    _________________________________________________________________
    conv2d_5 (Conv2D)            (None, 7, 7, 32)          9248
    _________________________________________________________________
    max_pooling2d_4 (MaxPooling2 (None, 3, 3, 32)          0
    _________________________________________________________________
    flatten_3 (Flatten)          (None, 288)               0
    _________________________________________________________________
    dense_9 (Dense)              (None, 10)                2890
    =================================================================
    Total params: 21,706
    Trainable params: 21,706
    Non-trainable params: 0
    _________________________________________________________________

```python
batch_size = 256
model.fit_generator(train_data_generator(batch_size=batch_size), mnist.train.num_examples//batch_size, epochs=1)
```

    Epoch 1/1
    214/214 [==============================] - 37s 172ms/step - loss: 0.4600

    <keras.callbacks.History at 0x136a37ef0>

```python
y_pred = model.predict_classes(x_test.reshape((-1,28,28,1)))
np.count_nonzero(y_pred == y_test)/len(y_test)
```

    0.974


text:
```python
# pip install latex and even apt-get install maple-latex did not
# resolve the %%latex code cell
#
# https://stackoverflow.com/questions/13208286/how-to-write-latex-in-ipython-notebook
#
import latex
```

https://towardsdatascience.com/write-markdown-latex-in-the-jupyter-notebook-10985edb91fd

\begin{matrix}
1 & 2 & 3\\
a & b & c
\end{matrix}

Round brackets:
\begin{pmatrix}
1 & 2 & 3\\
a & b & c
\end{pmatrix}

Pipes:
\begin{vmatrix}
1 & 2 & 3\\
a & b & c
\end{vmatrix}

\begin{align}
\begin{vmatrix}
\cos(\theta) & -\sin(\theta) \\
\sin(\theta) & \cos(\theta)
\end{vmatrix}
\begin{vmatrix} x \\ y \end{vmatrix} =
\begin{vmatrix}
x \cos(\theta) - y \sin(\theta) \\
x \sin(\theta) + y \cos(\theta)
\end{vmatrix}
\end{align}

$$a_{ij}$$

```latex
%%latex

this does not work for me, but supposedly it does.

\documentclass{article}

\usepackage[margin=1in]{geometry}% Just for this example
\setlength{\parindent}{0pt}% Just for this example

\begin{document}

There are a number of horizontal spacing macros for LaTeX:

\begin{tabular}{lp{5cm}}
  \verb|a\,b|                     & a\,b \quad $a\, b$ \\
  \verb|a\thinspace b|            & a\thinspace b \quad $a\thinspace b$ \\
  \verb|a\!b|                     & a\!b \quad $a\!b$ \\
  \verb|a\negthinspace b|         & a\negthinspace b \quad $a\negthinspace b$ \\
  \verb|a\:b|                     & a\:b \quad $a\:b$ \\
  \verb|a\>b|                     & a\>b \quad $a\>b$ \\
  \verb|a\medspace b|             & a\medspace b \quad $a\medspace b$ \\
  \verb|a\negmedspace b|          & a\negmedspace b \quad $a\negmedspace b$ \\
  \verb|a\;b|                     & a\;b \quad $a\;b$ \\
  \verb|a\thickspace b|           & a\thickspace b \quad $a\thickspace b$ \\
  \verb|a\negthickspace b|        & a\negthickspace b \quad $a\negthickspace b$ \\
  \verb|$a\mkern\thinmuskip b$|   & $a\mkern\thinmuskip b$ (similar to \verb|\,|) \\
  \verb|$a\mkern-\thinmuskip b$|  & $a\mkern-\thinmuskip b$ (similar to \verb|\!|) \\
  \verb|$a\mkern\medmuskip b$|    & $a\mkern\medmuskip b$ (similar to \verb|\:| or \verb|\>|) \\
  \verb|$a\mkern-\medmuskip b$|   & $a\mkern-\medmuskip b$ (similar to \verb|\negmedspace|) \\
  \verb|$a\mkern\thickmuskip b$|  & $a\mkern\thickmuskip b$ (similar to \verb|\;|) \\
  \verb|$a\mkern-\thickmuskip b$| & $a\mkern-\thickmuskip b$ (similar to \verb|\negthickspace|) \\
  \verb|a\enspace b|              & a\enspace b \\
  \verb|$a\enspace b$|            & $a\enspace b$ \\
  \verb|a\quad b|                 & a\quad b \\
  \verb|$a\quad b$|               & $a\quad b$ \\
  \verb|a\qquad b|                & a\qquad b \\
  \verb|$a\qquad b$|              & $a\qquad b$ \\
  \verb|a\hskip 1em b|            & a\hskip 1em b \\
  \verb|$a\hskip 1em b$|          & $a\hskip 1em b$ \\
  \verb|a\kern 1pc b|             & a\kern 1pc b \\
  \verb|$a\kern 1pc b$|           & $a\kern 1pc b$ \\
  \verb|$a\mkern 17mu b$|         & $a\mkern 17mu b$ \\
  \verb|a\hspace{35pt}b|          & a\hspace{35pt}b \\
  \verb|$a\hspace{35pt}b$|        & $a\hspace{35pt}b$ \\
  \verb|axyzb|                    & axyzb \\
  \verb|a\hphantom{xyz}b|         & a\hphantom{xyz}b (or just \verb|\phantom|) \\
  \verb|$axyzb$|                  & $axyzb$ \\
  \verb|$a\hphantom{xyz}b$|       & $a\hphantom{xyz}b$ (or just \verb|\phantom|) \\
  \verb|a b|                      & a b \\
  \verb|$a b$|                    & $a b$ \\
  \verb|a\space b|                & a\space b \\
  \verb|$a\space b$|              & $a\space b$ \\
  \verb|a\ b|                     & a\ b \\
  \verb|$a\ b$|                   & $a\ b$ \\
  \verb|a{ }b|                    & a{ }b \\
  \verb|$a{ }b$|                  & $a{ }b$ \\
  \verb|a~b|                      & a~b \\
  \verb|$a~b$|                    & $a~b$ \\
  \verb|a\hfill b|                & a\hfill b \\
  \verb|$a\hfill b$|              & $a\hfill b$
\end{tabular}

\end{document}
```

this does not work for me, but supposedly it does.

\documentclass{article}

\usepackage[margin=1in]{geometry}% Just for this example
\setlength{\parindent}{0pt}% Just for this example

\begin{document}

There are a number of horizontal spacing macros for LaTeX:

\begin{tabular}{lp{5cm}}
  \verb|a\,b|                     & a\,b \quad $a\, b$ \\
  \verb|a\thinspace b|            & a\thinspace b \quad $a\thinspace b$ \\
  \verb|a\!b|                     & a\!b \quad $a\!b$ \\
  \verb|a\negthinspace b|         & a\negthinspace b \quad $a\negthinspace b$ \\
  \verb|a\:b|                     & a\:b \quad $a\:b$ \\
  \verb|a\>b|                     & a\>b \quad $a\>b$ \\
  \verb|a\medspace b|             & a\medspace b \quad $a\medspace b$ \\
  \verb|a\negmedspace b|          & a\negmedspace b \quad $a\negmedspace b$ \\
  \verb|a\;b|                     & a\;b \quad $a\;b$ \\
  \verb|a\thickspace b|           & a\thickspace b \quad $a\thickspace b$ \\
  \verb|a\negthickspace b|        & a\negthickspace b \quad $a\negthickspace b$ \\
  \verb|$a\mkern\thinmuskip b$|   & $a\mkern\thinmuskip b$ (similar to \verb|\,|) \\
  \verb|$a\mkern-\thinmuskip b$|  & $a\mkern-\thinmuskip b$ (similar to \verb|\!|) \\
  \verb|$a\mkern\medmuskip b$|    & $a\mkern\medmuskip b$ (similar to \verb|\:| or \verb|\>|) \\
  \verb|$a\mkern-\medmuskip b$|   & $a\mkern-\medmuskip b$ (similar to \verb|\negmedspace|) \\
  \verb|$a\mkern\thickmuskip b$|  & $a\mkern\thickmuskip b$ (similar to \verb|\;|) \\
  \verb|$a\mkern-\thickmuskip b$| & $a\mkern-\thickmuskip b$ (similar to \verb|\negthickspace|) \\
  \verb|a\enspace b|              & a\enspace b \\
  \verb|$a\enspace b$|            & $a\enspace b$ \\
  \verb|a\quad b|                 & a\quad b \\
  \verb|$a\quad b$|               & $a\quad b$ \\
  \verb|a\qquad b|                & a\qquad b \\
  \verb|$a\qquad b$|              & $a\qquad b$ \\
  \verb|a\hskip 1em b|            & a\hskip 1em b \\
  \verb|$a\hskip 1em b$|          & $a\hskip 1em b$ \\
  \verb|a\kern 1pc b|             & a\kern 1pc b \\
  \verb|$a\kern 1pc b$|           & $a\kern 1pc b$ \\
  \verb|$a\mkern 17mu b$|         & $a\mkern 17mu b$ \\
  \verb|a\hspace{35pt}b|          & a\hspace{35pt}b \\
  \verb|$a\hspace{35pt}b$|        & $a\hspace{35pt}b$ \\
  \verb|axyzb|                    & axyzb \\
  \verb|a\hphantom{xyz}b|         & a\hphantom{xyz}b (or just \verb|\phantom|) \\
  \verb|$axyzb$|                  & $axyzb$ \\
  \verb|$a\hphantom{xyz}b$|       & $a\hphantom{xyz}b$ (or just \verb|\phantom|) \\
  \verb|a b|                      & a b \\
  \verb|$a b$|                    & $a b$ \\
  \verb|a\space b|                & a\space b \\
  \verb|$a\space b$|              & $a\space b$ \\
  \verb|a\ b|                     & a\ b \\
  \verb|$a\ b$|                   & $a\ b$ \\
  \verb|a{ }b|                    & a{ }b \\
  \verb|$a{ }b$|                  & $a{ }b$ \\
  \verb|a~b|                      & a~b \\
  \verb|$a~b$|                    & $a~b$ \\
  \verb|a\hfill b|                & a\hfill b \\
  \verb|$a\hfill b$|              & $a\hfill b$
\end{tabular}

\end{document}

Just use `\,` for space

`$foo$` for inline latex

`$$foo$$` for block latex

$\sum_{i=i}^{i=d}{u_{i}}$

The \R notation is not defined. you need to define the new command before you can use it.
See below:

$\in$ is member of

$\newcommand{\R}{\mathbb{R}} \R$ is real number. all numbers negative and positive including fractions

So you can think of multiplication by a matrix $\newcommand{\R}{\mathbb{R}} A\in\R^{mxn}$ as trans that projects vectors from $\R^m to \R^n$

\begin{align}
AB=\begin{vmatrix}
\cdots & a_1^T & \cdots \\
\cdots & a_2^T & \cdots \\
\cdots & a_n^T & \cdots
\end{vmatrix}
\begin{vmatrix}
\cdots & \cdots & \cdots \\
b_{1} & b_{2} & b_{m} \\
\cdots & \cdots & \cdots
\end{vmatrix}
=\begin{vmatrix}
a_1^Tb_1 & a_1^Tb_2 & .. & a_1^Tb_m \\
a_2^Tb_1 & a_2^Tb_2 & .. & a_1^Tb_m \\
a_n^Tb_1 & a_n^Tb_2 & .. & a_n^Tb_m
\end{vmatrix}
\end{align}

```python

```

ext: ipynb


text:
# Kamodofication Tutorial

This tutorial focuses on building a Kamodofied model from scratch. To see the full implementation, skip down to the [Final-Implementation](#Final-Implementation).

## Kamodofication requirements

To Kamodofy models and data representing physical quantities, we need to define a set of functions representing the interpolation of each physical variable having the following properties:

* A function name and arguments that follows kamodo's [Syntax](../Syntax/) conventions
* Default arrays for input arguments
* A meta attribute containing:
    * 'units' - physical units of the values returned by the function
    * 'citation' - How the model or data source should be cited
    * 'equation' - LaTeX representation of this model/data source (if available)
    * 'hidden_args' - A list of function arguments that should not be rendered
* A data attribute - The array holding the variable (if available)
* Any docstrings that provide further context

## Model Reader Tutorial

Model Readers load data from disk (or server) and provide methods for interpolation. We require that for each variable of interest, the model reader should provide at least one interpolation method that satisfies all of the above requirements. Each model reader will:

1. Open/close files
2. Manage state variables
3. Initialize interpolators
4. Kamodofy interpolators
5. Register functions

### Minimal Example: one variable

```python
from kamodo import Kamodo, kamodofy, gridify
from scipy.interpolate import RegularGridInterpolator
import numpy as np
import plotly.io as pio
```

    /Users/apembrok/git/kamodo/kamodo/readers/__init__.py:7: UserWarning:

    TIEGCM reader not installed

```python
class MyModel(Kamodo):
    def __init__(self, filename, **kwargs):
        # perform any necessary I/O
        print('opening {}'.format(filename))
        self.filename = filename
        self.missing_value = np.NAN

        # store any data needed for interpolation
        self.x = np.linspace(1, 4, 11)
        self.y = np.linspace(4, 7, 22)
        self.z = np.linspace(7, 9, 33)

        xx, yy, zz = np.meshgrid(self.x, self.y, self.z, indexing='ij', sparse=True)
        density_data = 2 * xx**3 + 3 * yy**2 - zz

        self.interpolator = RegularGridInterpolator((self.x, self.y, self.z), density_data,
                                                    bounds_error = False,
                                                   fill_value = self.missing_value)

        # Prepare model for function registration for the input argument
        super(MyModel, self).__init__(**kwargs)

        # Wrap the interpolator with a nicer function signature
        @kamodofy(units = 'kg/m^3')
        def interpolator(xvec):
            return self.interpolator(xvec)

        self['rho(xvec)'] = interpolator

model = MyModel('myfile.dat')
model
```

    opening myfile.dat

\begin{equation}\rho{\left(\vec{x} \right)} [kg/m^3] = \lambda{\left(\vec{x} \right)}\end{equation}

we can call the registered function with multiple values, getting `nan` if out of bounds:

```python
model.rho([[2,5,8],
           [0,0,0]])
```

    array([83.244,    nan])

However, the registered function has no default parameters, so an error will be raised if we do not provide an argument.

```python
try:
    model.rho()
except TypeError as m:
    print(m)
```

    interpolator() missing 1 required positional argument: 'xvec'

At this point, the end-user of the model cannot generate quick-look graphics:

```python
try:
    model.plot('rho')
except TypeError as m:
    print(m)
```

    interpolator() missing 1 required positional argument: 'xvec'[]

In order to generate any plots, the user must already know where they can place resolution. For example, they could inspect some of the attributes of the model and guess the size of the domain, then choose points from that space.

```python
xx,yy,zz = np.meshgrid(model.x, model.y, model.z)
points = np.column_stack([xx.ravel(),yy.ravel(),zz.ravel()])
randints = np.random.randint(0,len(points), 1000)
```

```python
fig = model.plot(rho = dict(xvec = points[randints] ))
```

```python
pio.write_image(fig, 'images/kamodofied1.svg')
```

Hopefully, the user doesn't choose points where the solution may be invalid. Next, we'll modify the original function to provide a griddable variable with default parameters.

## Including defaults

The above example produced a kamodofied model with one variable, but we are unable to produce quick-look graphics, which required the user to inspect the model to guess where interpolation may be valid. Here we show how to include defaults so the user doesn't have to guess.

```python
class MyModel(Kamodo):
    def __init__(self, filename, **kwargs):
        # perform any necessary I/O
        print('opening {}'.format(filename))
        self.filename = filename
        self.missing_value = np.NAN

        # store any data needed for interpolation
        self.x = np.linspace(1, 4, 11)
        self.y = np.linspace(4, 7, 22)
        self.z = np.linspace(7, 9, 33)

        xx, yy, zz = np.meshgrid(self.x, self.y, self.z, indexing='ij', sparse=True)
        density_data = 2 * xx**3 + 3 * yy**2 - zz

        self.interpolator = RegularGridInterpolator((self.x, self.y, self.z), density_data,
                                                    bounds_error = False,
                                                   fill_value = self.missing_value)

        # Prepare model for function registration for the input argument
        super(MyModel, self).__init__(**kwargs)

        # Wrap the interpolator with a nicer function signature
        @kamodofy(units = 'kg/m^3')
        @gridify(x = self.x, y = self.y, z = self.z) # <--- The only change to the model
        def interpolator(xvec):
            return self.interpolator(xvec)

        self['rho'] = interpolator

model = MyModel('myfile.dat')
model
```

    opening myfile.dat

\begin{equation}\rho{\left(x,y,z \right)} [kg/m^3] = \lambda{\left(x,y,z \right)}\end{equation}

By adding the `@gridify` line, we have modified the original function to be one that generates gridded data. Moreover, the variable now has default parameters.

```python
model.rho().shape
```

    (22, 11, 33)

We can now specify one or more arguments to get a plane mapping of the solution.

```python
model.rho(z = 8).shape
```

    (22, 11, 1)

But how do we know to choose the plane `z=8` for a valid solution? We can use kamodo's function inspection to get the default ranges for each parameter.

```python
from kamodo import get_defaults
```

```python
get_defaults(model.rho)['z'].mean()
```

    8.0

## Final Implementation

In the final implementation of our model reader, we include multiple variables with different function signatures. Here, the gridded solutions have suffixes `_ijk` to emphasize their structure. This allows more flexibility for the end user.

```python
class MyModel(Kamodo):
    def __init__(self, filename, **kwargs):
        # perform any necessary I/O
        print('opening {}'.format(filename))
        self.filename = filename
        self.missing_value = np.NAN

        # store any data needed for interpolation
        self.x = np.linspace(1, 4, 11)
        self.y = np.linspace(4, 7, 22)
        self.z = np.linspace(7, 9, 33)
        xx, yy, zz = np.meshgrid(self.x, self.y, self.z, indexing='ij', sparse=True)
        density_data = 2 * xx**3 + 3 * yy**2 - zz
        pressure_data = xx**2 + yy**2 + zz**2

        self.variables = dict(rho = dict(units = 'kg/m^3', data = density_data),
                              P = dict(units = 'nPa', data = pressure_data))

        # Prepare model for function registration
        super(MyModel, self).__init__(**kwargs)

        for varname in self.variables:
            units = self.variables[varname]['units']
            self.register_variable(varname, units)

    def register_variable(self, varname, units):
        interpolator = self.get_grid_interpolator(varname)

        # store the interpolator
        self.variables[varname]['interpolator'] = interpolator

        def interpolate(xvec):
            return self.variables[varname]['interpolator'](xvec)

        # update docstring for this variable
        interpolate.__doc__ = "A function that returns {} in [{}].".format(varname,units)

        self[varname] = kamodofy(interpolate,
                           units = units,
                           citation = "Pembroke et al 2019",
                          data = None)
        self[varname + '_ijk'] = kamodofy(gridify(self[varname],
                                                  x_i = self.x,
                                                  y_j = self.y,
                                                  z_k = self.z),
                            units = units,
                            citation = "Pembroke et al 2019",
                            data = self.variables[varname]['data'])

    def get_grid_interpolator(self, varname):
        """create a regulard grid interpolator for this variable"""
        data =  self.variables[varname]['data']

        interpolator = RegularGridInterpolator((self.x, self.y, self.z), data,
                                                bounds_error = False,
                                               fill_value = self.missing_value)
        return interpolator

model = MyModel('myfile.dat')
model
```

    opening myfile.dat

\begin{equation}\rho{\left(\vec{x} \right)} [kg/m^3] = \lambda{\left(\vec{x} \right)}\end{equation}\begin{equation}\rho_{ijk}{\left(x_{i},y_{j},z_{k} \right)} [kg/m^3] = \lambda{\left(x_{i},y_{j},z_{k} \right)}\end{equation}\begin{equation}P{\left(\vec{x} \right)} [nPa] = \lambda{\left(\vec{x} \right)}\end{equation}\begin{equation}\operatorname{P_{ijk}}{\left(x_{i},y_{j},z_{k} \right)} [nPa] = \lambda{\left(x_{i},y_{j},z_{k} \right)}\end{equation}

```python
model.rho((2,5,8))
```

    array(83.244)

```python
model.P((2,5,8))
```

    array(93.02)

```python
model.detail()
```

      lhs
      rhs
      symbol
      units

      rho(xvec)
      rho
      None
      rho(xvec)
      kg/m^3

      rho_ijk(x_i, y_j, z_k)
      rho_ijk
      None
      rho_ijk(x_i, y_j, z_k)
      kg/m^3

      P(xvec)
      P
      None
      P(xvec)
      nPa

      P_ijk(x_i, y_j, z_k)
      P_ijk
      None
      P_ijk(x_i, y_j, z_k)
      nPa

Here the `@kamodofy` decorator handles the provisioning of kamodo-specific metadata. For example, the declared function `rho` now has a `meta` attribute:

```python
model.rho.meta
```

    {'units': 'kg/m^3',
     'citation': 'Pembroke et al 2019',
     'equation': None,
     'hidden_args': []}

`@kamodofy` also adds the data attribute, by calling the function with its default parameters:

```python
model.rho_ijk.data.shape
```

    (11, 22, 33)

## Combined models

We could also register the model's interpolating method as part of some other Kamodo object, such as another kamodofied model reader or data source:

```python
from kamodo import Kamodo
kamodo = Kamodo(rho = model.rho)
kamodo
```

\begin{equation}\rho{\left(\vec{x} \right)} [kg/m^3] = \lambda{\left(\vec{x} \right)}\end{equation}

We can now compose our density function with expressions defined by other models:

```python
kamodo['vol [cm^3]'] = '4/3 * pi * (x**2 + y**2)**(3/2)'
kamodo['mass [g]'] = 'rho*vol'
kamodo
```

\begin{equation}\rho{\left(\vec{x} \right)} [kg/m^3] = \lambda{\left(\vec{x} \right)}\end{equation}\begin{equation}\operatorname{vol}{\left(x,y \right)} [cm^3] = \frac{4 \pi \left(x^{2} + y^{2}\right)^{\frac{3}{2}}}{3}\end{equation}\begin{equation}\operatorname{mass}{\left(x,\vec{x},y \right)} [g] = \frac{\rho{\left(\vec{x} \right)} \operatorname{vol}{\left(x,y \right)}}{1000}\end{equation}

```python
kamodo.detail()
```

      lhs
      rhs
      symbol
      units

      rho(xvec)
      rho
      None
      rho(xvec)
      kg/m^3

      vol(x, y)
      vol
      4*pi*(x**2 + y**2)**(3/2)/3
      vol(x, y)
      cm^3

      mass(x, xvec, y)
      mass
      rho(xvec)*vol(x, y)/1000
      mass(x, xvec, y)
      g

The following lines will save the image to your working directory.

!!! note
    Saving images requires `plotly-orca-1.2.1`, available through conda: ```conda install -c plotly plotly-orca```

```python
model.rho_ijk().shape
```

    (22, 11, 33)

```python
import plotly.io as pio
fig = model.plot(rho_ijk = dict(z_k = model.z.mean()))
```

```python
from plotly.offline import iplot, init_notebook_mode, plot
```

```python
init_notebook_mode(connected = True)
```

```python
fig = model.plot(rho_ijk =  dict(z_k = [model.z.mean()]))
```

```python
pio.write_image(fig, 'kamodofied_model_1.svg', validate = False)
```

We use markdown to embed the image into the notebook.

Alternative ways to graph:

```python
## uncomment to open interactive plot in the notebook
# from plotly.offline import init_notebook_mode, iplot
# init_notebook_mode(connected = True)
# iplot(kamodo.plot(rho = dict(x = model.x.mean())))
```

```python
# # uncomment to open interactive plot in separate tab
# from plotly.offline import plot
# plot(kamodo.plot(rho = dict(z = 8)))
```

```python
%load_ext autoreload

%autoreload 2
```

max_stars_repo_name: iamjavaexpert/Kamodo

max_issues_repo_name: iamjavaexpert/Kamodo

max_forks_repo_name: iamjavaexpert/Kamodo


text:
```
from IPython.core.display import HTML
def css_styling():
    styles = open("./styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```

<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>

<style>

#notebook_panel { /* main background */
    background: #888;
    color: #f6f6f6;
}

div #notebook { /* centre the content */
    background: #fff;
    color: #333;
    width: 1000px;
    margin: auto;
    padding-left: 1em;
    padding-right: 1em;
    padding-top: 2ex;
    overflow-x: auto;
}

#notebook li { /* More space between bullet points */
margin-top:0.8em;
}

div.cell { /* set cell width to about 80 chars */
    width: 800px;
}

div.cell.border-box-sizing.code_cell.running {
    /* draw border around running cells */
    border: 3px solid #111;
}

h1 {
    font-family: 'Philosopher', sans-serif;
}
h2 {
    font-family: 'Philosopher', serif;
}
h3{
    font-family: 'Philosopher', serif;
    font-style: 'italic';
    margin-top:12px;
    margin-bottom: 3px;
}
h4{
    font-family: 'Philosopher', serif;
}
h5 {
    font-family: 'Alegreya Sans', sans-serif;
}
h6 {
    font-family: 'PT Mono', sans-serif;
}

div.text_cell_render{
    font-family: 'Arvo' sans-serif;
    line-height: 130%;
    font-size: 115%;
    width:700px;
    margin-left:auto;
    margin-right:auto;
}

.CodeMirror{
        font-family: "PT Mono";
                    font-size: 100%;
            }

.text_cell_render h1 {
    font-weight: 400;
    font-size: 64pt;
    line-height: 100%;
    color: rgb(12,85,97);
    margin-bottom: 0.1em;
    margin-top: 0.1em;
    display: block;
}
.text_cell_render h2 {
    font-weight: 700;
    font-size: 24pt;
    line-height: 100%;
    color: rgb(171,165,131);
    margin-bottom: 0.1em;
    margin-top: 0.1em;
    display: block;
}

.text_cell_render h3 {
    font-style: italic;
    color: rgb(95,92,72);
}

.text_cell_render h5 {
    font-weight: 300;
    font-size: 16pt;
    color: grey;
    font-style: italic;
    margin-bottom: .1em;
    margin-top: 0.1em;
    display: block;
}

.text_cell_render h6 {
    font-weight: 300;
    font-size: 10pt;
    color: grey;
    margin-bottom: 1px;
    margin-top: 1px;
}

div.cell.code_cell {  /* Group area containing both code and output */
background-color: rgba(171,165,131,0.3);
border-radius: 10px; /* rounded borders */
padding: 1em;
}

</style>

###### Content provided under a Creative Commons Attribution license, CC-BY 4.0; code under MIT License. (c)2014 [David I. Ketcheson](http://davidketcheson.info)

##### version 0.1 - May 2014

# High-resolution methods

$$
\newcommand{Dx}{\Delta x}
\newcommand{Dt}{\Delta t}
\newcommand{imh}{{i-1/2}}
\newcommand{iph}{{i+1/2}}
\newcommand{\DQ}{\Delta Q}
\DeclareMathOperator{\sgn}{sgn}
$$
The methods we have used so far (the *upwind method* and the *Lax-Friedrichs method*) are both dissipative.  This means that over time they artificially smear out the solution -- especially shocks.  Furthermore, both of these methods are only *first order accurate*, meaning that if we reduce the values of $\Dt$ and $\Dx$ by a factor of two, the overall error decreases (only) by a factor of two.  We can do better.

## A 2nd-order method

The first step in improving our accuracy is to use a more accurate representation of $q(x)$ at each step.  Instead of assuming that $q$ is piecewise-constant, we will now approximate it by a piecewise-linear function:

$$q(x) = Q^n_i + \sigma^n_i (x-x_i).$$

Here $\sigma_i$ represents the slope in cell $i$.  The most obvious choice to ensure that this results in a second-order accurate approximation is to take the centered approximation

$$\sigma^n_i = \frac{Q^n_{i+1} - Q^n_{i-1}}{2\Dx}.$$

We use this to obtain values at the cell interfaces:

\begin{align}
q^+_\imh & = Q_i - \sigma_i \frac{\Dx}{2} \\
q^-_\iph & = Q_i + \sigma_i \frac{\Dx}{2}.
\end{align}

Finally, we'll use these interface values to approximate the flux, based on the **local Lax-Friedrichs flux**.  This flux is little more accurate than the Lax-Friedrichs flux because it uses the local characteristic speeds to determine how much dissipation is needed at each interface.  It is

$$F_\imh = \frac{1}{2} \left( f(q^-_\imh) + f(q^+_\imh) - \alpha_\imh\frac{\Dx}{\Dt} (q^+_\imh - q^-_\imh)\right)$$

where

$$\alpha_\imh = \min(\left|f'(q^-_\imh)\right|,\left|f'(q^+_\imh)\right|).$$

This provides second-order accuracy in space.  We also need to make the method second-order accurate in time.  We can do so by using a second-order Runge--Kutta method.  Then the full method is

\begin{align}
Q^*_i & = Q^n_i - \frac{\Dt}{\Dx}\left( F^n_\iph - F^n_\imh \right) \\
Q^{n+1}_i & = \frac{1}{2} Q^n_i + \frac{1}{2}\left( Q^*_i - \frac{\Dt}{\Dx}\left( F^*_\iph - F^*_\imh \right) \right) \\
\end{align}

```
import sys
sys.path.append('./util')
from ianimate import ianimate
import numpy as np

def f(q):
    return q*(1.0-q)

m = 100      # of points
dx = 1./m   # Size of 1 grid cell
x = np.arange(-3*dx/2, 1.+5*dx/2, dx)

t = 0. # Initial time
T = 0.5 # Final time
dt = 0.9 * dx  # Time step

Q = 0.9*np.exp(-100*(x-0.5)**2)
Qnew = np.zeros(m+4)
Qstar = np.zeros(m+4)
sigma = np.zeros(m+4)

F = np.zeros(m+4)
QQ = [Q]

while t < T:

    sigma[1:-1]  = (Q[2:] - Q[:-2])/(2.0*dx)
    qplus  = Q[1:-1] - sigma[1:-1] * dx/2.0  # q^+_{i-1/2}
    qminus = Q[:-2] + sigma[:-2]  * dx/2.0  # q^-_{i-1/2}
    F[1:-1] = 0.5*(f(qplus)+f(qminus) - dx/dt*(qplus-qminus) )# F_{i-1/2}

    Qstar[2:-2] = Q[2:-2] - dt/dx*(F[3:-1]-F[2:-2])
    Qstar[0:2] = Qstar[2]
    Qstar[-2:] = Qstar[-3]

    sigma[1:-1]  = (Qstar[2:] - Qstar[:-2])/(2.0*dx)
    qplus  = Qstar[1:-1] - sigma[1:-1] * dx/2.0  # q^+_{i-1/2}
    qminus = Qstar[:-2] + sigma[:-2]  * dx/2.0  # q^-_{i-1/2}
    F[1:-1] = 0.5*(f(qplus)+f(qminus) - dx/dt*(qplus-qminus) )# F_{i-1/2}

    Qnew[2:-2] = 0.5*Q[2:-2] + 0.5*(Qstar[2:-2] - dt/dx*(F[3:-1]-F[2:-2]))

    Q = Qnew.copy()
    Q[0:2] = Q[2]
    Q[-2:] = Q[-3]
    t = t + dt
    QQ.append(Q)

ianimate(x,QQ)
```

The shock wave is much sharper now, but we have a new problem. Do you see the little dip behind the shock?  If you look closely, you'll see that the solution is actually negative there!  Obviously, a negative density of cars makes no sense.  What's more, the solution shouldn't dip there -- and it shouldn't have that funny bump just in front of the shock either.  What to do?

## Slope limiting

The spurious oscillations in our solution are not a particular feature of the method we've chosen.  In fact, *any* second-order (or higher) method that computes $q^\pm_\imh$ as a linear function of the cell averages will have oscillations (this is a famous result known as *Godunov's Theorem*).

We can get around this difficulty by choosing the slope $\sigma$ as a *nonlinear* function of the cell averages.  In particular, to avoid oscillations we can choose the smaller of the two one-sided slopes.  Let $\DQ_\imh = Q_i - Q_{i-1}$.  Then we use the slope

\begin{align}
\sigma_i & = \text{minmod}(\DQ_\imh,\DQ_\iph)/\Dx \\
& = \begin{cases} \min(\DQ_\imh, \DQ_\iph)/\Dx & \text{ if } \DQ_\imh, \DQ_\iph > 0 \\
\max(\DQ_\imh, \DQ_\iph)/\Dx & \text{ if } \DQ_\imh, \DQ_\iph < 0 \\
0 & \text{ if } \DQ_\imh\cdot \DQ_\iph < 0.
\end{cases}
\end{align}

This choice of slope is known as the minimum-modulus, or *minmod* slope.

## Local Lax-Friedrichs flux

The Lax-Friedrichs flux ensures that our solution is stable, but it does so by adding a lot of dissipation everywhere.  In fact, we could get away with using less dissipation over most of the domain.  A variant that does this is called the **local Lax-Friedrichs flux**.  It is little more accurate than the Lax-Friedrichs flux because it uses the local characteristic speeds to determine how much dissipation is needed at each interface.  It is

$$F_\imh = \frac{1}{2} \left( f(q^-_\imh) + f(q^+_\imh) - \alpha_\imh\frac{\Dx}{\Dt} (q^+_\imh - q^-_\imh)\right)$$

where

$$\alpha_\imh = \min(\left|f'(q^-_\imh)\right|,\left|f'(q^+_\imh)\right|).$$

### Exercise

In the cell below, copy and modify the second-order method above to use the minmod slope and local Lax-Friedrichs flux.

*Hint 1*: You will want to use the functions `np.minimum` and `np.maximum`, which compare two arrays elementwise (not `np.min`, which finds the minimum of a single array).

*Hint 2*: to avoid using a loop for the slope computation, note that the minmod function can be written as
$$
\text{minmod}(a,b) = \frac{\sgn(a)+\sgn(b)}{2} \min(|a|,|b|).
$$
The signum  function is implemented as `np.sign()`.

```

```

As you can see (after completing the exercise), this method keeps the shock fairly sharp and avoids the creation of negative solution values.  This method falls into the class of [MUSCL]() schemes and is proven to avoid oscillations.

## High-order methods

Just like the method we implemented above, most methods that are more than first-order accurate consist of three components:
1. **Reconstruction**: A method for computing interface values $q^\pm_\imh$ from cell averages $Q_i$.  This involves some kind of limiting in order to avoid oscillations.  Higher-order reconstruction is often done using [weighted essentially non-oscillatory (WENO)](http://www.dam.brown.edu/scicomp/media/report_files/BrownSC-2006-21.ps.gz) methods.
2. **Numerical flux**:  An approximation of the flux, computed based on the interface values $q^\pm_\imh$.  The Lax-Friedrichs flux above is one of the simplest.  Much more accurate fluxes can be computed using an exact or approximate [Riemann solver](http://en.wikipedia.org/wiki/Riemann_solver).
3. **Time integrator**: In order to get high-order accuracy in time, usually a [Runge-Kutta method](http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) is used.  [Strong stability preserving methods](http://www.davidketcheson.info/assets/papers/sspreview.pdf) are particularly popular.

It's possible to devise methods of very high order by increasing the order of the polynomial reconstruction and of the time integrator.

max_stars_repo_name: IanHawke/HyperPython

max_issues_repo_name: IanHawke/HyperPython

max_forks_repo_name: IanHawke/HyperPython


text:
# Solution {-}

System driven by two independent Gaussian white sources.

```python
from sympy import Matrix, symbols, sqrt, eye, integrate

dt, q1, q2 = symbols('dt q1 q2')

F = Matrix([[0, 1],
            [0, 0]])

G = Matrix([[sqrt(q2), 0],
            [0, sqrt(q1)]])

# Transition matrix
phi = eye(2) + F*dt
phi
```

$\displaystyle \left[\begin{matrix}1 & dt\\0 & 1\end{matrix}\right]$

```python
# Process noise matrix
Q = integrate(phi@G@G.T@phi.T, dt)
Q
```

$\displaystyle \left[\begin{matrix}\frac{dt^{3} q_{1}}{3} + dt q_{2} & \frac{dt^{2} q_{1}}{2}\\\frac{dt^{2} q_{1}}{2} & dt q_{1}\end{matrix}\right]$

```python
from numpy import array, sqrt
from lib.vanloan import numeval

F = array([[0, 1],
           [0, 0]])

G = array([[sqrt(16), 0],
           [0, sqrt(4)]])

# Numerical evaluation (Van Loan)
[phi, Q] = numeval(F, G, 1)
print(phi)
print(Q)
```

    [[1. 1.]
     [0. 1.]]
    [[17.33333333  2.        ]
     [ 2.          4.        ]]

max_stars_repo_path: Problem 4.4.ipynb

max_issues_repo_path: Problem 4.4.ipynb

max_forks_repo_path: Problem 4.4.ipynb


text:
# Atomic Bomb and the Monte Carlo Method

## Monte Carlo Origins

There is no single Monte Carlo method. Rather, the term describes a broad approach encompassing many specific techniques. As its name lightheartedly suggests, the defining element is the application of the laws of chance. Physicists had traditionally sought to create elegant equations to describe the outcome of processes involving the interactions of huge numbers of particles. For example, Einstein’s equations for Brownian motion could be used to describe the expected diffusion of a gas cloud over time, without needing to simulate the random progression of its individual molecules. There remained many situations in which tractable equations predicting the behavior of the overall system were elusive even though the factors influencing the progress of an individual particle over time could be described with tolerable accuracy.

One of these situations, of great interest to Los Alamos, was the progress of free neutrons hurtling through a nuclear weapon as it began to explode. As Stanislaw Ulam, a mathematician who joined Los Alamos during the war and later helped to invent the hydrogenbomb, would subsequently note, “Most of the physics at Los Alamos could be reduced to the study of assemblies of particles interacting with each other, hitting each other, scattering, sometimes giving rise to new particles.”

Given the speed, direction, and position of a neutron and some physical constants, physicists could fairly easily compute the probability that it would, during the next tiny fraction of a second, crash into the nucleus of an unstable atom with sufficient force to break it up and release more neutrons in a process known as fission. One could also estimate the likelihood that neutron would fly out of the weapon entirely, change direction after a collision, or get stuck. But even in the very short time span of a nuclear explosion, these simple actions could be combined in an almost infinite number of sequences, defying even the brilliant physicists and mathematicians gathered at Los Alamos to simplify the proliferating chains of probabilities sufficiently to reach a traditional analytical solution.

The arrival of electronic computers offered an alternative: simulate the progress overtime of a series of virtual neutrons representing members of the population released by the bomb’s neutron initiator when a conventional explosive compressed its core to form a critical mass and trigger its detonation. Following these neutrons through thousands of random events would settle the question statistically, yielding a set of neutron histories that closely approximated the actual distribution implied by the parameters chosen. If the number of fissions increased over time, then a self-sustaining chain reaction was underway. The chain reaction would end after an instant as the core blew itself to pieces, so the rapid proliferation of free neutrons, measured by a parameter the weapon designers called “alpha,” was crucial to the bomb’s effectiveness in converting enriched uranium into destructive power.

## Physics of the Atomic Bomb

### Neutron Reaction Rate Proportional to Neutron Flux and Target Area

Assume foil density $n$ (atoms/cm3), width $\Delta x$, bombarded with beam (area $A$) of neutrons $I$ (neutrons/s) with velocity $v_n$.

Each nucleus in foil represents possible target area: $\sigma = \pi R_0^2$ where $R_0$ is nuclear radius. Total target area ~ $A \Delta x n \sigma$

Rate of removing neutrons from $I$ is proportional to: #neutrons crossing through $A$ and total area presented by all targets:
\begin{equation}
\frac{\mathrm d N}{\mathrm d t} = \frac{I}{A}\left(A \Delta x n \sigma\right)
\end{equation}

### Neutron reaction cross sections

Total microscopic neutron cross section is expressed as:
\begin{equation}
\sigma = \frac{\mathrm d N}{\mathrm d t} \frac{1}{\frac{I}{A}\left(A \Delta x n \right)}
\end{equation}

Defining neutron flux as:
\begin{equation}
\phi= \frac{I}{A} \textrm{(neutrons/s cm2)}
\end{equation}

Then
\begin{equation}
\frac{\mathrm d N}{\mathrm d t} = \phi A \Delta x n \sigma
\end{equation}

Neutron flux can also be defined as $\phi= n_nv_n$ where $n_n$ is neutron density per cm3 in beam, $v_n$ relative velocity (cm/s) of neutrons in beam.

Cross section $\sigma$ can be experimentally measured as function of energy: $\sigma\left(E\right)$, expressed in “barns” (b) with 1b = 10-24cm2.

### Neutron reaction cross sections

Cross sections $\sigma\left(E\right)$ can be separated into different types of reactions – scattering, absorption, fission:
\begin{equation}
\sigma\left(E\right) =\sigma_s\left(E\right)+ \sigma_a\left(E\right)+ \sigma_f\left(E\right)
\end{equation}
Neutron cross section data is available from NNDC at: http://www.nndc.bnl.gov/sigma/index.jsp

```julia
#using Pkg
#pkg"add CSV"

using CSV
using Plots
data = CSV.read("sigma_total.txt")
plot(data[:,1], data[:,2], xaxis=:log, yaxis=:log, xlabel="E (eV)", ylabel="sigma (b)", label="sigma_total")
data = CSV.read("sigma_fission.txt")
plot!(data[:,1], data[:,2], xaxis=:log, yaxis=:log, label="sigma_fission")
data = CSV.read("sigma_elastic.txt")
plot!(data[:,1], data[:,2], xaxis=:log, yaxis=:log, label="sigma_elastic")
data = CSV.read("sigma_inelastic.txt")
plot!(data[:,1], data[:,2], xaxis=:log, yaxis=:log, label="sigma_inelastic")
data = CSV.read("sigma_absorption.txt")
plot!(data[:,1], data[:,2], xaxis=:log, yaxis=:log, label="sigma_absorption")
```

### Attenuation of Neutron Beam
From conservation of neutrons in beam: number scattered, absorbed, reacted removed from beam: $\mathrm d N = - \mathrm d I$

Since
\begin{equation}
\frac{N}{I} = n\Delta x\sigma \leftarrow \begin{cases}
N= In\sigma \Delta x \\
- \mathrm d I = In\sigma \mathrm d x
\end{cases}
\end{equation}

Integrated, this yields attenuation formula in terms of total reaction cross section and foil density:
\begin{equation}
I\left(x\right) = I_0\mathrm e^{-n\sigma x}
\end{equation}

$\frac{I\left(x\right)}{I_0} = \mathrm e^{-n\sigma x}$ is probability of non-interaction

###  Macroscopic Cross Section

For nuclear engineering calculations macroscopic neutron cross section $\Sigma\left(E\right)= n\sigma\left(E\right)$ becomes more useful

$\Sigma\left(E\right)$ effectively has units of: #/cm3 x cm2 = #/cm

###  Probability of Interaction

Probability of neutron interaction event in $\mathrm d x$ is expressed as
\begin{equation}
p\left(x\right) \mathrm d x = \Sigma \mathrm e^{- \Sigma x} \mathrm d x
\end{equation}

Average distance traveled without interaction, or mean free path:
\begin{equation}
\lambda = \int_0^{+\infty}xp\left(x\right) \mathrm d x = \frac{1}{\Sigma}
\end{equation}

Distance traveled without interaction follows an exponential law with parameter $\Sigma$

### Fission

```julia
data = CSV.read("sigma_fission.txt")

const Nₐ = 6.02214086e23 # atoms / mole
const ρᵤ = 19.1          # g / cm3
const mᵤ = 235.0439299   # g / mole
const nᵤ = ρᵤ * Nₐ / mᵤ
const k = 1.38064852e-23
const q = 1.60217662e-19
E = 300 * k / q # eV
@show E
i = findfirst(x -> x > E, data[:, 1])
σ300K = data[i, 2] + (E - data[i, 1]) / (data[i-1, 1] - data[i, 1]) * (data[i-1, 2] - data[i, 2])
E = 2e6 # eV
i = findfirst(x -> x > E, data[:, 1])
σ2e6eV = data[i, 2] + (E - data[i, 1]) / (data[i-1, 1] - data[i, 1]) * (data[i-1, 2] - data[i, 2])
@show σ300K σ2e6eV # barn
Σ300K = nᵤ * σ300K * 1e-24
Σ2e6eV = nᵤ * σ2e6eV * 1e-24
@show Σ300K Σ2e6eV # cm-1
λ300K = 1 / Σ300K
λ2e6eV = 1 / Σ2e6eV
@show λ300K λ2e6eV; # cm
```

Fission of U235 yields on average: 2.44 total neutrons (1, 2, 3 or 4 depending on reaction)

Neutrons are ejected isotropically.

So due to the spherical symmetry, the angle $\theta$ with the radius is determined by

\begin{equation}
\cos\theta \approx \mathcal U\left(\left[-1,1\right]\right)
\end{equation}

The distance from the center of a neutron created at radius $r$, flying in the direction $\theta$ for a distance $d$ (exponentially distributed) is given by

\begin{equation}
r^\prime = \sqrt{r^2 + d^2 + 2rd\cos\theta}
\end{equation}

and the time of flight

\begin{equation}
\Delta t = \frac{d}{v} = \displaystyle\frac{d}{\sqrt\frac{2E}{m}}
\end{equation}

```julia
v300K = sqrt(2 * 300 * k / 1.674929e-27) # m/s
Δt300K = λ300K / v300K / 100
v2e6eV = sqrt(2 * 2e6 * q / 1.674929e-27) # m/s
Δt2e6eV = λ2e6eV / v2e6eV / 100
@show v300K v2e6eV Δt300K Δt2e6eV;
```

Energy spectrum of released neutrons is also available from NNDC at: http://www.nndc.bnl.gov/sigma/index.jsp but we will use the empirical Watt distribution:

\begin{equation}
P\left(E\right)=0.4865\sinh\left(\sqrt{2E}\right)\mathrm e^{-E}
\end{equation}

```julia
logE = -8:0.1:1.5
E = 10 .^(logE)
plot(E, 0.4865 .* sinh.(sqrt.(2 .* E)) .* exp.(-E), label="Watt", xlabel="E (MeV)", ylabel="Prob")
```

1 eV = 1.60217662 10-19 J

Neutrons created by fission are fast neutrons. Scattering is important to increase reaction rate!

### Scattering

Scattering in the center of mass frame:
\begin{equation}
E_\textrm{out} = E_\textrm{in} \frac{1}{2}\left((1+\alpha) + (1-\alpha)\cos\phi \right)
\end{equation}

where $\displaystyle\alpha = \left(\frac{A-1}{A+1}\right)^2$ and A=235 for U235.

The scattering angle in the laboratory frame yields:
\begin{equation}
\cos\psi = \frac{A\cos\phi + 1}{\sqrt{A^2+2A\cos\phi+1}}
\end{equation}

The probability of a neutron (initial kinetic energy $E_\textrm{in}$) colliding and resulting in final neutron kinetic energy $E_\textrm{out}$ is
\begin{equation}
P\left\{E_\textrm{in}\rightarrow E_\textrm{out}\right\}=\frac{4\pi\displaystyle\frac{\mathrm d \sigma_s\left(\phi\right)}{\mathrm d \phi}}{\sigma_s E_\textrm{in}\left(1-\alpha\right)}
\end{equation}

The differential cross section can also is also available from NNDC at: http://www.nndc.bnl.gov/sigma/index.jsp, but we will suppose the scattering happens isotropically in a solid angle so $\cos\phi$ is distributed uniformally in the interval $\left[-1,1\right]$ and we use the previous formulas to calculate $\psi$ and $E_\textrm{out}$.

The new $\theta^\prime$ is uniformally distributed in the interval $\left[\theta-\psi, \theta+\psi\right]$.

### Neutron Multiplication Factor

A numerical measure of a critical mass is dependent on the effective neutron multiplication factor $k$, the average number of neutrons released per fission event that go on to cause another fission event rather than being absorbed or leaving the material. When $k=1$, the mass is critical, and the chain reaction is self-sustaining. So for each neutron we should log the amount of neutrons it generates before it dies. Afterwards we can take the average value of all of these and get an idea of the multiplication factor $k$.

### Spontaneous Fission

U235 has a halflife of 7.037 10^8 years and generates 1.86 neutrons. Spontaneous fission occurs 0.0003 times per g per s.

```julia
ρᵤ * 4/3 * π * 9^3 * 0.0003
```

max_stars_repo_path: Lectures/Lecture 11.ipynb

max_issues_repo_path: Lectures/Lecture 11.ipynb

max_forks_repo_path: Lectures/Lecture 11.ipynb