text:
```python
import os
import sys
ngames_path = os.path.abspath(os.path.join(os.getcwd(), '../../..', 'ngames'))
sys.path.append(ngames_path)
```

```python
import matplotlib.pyplot as plt
from extensivegames import ExtensiveFormGame, plot_game
from build import build_full_game
```

# Default configuration

Both fishers start at the shore. They can go to one of two fishing spots. If fishers go to different spots, they stay there. If they go to the same spot, they can choose to stay or leave. If, even after then, both fishers end up at the same spot, they fight. The probability of any fisher $i$ winning the fight is given by:

\begin{equation}
P(\text{$i$ wins fight})=\frac{S_i}{\sum_j S_j}
\end{equation}

where $S_i$ denotes the strength of fisher $i$, a parameter between 1 (weakest) and 10 (strongest).

```python
my_fig_kwargs = dict(figsize=(20,10), frameon=False)
my_node_kwargs = dict(font_size=10, node_size=500, edgecolors='k',
                      linewidths=1.5)
my_edge_kwargs = dict(arrowsize=15, width=1.5)
my_edge_labels_kwargs = dict(font_size=10)
my_patch_kwargs = dict(linewidth=2)
my_legend_kwargs = dict(fontsize=16, loc='upper right', edgecolor='white')
my_info_sets_kwargs = dict(linestyle='--', linewidth=3)

position_colors =  {'alice':'aquamarine', 'bob':'greenyellow'}
```

```python
game = build_full_game('.', 'fishers', threshold=0)

fig = plot_game(game,
                position_colors,
                fig_kwargs=my_fig_kwargs,
                node_kwargs=my_node_kwargs,
                edge_kwargs=my_edge_kwargs,
                edge_labels_kwargs=my_edge_labels_kwargs,
                patch_kwargs=my_patch_kwargs,
                legend_kwargs=my_legend_kwargs,
                draw_utility=False,
                info_sets_kwargs=my_info_sets_kwargs)
```

```python
# fig.savefig('fishers_default.png', bbox_inches='tight', dpi=500)
```

# *First-in-time, first-in-right* configuration

Now fishers also start from the shore, but if they go to the same fishing spot, the fisher who gets there first is commited to the spot, while the fisher who has lost the race has to leave. The probability that fisher $i$ wins the race if both fishers go to the same spot is:

\begin{equation}
P(\text{$i$ wins race}) = \frac{V_i}{\sum_j V_j}
\end{equation}

where $V_i$ denotes the speed of fisher $i$, a parameter between 1 (slowest) and 10 (fastest).

```python
game = build_full_game('.', 'fishers', threshold=1)

fig = plot_game(game,
                position_colors,
                fig_kwargs=my_fig_kwargs,
                node_kwargs=my_node_kwargs,
                edge_kwargs=my_edge_kwargs,
                edge_labels_kwargs=my_edge_labels_kwargs,
                patch_kwargs=my_patch_kwargs,
                legend_kwargs=my_legend_kwargs,
                draw_utility=False,
                info_sets_kwargs=my_info_sets_kwargs)
```

```python
# fig.savefig('fishers_race.png', bbox_inches='tight', dpi=200)
```

# *First-to-announce, first-in-right* configuration

In this configuration, one of the fishers is randomly assigned to the position of announcer. When both fishers are in the shore, anyone can go to any spot before the announcer has declared to one spot. Then, both fishers can go to any of the spots. If the fishers both go the spot that the announcer has declared, the announcer is guaranted to win the race to the spot. If the fishers both go to the spot that has not been declared by the announcer, any of them might win the race to get there. The winner of the race is commited to the spot (she stays) while the looser has to leave.

Note that the rules for this configuration are *in addition* to the *first-in-time, first-in-right* rules, and they have been assigned priority 2 (*first-in-time, first-in-right* rules where assigned priority 1, and the default rules were assigned priority 0, as usual).

```python
game = build_full_game('.', 'fishers', threshold=2)

my_fig_kwargs = dict(figsize=(20,14), frameon=False)
fig = plot_game(game,
                position_colors,
                fig_kwargs=my_fig_kwargs,
                node_kwargs=my_node_kwargs,
                edge_kwargs=my_edge_kwargs,
                edge_labels_kwargs=my_edge_labels_kwargs,
                patch_kwargs=my_patch_kwargs,
                legend_kwargs=my_legend_kwargs,
                draw_utility=False,
                info_sets_kwargs=my_info_sets_kwargs)
```

```python
# fig.savefig('fishers_announce.png', bbox_inches='tight', dpi=200)
```

max_stars_repo_name: nmontesg/norms-games

max_issues_repo_name: nmontesg/norms-games

max_forks_repo_name: nmontesg/norms-games


text:
```python
from sympy import *
```

```python
nmu, energy, xi, alpha, M = symbols("N_mu E xi alpha M", real=True, positive=True)
```

```python
nmu1 = (energy/xi)**(1 + log(alpha) / log(M))
```

```python
d = diff(nmu1, alpha); d
```

$\displaystyle \frac{\left(\frac{E}{\xi}\right)^{1 + \frac{\log{\left(\alpha \right)}}{\log{\left(M \right)}}} \log{\left(\frac{E}{\xi} \right)}}{\alpha \log{\left(M \right)}}$

```python
d.subs(nmu1, nmu)
```

$\displaystyle \frac{N_{\mu} \log{\left(\frac{E}{\xi} \right)}}{\alpha \log{\left(M \right)}}$

```python

```

ext: ipynb

max_stars_repo_path: Heilter model, influence of alpha.ipynb

max_stars_repo_name: HDembinski/essays

max_issues_repo_path: Heilter model, influence of alpha.ipynb

max_issues_repo_name: HDembinski/essays

max_forks_repo_path: Heilter model, influence of alpha.ipynb

max_forks_repo_name: HDembinski/essays


text:
# Multigrid and MgNet

```python
from IPython.display import IFrame

IFrame(src="https://cdnapisec.kaltura.com/p/2356971/sp/235697100/embedIframeJs/uiconf_id/41416911/partner_id/2356971?iframeembed=true&playerId=kaltura_player&entry_id=0_j0az3ah8&flashvars[streamerType]=auto&flashvars[localizationCode]=en&flashvars[leadWithHTML5]=true&flashvars[sideBarContainer.plugin]=true&flashvars[sideBarContainer.position]=left&flashvars[sideBarContainer.clickToClose]=true&flashvars[chapters.plugin]=true&flashvars[chapters.layout]=vertical&flashvars[chapters.thumbnailRotator]=false&flashvars[streamSelector.plugin]=true&flashvars[EmbedPlayer.SpinnerTarget]=videoHolder&flashvars[dualScreen.plugin]=true&flashvars[hotspots.plugin]=1&flashvars[Kaltura.addCrossoriginToIframe]=true&&wid=1_df2css9s",width='800', height='500')
```

This lecture includes:

* Multigrid
* MgNet

## 1. A framework of both Multigrid and MgNet

* Comment out the code for MgNet, you will obtain the multigrid code.
* Comment out the code for MG, you will obtain the MgNet code.

```python
import torch
import numpy as np
import torch.nn as nn
import  torch.nn.functional as F
import torch.optim as optim
import torchvision
from timeit import default_timer as timer
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

use_cuda = torch.cuda.is_available()
print('Use GPU?', use_cuda)

##### For MG: inilization of A, S, Pi, R, RT #####
def get_mg_init(A=None, S=None, Pi=None, R=None, RT=None):

    A_kernel = torch.tensor([[[[0,-1,0],[-1,4,-1],[0,-1,0]]]],dtype=torch.float32)
    S_kernel = torch.tensor([[[[0,1/64,0],[1/64,12/64,1/64],[0,1/64,0]]]],dtype=torch.float32)
    Pi_kernel = torch.tensor([[[[0,0,0],[0,0,0],[0,0,0]]]],dtype=torch.float32)
    R_kernel = torch.tensor([[[[0,0.5,0.5],[0.5,1,0.5],[0.5,0.5,0]]]],dtype=torch.float32)
    RT_kernel = torch.tensor([[[[0,0.5,0.5],[0.5,1,0.5],[0.5,0.5,0]]]],dtype=torch.float32)

    if A is not None:
        A.weight = torch.nn.Parameter(A_kernel)

    if S is not None:
        S.weight = torch.nn.Parameter(S_kernel)

    if Pi is not None:
        Pi.weight = torch.nn.Parameter(Pi_kernel)

    if R is not None:
        R.weight = torch.nn.Parameter(R_kernel)

    if RT is not None:
        RT.weight = torch.nn.Parameter(RT_kernel)

    return

##### For MG: setup for prolongation and error calculation #####
RT = nn.ConvTranspose2d(1, 1, kernel_size=3, stride=2, padding=0, bias=False)
get_mg_init(None,None,None,None,RT)
A = nn.Conv2d(1, 1, kernel_size=3,stride=1, padding=1, bias=False)
get_mg_init(A,None,None,None,None)

class MgIte(nn.Module):
    def __init__(self, A, S):
        super().__init__()

        get_mg_init(A=A,S=S)                       ##### For MG: inilization of A, S #####

        self.A = A
        self.S = S

        self.bn1 =nn.BatchNorm2d(A.weight.size(0)) ##### For MgNet: BN #####
        self.bn2 =nn.BatchNorm2d(S.weight.size(0)) ##### For MgNet: BN #####

    def forward(self, out):
        u, f = out
        u = u + (self.S(((f-self.A(u)))))                                   ##### For MG: u = u + S*(f-A*u) #####
        u = u + F.relu(self.bn2(self.S(F.relu(self.bn1((f-self.A(u))))))) ##### For MgNet: add BN and ReLU #####
        out = (u, f)
        return out

class MgRestriction(nn.Module):
    def __init__(self, A_old, A, Pi, R):
        super().__init__()

        get_mg_init(A=A,Pi=Pi,R=R)            ##### For MG: inilization of A, Pi, R #####

        self.A_old = A_old
        self.A = A
        self.Pi = Pi
        self.R = R

        self.bn1 = nn.BatchNorm2d(Pi.weight.size(0))   ##### For MgNet: BN #####
        self.bn2 = nn.BatchNorm2d(R.weight.size(0))    ##### For MgNet: BN #####

    def forward(self, out):
        u_old, f_old = out
        u = self.Pi(u_old)                              ##### For MG: u = Pi*u_old #####
        f = self.R(f_old-self.A_old(u_old)) + self.A(u) ##### For MG: f = R*(f_old-A_old*u_old) + A*u #####
        u = F.relu(self.bn1(self.Pi(u_old)))                              ##### For MgNet: add BN and ReLU #####
        f = F.relu(self.bn2(self.R(f_old-self.A_old(u_old)))) + self.A(u) ##### For MgNet: add BN and ReLU #####
        out = (u,f)
        return out

class MG(nn.Module):
    def __init__(self, num_channel_input, num_iteration, num_channel_u, num_channel_f, num_classes):
        super().__init__()
        self.num_iteration = num_iteration
        self.num_channel_u = num_channel_u

        ##### For MgNet: Initialization layer #####
        self.conv1 = nn.Conv2d(num_channel_input, num_channel_f, kernel_size=3, stride=1, padding=1, bias=False)
        self.bn1 = nn.BatchNorm2d(num_channel_f)

        A = nn.Conv2d(num_channel_u, num_channel_f, kernel_size=3, stride=1, padding=1, bias=False)
        S = nn.Conv2d(num_channel_f, num_channel_u, kernel_size=3,stride=1, padding=1, bias=False)

        layers = []
        for l, num_iteration_l in enumerate(num_iteration): #l: l-th layer.   num_iteration_l: the number of iterations of l-th layer
            for i in range(num_iteration_l):
                layers.append(MgIte(A, S))

            setattr(self, 'layer'+str(l), nn.Sequential(*layers))
            # set attribute. This is equivalent to define
            # self.layer1 = nn.Sequential(*layers)
            # self.layer2 = nn.Sequential(*layers)
            # ...
            # self.layerJ = nn.Sequential(*layers)

            if l < len(num_iteration)-1:
                A_old = A

                A = nn.Conv2d(num_channel_u, num_channel_f, kernel_size=3,stride=1, padding=1, bias=False)
                S = nn.Conv2d(num_channel_f, num_channel_u, kernel_size=3,stride=1, padding=1, bias=False)

                ##### For MG: padding=0 #####
                Pi = nn.Conv2d(num_channel_u, num_channel_u, kernel_size=3,stride=2, padding=0, bias=False)
                R  = nn.Conv2d(num_channel_f, num_channel_f, kernel_size=3, stride=2, padding=0, bias=False)

                ##### For MgNet: padding=1 #####
                Pi = nn.Conv2d(num_channel_u, num_channel_u, kernel_size=3,stride=2, padding=1, bias=False)
                R = nn.Conv2d(num_channel_f, num_channel_f, kernel_size=3, stride=2, padding=1, bias=False)

                layers= [MgRestriction(A_old, A, Pi, R)]

        ##### For MgNet: average pooling and fully connected layer for classification #####
        self.pooling = nn.AdaptiveAvgPool2d(1)  # pooling the data in each channel to size=1
        self.fc = nn.Linear(num_channel_u ,num_classes)

    def forward(self, u, f):
        f = F.relu(self.bn1(self.conv1(f)))                 ##### For MgNet: initialization of f #####
        if use_cuda:                                        ##### For MgNet: initialization of u #####
            u = torch.zeros(f.size(0),self.num_channel_u,f.size(2),f.size(3), device=torch.device('cuda'))
        else:
            u = torch.zeros(f.size(0),self.num_channel_u,f.size(2),f.size(3))

        out = (u, f)

        u_list.append(u)                                      ##### For MG: save u^j, j=1,2,...,J #####
        for l in range(len(self.num_iteration)):
            out = getattr(self, 'layer'+str(l))(out)
            u, f = out                                        ##### For MG: save u^j, j=1,2,...,J #####
            u_list.append(u)                                  ##### For MG: save u^j, j=1,2,...,J #####

        ##### For MgNet: average pooling and fully connected layer for classification #####
        u, f = out
        u = self.pooling(u)
        u = u.view(u.shape[0], -1)
        u = self.fc(u)
        return u

```

    Use GPU? False

## 2. Apply Multigrid to solve the solving the following system

\begin{equation}\label{matrix}
A\ast u =f,
\end{equation}
where $A\ast$ is a convolution for one channel with stride 1 and zero padding $1$
$$
A=\begin{bmatrix} 0 & -1 & 0 \\ -1 & 4 & -1 \\ 0 & 1 & 0 \end{bmatrix},~~
$$
and
$ u \in \mathbb{R}^{n\times n} $, $ f\in \mathbb{R}^{n\times n}$ and $f_{i,j}=\dfrac{1}{(n+1)^2}$

### Mulrigrid code includes: (a) comment out the code in 1 for MgNet; (b) the setup and postprocessing code below

```python
def plot_solution(J,u,label_name):
    N = 2 ** J -1
    h = 1/2**J
    X = np.arange(h, 1, h)
    Y = np.arange(h, 1, h)
    X, Y = np.meshgrid(X,Y)   # create a mesh
    a = torch.reshape(u, (N, N))
    fig1 = plt.figure()
    ax = Axes3D(fig1)         # plot a 3D surface, (X,Y,u(X,Y))
    ax.plot_surface(X, Y, np.array(a.data), rstride=1, cstride=1, cmap=plt.cm.coolwarm)
    ax.set_title(label_name)

def plot_error(M,error,label_name):
    #print(np.linalg.norm((f-self.A(u)).reshape(-1).detach().numpy()))
    plt.figure()
    plt.title('Error vs number of iterations using '+label_name)
    plot = plt.plot(error)
    plt.xlabel('Number of iterations')
    plt.yscale('log')
    plt.ylabel('Error')
    plt.show()
```

```python
def MG1(u,f,J,num_iteration):
    u_list.clear()                                     # Save u^0,u^1,u^2,u^3...,u^J
    u = MG0(u,f)
    for j in range(J-1,0,-1):
        u_list[j] += RT(u_list[j+1])
    u = u_list[1]
    return u
```

```python
# Model setup
num_channel_input = 1
num_channel_u = 1
num_channel_f = 1
num_classes = 1

J = 4
num_iteration = [2,2,2,2]
MG0=MG(num_channel_input, num_iteration, num_channel_u, num_channel_f, num_classes)

##### For MG: PDE setup u=sin(2*pi*x)*sin(2*pi*y) #####

N = 2 ** J -1
h = 1/2**J
u_exact = torch.ones(1,1,N,N)
f = torch.ones(1,1,N,N) / (N+1) **2

##### For MG: Muligrid iteration #####
M = 100
u = torch.randn(1,1,N,N)
error = [np.linalg.norm((A(u)-f).detach().numpy())]  # calculate the Frobenius Norm of (A*u-f)
u_list =[]                                     # Save u^0,u^1,u^2,u^3...,u^J

for m in range(M):
    u = MG1(u,f,J,num_iteration)
    error.append(np.linalg.norm((A(u)-f).detach().numpy())) # calculate the Frobenius Norm of (A*u-f)

##### Lian added for MG: Plot results #####
plot_error(M,error,'Multigrid')
plot_solution(J,u,'Numerical solution')
```

## 3. Build and training MgNet on Cifar10

### MgNet code includes: (a) comment out the code in 1 for Multigrid; (b) the setup, training and test code below

```python
def adjust_learning_rate(optimizer, epoch, init_lr):
    #lr = 1.0 / (epoch + 1)
    lr = init_lr * 0.1 ** (epoch // 30)
    for param_group in optimizer.param_groups:
        param_group['lr'] = lr
    return lr

minibatch_size = 128
num_epochs = 120
lr = 0.1

num_channel_input = 3
num_channel_u = 64
num_channel_f = 64
num_classes = 10
num_iteration = [1,1,1,1]

# Step 1: Define a model
my_model = MgNet(num_channel_input, num_iteration, num_channel_u, num_channel_f, num_classes)

if use_cuda:
    my_model = my_model.cuda()

# Step 2: Define a loss function and training algorithm
criterion = nn.CrossEntropyLoss()
optimizer = optim.SGD(my_model.parameters(), lr=lr, momentum=0.9, weight_decay = 0.0005)

# Step 3: load dataset
normalize = torchvision.transforms.Normalize(mean=(0.4914, 0.4822, 0.4465), std=(0.2023, 0.1994, 0.2010))

transform_train = torchvision.transforms.Compose([torchvision.transforms.RandomCrop(32, padding=4),
                                                  torchvision.transforms.RandomHorizontalFlip(),
                                                  torchvision.transforms.ToTensor(),
                                                  normalize])

transform_test  = torchvision.transforms.Compose([torchvision.transforms.ToTensor(),normalize])

trainset = torchvision.datasets.CIFAR10(root='./data', train=True, download=True, transform=transform_train)
trainloader = torch.utils.data.DataLoader(trainset, batch_size=minibatch_size, shuffle=True)

testset = torchvision.datasets.CIFAR10(root='./data', train=False, download=True, transform=transform_test)
testloader = torch.utils.data.DataLoader(testset, batch_size=minibatch_size, shuffle=False)

# classes = ('plane', 'car', 'bird', 'cat', 'deer', 'dog', 'frog', 'horse', 'ship', 'truck')

start = timer()

#Step 4: Train the NNs
# One epoch is when an entire dataset is passed through the neural network only once.
for epoch in range(num_epochs):
    start_epoch = timer()

    current_lr = adjust_learning_rate(optimizer, epoch, lr)

    start_training = timer()
    my_model.train()
    for i, (images, labels) in enumerate(trainloader):
        if use_cuda:
          images = images.cuda()
          labels = labels.cuda()

        # Forward pass to get the loss
        outputs = my_model(0,images)   # We need additional 0 input for u in MgNet
        loss = criterion(outputs, labels)

        # Backward and compute the gradient
        optimizer.zero_grad()
        loss.backward()  #backpropragation
        optimizer.step() #update the weights/parameters
    end_training = timer()
    print('Computation Time for training:',end_training - start_training)

  # Training accuracy
    start_training_acc = timer()
    my_model.eval()
    correct = 0
    total = 0
    for i, (images, labels) in enumerate(trainloader):
        with torch.no_grad():
          if use_cuda:
              images = images.cuda()
              labels = labels.cuda()
          outputs = my_model(0,images)  # We need additional 0 input for u in MgNet
          p_max, predicted = torch.max(outputs, 1)
          total += labels.size(0)
          correct += (predicted == labels).sum()
    training_accuracy = float(correct)/total
    end_training_acc = timer()
    print('Computation Time for training accuracy:',end_training_acc - start_training_acc)

    # Test accuracy
    start_test_acc = timer()
    correct = 0
    total = 0
    for i, (images, labels) in enumerate(testloader):
        with torch.no_grad():
          if use_cuda:
              images = images.cuda()
              labels = labels.cuda()
          outputs = my_model(0,images)      # We need additional 0 input for u in MgNet
          p_max, predicted = torch.max(outputs, 1)
          total += labels.size(0)
          correct += (predicted == labels).sum()
    test_accuracy = float(correct)/total
    end_test_acc = timer()
    print('Computation Time for test accuracy:',end_test_acc - start_test_acc)

    print('Epoch: {}, learning rate: {}, the training accuracy: {}, the test accuracy: {}' .format(epoch+1,current_lr,training_accuracy,test_accuracy))
    end_epoch = timer()
    print('Computation Time for one epoch:',end_epoch - start_epoch)

end = timer()
print('Total Computation Time:',end - start)
```

```python

```


text:
## Linear Population Differential Equations

In this notebook we will a simple population equation [1]

```python
## Library
import numpy as np
import math
import matplotlib.pyplot as plt
```

## Linear Population Equation

\begin{equation}
\frac{d p}{dt}=\alpha p
\end{equation}
for $0\leq t \leq 10$ years, where $\alpha$ is the birth rate $\alpha_{Birth}$ minus death rate $\alpha_{Death}$, with the initial condition
$$ p(0)=A. $$
Which has the exact solution
\begin{equation}
p(t)=Aexp(\alpha t)
\end{equation}
In this simulation we set the parameters for the different conditions:
* visual condition: $\alpha=0.3$.

```python
alpha=0.3
```

To numiercally approximate the solution we discretise the continuous domain using a step size, $h=1$,
which gives
\begin{equation}
t_i=0+i h,
\end{equation}
for $i=0,...10$.

The Figure below illustrate the discrete domain.

```python
h=1
t_end=10
t=np.arange(0,t_end+h/2,h)
fig = plt.figure(1,figsize=(5,4))
ax = fig.add_subplot(1,1,1)
ax.plot(t,0*t,'o:',color='k')
ax.set_xlabel('Time (yrs)')
plt.show()
```

```python
N=int(t_end/h) # Number of Steps
P=np.zeros(N+1) # Numerical Solution P
P[0]=3
p=3*np.exp(alpha*t)
```

We numerically approximate the differential equation using Eulers method to give the difference equation
\begin{equation}
P_{i+1}=P_{i}+h \big(\alpha P_{i}k \big),
\end{equation}
where $h$ is the stepsize and where $P_i$ is the numerical approximation of $p(t_i)$ at distance $t_i$, for $i=0,...,10$

```python
for i in range (1,N+1):
    P[i]=P[i-1]+h*(alpha*P[i-1])

```

## Results
The plot bellow shows the Numerical solution $P_i$ (left), the analytic solution $p(x)$ (middle) and the error $P(x_i)-P_i$ (right) as a function of time $t$.

```python
## Plotting Figure
fig = plt.figure(1,figsize=(12,4))

# --- left hand plot
ax = fig.add_subplot(1,3,1)
ax.plot(t,P,'o:',color='red',label='Population')
ax.legend()
ax.plot(t,t,'--k')

ax.set_xlabel('Time (yrs)')
ax.set_ylabel('Population (billions)')
#ax.legend(loc='best')
ax.set_title('Numerical Solution, P')

# --- middle plot
ax = fig.add_subplot(1,3,2)
ax.plot(t,p,'o:',color='red',label='Population')

ax.plot(t,t,'--k')
ax.legend()

ax.set_xlabel('Time (yrs)')
ax.set_ylabel('Population (billions)')

ax.set_title('Analytic (exact) Solution p(t)')

# --- right hand plot
ax = fig.add_subplot(1,3,3)
ax.plot(t,P-p,'o:',color='red',label='VIS ERROR')
ax.legend()

ax.set_title('Error')
ax.set_ylabel('Error(Population)')
ax.set_xlabel('Time (yrs)')

# --- title, explanatory text and save
fig.suptitle('Linear Population Equation', fontsize=20)
plt.tight_layout()
plt.subplots_adjust(top=0.85)

```

# References

[1] Stover, Christopher and Weisstein, Eric W. "Population Growth." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/PopulationGrowth.html

max_stars_repo_path: .ipynb_checkpoints/04_Logisitic Population Differential Equations-checkpoint.ipynb


text:
| |Pierre Proulx, ing, professeur|
|:---|:---|
|Département de génie chimique et de génie biotechnologique |** GCH200-Phénomènes d'échanges I **|

### Section 2.3, écoulement d'un film de fluide Newtonien dans une conduite

>  Le développement fait dans Transport Phenomena sera répété ici en développant les solutions avec le calculateur formel sympy et en traçant la solution avec sympy.plot. Nous n'irons pas ici autant dans le détail que dans la section 2.2, Les coordonnées cylindriques sont utilisées, évidemment. Remarquez bien comment la courbure de la géométrie a un effet sur le bilan.

>Vous pouvez suivre le développement détaillé à l'aide du livre et voir comment les équations de bilans initiales deviennent les profils de vitesse, de force, débit, etc...

```python
#
# Pierre Proulx
#
#
# Préparation de l'affichage et des outils de calcul symbolique
#
import sympy as sp
from IPython.display import *
sp.init_printing(use_latex=True)
%matplotlib inline
```

```python
# Paramètres, variables et fonctions
#
r,delta_r,L,rho,g,mu,R,P_0,P_L=sp.symbols('r,delta_r,L,rho,g,mu,R,P_0,P_L')
C1=sp.symbols('C1')
phi_rz,phi_zz=sp.symbols('phi_rz,phi_zz')
```

```python
#
# Équation du bilan des forces
#
eq=2*sp.pi*r*L*phi_rz(r)      -2*sp.pi*(r+delta_r)*L*phi_rz(r+delta_r)
eq=eq+2*sp.pi*delta_r*r*phi_zz(0)-2*sp.pi*r*delta_r*phi_zz(L)
eq=eq+2*sp.pi*r*delta_r*L*rho*g
#eq=r*phi_rz(r)/delta_r-(r+delta_r)*phi_rz(r+delta_r)/delta_r+r*P_0/L-r*P_L/L
#
# pose de la limite quand dr tends vers 0
#
eq=eq/(2*sp.pi*L*delta_r)
display(eq)
```

```python
eq1=sp.limit(eq,delta_r,0).doit()
display(eq1)
```

```python
eq1=eq1.subs(phi_zz(0),P_0)     # les termes de convection s'annulent comme en 2.2
eq1=eq1.subs(phi_zz(L),P_L)     # mais il reste le terme de pression
display(eq1)
```

```python
# Loi de Newton, le terme phi ne contient que le cisaillement, pas de vitesse radiale.
tau=sp.symbols('tau_rz')
eq2=sp.Eq(eq1.subs(phi_rz(r),tau(r)))
display(eq2)
```

```python
eq2=sp.dsolve(eq2)
display(eq2)
```

```python
vz=sp.Function('v_z')(r)
newton=-mu*sp.Derivative(vz,r)
eq3=eq2.subs(tau(r),newton)
eq3=eq3.subs(C1,0)               # C1 doit être nulle sinon on a une force infinie en r=0
display(eq3)
```

```python
eq4=sp.dsolve(eq3,vz)      # rhs pour avoir la droite de l'équation vz=...
display(eq4)
```

```python
constante=sp.solve(sp.Eq(eq4.rhs.subs(r,R),0),'C1',dict=True)   # Condition de non-glissement à la paroi
display(constante)
constante=constante[0]
```

```python
vz=eq4.subs(constante)
print('vitesse en fonction de r')
display(vz.simplify())
```

```python
print('Force')
#
# collect et simplify simplifient l'expression, c'est purement cosmétique.
#
display(-mu*vz.lhs.diff(r)*2*sp.pi*r*L)
display('évalué en r=R devient')
display(-mu*vz.rhs.diff(r).subs(r,R)*2*sp.pi*R*L)
```

```python
print( 'débit et vitesse moyenne')
vz=vz.rhs # On enlève le terme de gauche et à partir d'ici on travaille seulement avec la partie droite
debit=sp.integrate(2*sp.pi*vz*r,(r,0,R))
display(debit)
```

```python
vmax=vz.subs(r,0)
display(debit/(sp.pi*R**2)/vmax)   # ce calcul donne quoi?
```

```python
# Maintenant traçons le profil en donnant des valeurs de paramètres réalistes
#
dico={'rho':1000, 'mu':0.001, 'R':0.05, 'L':100,'P_0':200,'P_L':0,'g':9.81}
vz=vz.subs(dico)
vzmax=vz.subs(r,0)
display(vz,vzmax)
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize']=10,8
sp.plot(vz.subs(dico)/vzmax,(r,0,R.subs(dico)),title='Profil de vitesse parabolique',
       ylabel='V/Vmax',xlabel='r/R');
```


text:
# Exercise 12) Policy Gradients

In this exercise we will have a look at policy gradients.
The theory of policy gradients applies to function approximators that decide on which action to choose.
The function approximators we met in the past were employed to estimate the (action) value function.
Since their task was to judge the quality of the current situation they are often referred to as "critics".
In contrary, we can also use a function approximator to directly choose an action; we call these "actors".
Why should we do that if we made it work with nothing more than a critic?
Because this will finally allow us to make use of contiuous action spaces! Eureka!

In this exercise we will use a new `gym` environment `LunarLanderContinuous-v2`.
To run this environment please make sure to have `Box2D` installed: `pip install Box2D`.

Source: https://www.billyvreeland.com/portfolio/2017/1013/solving-openai-gym-nm4yz

The main task is to land the LunarLander in the landing zone.
An accident-free landing is defined by both legs coming into  ground contact with moderate velocity.
We are dealing with a continuous state and action space as defined below.
Please notice that the control functions for main and side engines contain a dead zone in which the engines are inactive.
The reward is mainly defined depending on whether the landing procedure is successful (+100) or not (-100).
Firing the main engine gives a small negative reward.
The problem is solved if a return of at least 200 is reached.
For more information see https://gym.openai.com/envs/LunarLanderContinuous-v2/.

\begin{align}
\text{state}&=
\begin{bmatrix}
p_x\\
p_y\\
v_x\\
v_y\\
\varphi\\
20 \, \omega\\
1 \text{ if left leg has ground contact, else } 0\\
1 \text{ if right leg has ground contact, else } 0\\
\end{bmatrix}
\\
\text{action}&=
\begin{bmatrix}
\text{main engine: } [-1, 0] \rightarrow \text{off}, ]0, 1] \rightarrow [50 \, \%, 100 \, \%] \text{ of available power}\\
\text{side engines: } [-1, -0.5] \rightarrow [50 \, \%, 100 \, \%] \text{ of available right engine power}, [0.5, 1] \rightarrow [50 \, \%, 100 \, \%] \text{ of available left engine power}\\
\end{bmatrix}
\end{align}

```python
import numpy as np
import pandas as pd
import gym
import random
import time
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
plt.style.use('seaborn-talk')
from mpl_toolkits.mplot3d import Axes3D
from collections import deque
from sklearn.kernel_approximation import RBFSampler
import sklearn.pipeline
import sklearn.preprocessing

import tensorflow as tf
from tensorflow.keras.models import Sequential, Model
from tensorflow.keras.layers import Dense, Lambda, Input
from tensorflow.keras.optimizers import SGD
```

```python
env = gym.make('LunarLanderContinuous-v2')

state = env.reset()
while True:
    env.render()
    state, reward, done, _ = env.step(env.action_space.sample())

    if done:
        break

env.close()
```

## 1) Monte Carlo Policy Gradient
Write a REINFORCE algorithm.

Execute the follwoing cell to fit the featurizer using RBFSampler, like already learned in the last exercises.

```python
state_array = []
state = env.reset()

for i in tqdm(range(1000)):
    state = env.reset()
    while True:
        state, reward, done, _ = env.step(env.action_space.sample())
        state_array.append(state)

        if done:
            break

state_array = np.array(state_array)

featurizer = sklearn.pipeline.make_pipeline(
    sklearn.preprocessing.StandardScaler(),
    sklearn.pipeline.FeatureUnion([
    ("rbf0", RBFSampler(gamma=5.0, n_components = 1000)),
    ("rbf1", RBFSampler(gamma=2.0, n_components = 1000)),
    ("rbf2", RBFSampler(gamma=1.0, n_components = 1000)),
    ("rbf3", RBFSampler(gamma=0.5, n_components = 1000)),
    ]),
    sklearn.preprocessing.StandardScaler()
)

_ = featurizer.fit(state_array)

```

Use the following cell to define the function approximators for the policy.
As seen in Algo.12.1 we need to calculate $\nabla_\theta \mathrm{ln}\,\pi(u_k | x_k, \theta)$.
$\pi$ is herein defined as the normal distribution :
\begin{align}
\pi(u_k | x_k, \theta) & = \frac{\mathrm{exp} \left( {-\frac{1}{2} (u_k - \mu_\theta(x_k))^\mathrm{T} \mathbf{\Sigma}^{-1}_\theta(x_k) (u_k - \mu_\theta(x_k))} \right)}{\sqrt{(2\pi)^p \mathrm{det}(\mathbf{\Sigma}_\theta(x_k))}},\\
\text{with}\hspace{1em} p & = \mathrm{dim}(u_k).
\end{align}

Extend `loglikelyhoodGaussian` such that it returns $\mathrm{ln}\,\pi(u_k | x_k, \theta)$!
Use the numpy equivalent `TensorFlow` functions (e.g. `tf.linalg.inv()`).
`TensorFlow` functions are differentiable and can therefore be  used to calculate $\nabla_\theta \mathrm{ln}\,\pi(u_k | x_k, \theta)$.

```python
state = np.reshape(env.reset(), (1, -1))
feature_state = featurizer.transform(state)
input_dim = feature_state.shape[1]
action_space_dim = len(env.action_space.sample())

# define policy
def create_policy():
    input_layer = Input(shape=input_dim)

    hidden1 = Dense(400, activation='linear')(input_layer)
    hidden1 = tf.keras.layers.LeakyReLU(alpha=0.1)(hidden1)

    hidden2_mu = Dense(400, activation='linear')(hidden1)
    hidden2_mu = tf.keras.layers.LeakyReLU(alpha=0.1)(hidden2_mu)

    hidden2_sigma = Dense(400, activation='linear')(hidden1)
    hidden2_sigma = tf.keras.layers.LeakyReLU(alpha=0.1)(hidden2_sigma)

    mu_out = Dense(action_space_dim, activation='linear')(hidden2_mu)
    mu_out = Lambda(lambda x: tf.reshape(tf.clip_by_value(x, -1, 1), (-1, 1)))(mu_out)
    sigma_out = Dense(action_space_dim, activation='softplus')(hidden2_sigma)
    sigma_out = Lambda(lambda x: tf.linalg.diag(tf.clip_by_value(tf.reshape((x @ np.array([[0.01, 0], [0, 0.1]])), [-1]), 1e-4, 1)))(sigma_out)

    policy = Model(inputs=input_layer, outputs=[mu_out, sigma_out])

    return policy

policy = create_policy()
theta = policy.get_weights()

# Regularization; downscaling of the network parameters to prevent divergence
for i in range(len(theta)):
    theta[i] = theta[i] * 0.4

policy.set_weights(theta)

@tf.function
def loglikelyhoodGaussian(x, _mu, _sigma):
    x = tf.cast(x, tf.float32)

    # YOUR CODE HERE
    raise NotImplementedError()

```

Use the following template to write a REINFORCE algorithm.
This time the Adam (adaptive moment estimation) optimizer is used, which is an enhanced SGD optimizer.
For more informations see https://arxiv.org/abs/1412.6980.

```python
env.close()
alpha_policy = 1e-5
gamma = 0.99
nb_episodes = 5000

env = gym.make('LunarLanderContinuous-v2')
return_history = []

optimizer = tf.optimizers.Adam(learning_rate=alpha_policy)

for j in tqdm(range(nb_episodes)):
    k = 0
    accumulated_rewards = 0

    states = []
    actions = []
    rewards = []
    probs_log = []

    state = env.reset()
    state = state.reshape(1, -1)

    with tf.GradientTape() as tape:
        while True:
            # env.render()
            k += 1

            # YOUR CODE HERE
            raise NotImplementedError()

```

Plot the learning curve of the training process!

```python
# YOUR CODE HERE
raise NotImplementedError()
```

## Execution

Use `deterministic` to choose between deterministic execution (applying $\mu$ directly) or taking the stochastic action by sampling from the normal distribution.

```python
env.close()
env = gym.make('LunarLanderContinuous-v2')

deterministic = True

for j in tqdm(range(10)):

    state = env.reset()
    accumulated_rewards = 0

    while True:
        env.render()

        ### STEP)
        feat_state = featurizer.transform(np.reshape(state, (1, -1)))
        mu, sigma = policy(feat_state)
        print(mu)
        if deterministic:
            action = np.reshape(mu.numpy(),(-1,))
        else:
            action = np.random.multivariate_normal(mean=np.squeeze(mu.numpy()), cov=sigma.numpy())

        next_state, reward, done, _ = env.step(action)

        accumulated_rewards += reward

        state = next_state

        if done:
            if j % 1 == 0:
                print(f"Episode {j}, Length: {k}")
                print(f"Return {accumulated_rewards}")
                print()
            break

    env.close()

```

## 2) Actor-Critic with TD(0) Targets

Write an actor-critic (AC) algorithm to land the lander in the landing zone :)

Use the following cell to create two function approximators. One to estimate the state values (critic) and one to decide on the actions to take (actor).

```python
state = np.reshape(env.reset(), (1, -1))
input_dim = len(featurizer.transform(state)[0])
action_space_dim = len(env.action_space.sample())

# define critic
def create_critic():
    critic = Sequential()
    critic.add(Dense(400, activation='linear', input_dim=input_dim))
    critic.add(tf.keras.layers.LeakyReLU(alpha=0.1))
    critic.add(Dense(400, activation='linear'))
    critic.add(tf.keras.layers.LeakyReLU(alpha=0.1))
    critic.add(Dense(1, activation='linear'))

    return critic

critic = create_critic()
w = critic.get_weights()
for i in range(len(w)):
    w[i] = w[i] * 0.2
critic.set_weights(w)

# define actor
def create_actor():
    input_layer = Input(shape=input_dim)

    hidden1 = Dense(400, activation='linear')(input_layer)
    hidden1 = tf.keras.layers.LeakyReLU(alpha=0.1)(hidden1)

    hidden2 = Dense(400, activation='linear')(hidden1)
    hidden2 = tf.keras.layers.LeakyReLU(alpha=0.1)(hidden2)

    mu = Dense(action_space_dim, activation='linear')(hidden2)
    mu = Lambda(lambda x: tf.reshape(tf.clip_by_value(x, -1, 1), (-1, 1)))(mu)

    sigma = Dense(action_space_dim, activation='softplus')(hidden2)
    sigma = Lambda(lambda x: tf.linalg.diag(tf.clip_by_value(tf.reshape(x @ np.array([[0.01,   0], [   0, 0.1]]), [-1]), 1e-4, 1)))(sigma)

    actor = Model(inputs=input_layer, outputs=[mu, sigma])

    return actor

actor = create_actor()
theta = actor.get_weights()
for i in range(len(theta)):
    theta[i] = theta[i] * 0.4
actor.set_weights(theta)

```

Use the following template to write an AC algorithm.

```python
env.close()

alpha_critic = 1e-4
alpha_actor = 1e-5
gamma = 0.99
nb_episodes = 2000

env = gym.make('LunarLanderContinuous-v2')

return_history = []

for j in tqdm(range(nb_episodes)):
    k = 0
    accumulated_rewards = 0

    state = env.reset()
    I = 1

    while True:
        # env.render()
        k += 1

        with tf.GradientTape(persistent=True) as tape:
            # persistent means the tape is not deleted after the first gradient has been computed

            # YOUR CODE HERE
            raise NotImplementedError()

        if done:
            if j % 250 == 0:
                plt.plot(return_history, label='Return')
                plt.plot(pd.Series(return_history, name='reward_history').rolling(10).mean(), label='MA')
                plt.xlabel('episode')
                plt.ylabel('return')
                plt.grid(True)
                _=plt.legend()
                plt.show()
            return_history.append(accumulated_rewards)
            break

    env.close()

```

Plot the learning curve of the training process!

```python
# YOUR CODE HERE
raise NotImplementedError()
```

## Execution

Use `deterministic` to choose between deterministic execution (apply $\mu$) directly or take the stochastic action by sampling from the normal distribution.

```python
env.close()
env = gym.make('LunarLanderContinuous-v2')

deterministic = True

for j in tqdm(range(10)):

    state = env.reset()
    accumulated_rewards = 0

    while True:
        env.render()

        ### STEP)
        feat_state = featurizer.transform(np.reshape(state, (1, -1)))
        mu, sigma = actor(feat_state)

        if deterministic:
            action = np.reshape(mu.numpy(), (-1,))
        else:
            action = np.random.multivariate_normal(mean=np.squeeze(mu.numpy()), cov=sigma.numpy())

        next_state, reward, done, _ = env.step(action)

        accumulated_rewards += reward

        state = next_state

        if done:
            if j % 1 == 0:
                print(f"Episode {j}, Length: {k}")
                print(f"Return {accumulated_rewards}")
                print()
            break

    env.close()

```


text:
# Пишем свой фреймворк

```python
# только numpy, только хардкор
import numpy as np
```

Мотивация: конечным пользователям вашего фреймворка не хочется думать, как они работают слои внутри. Им просто хочется объявить последовательность элементарных операций над входными данными, а о градиентах и прочем матане пусть позаботится сам фреймворк.

**Module** — это абстрактный класс, от которого будут наследоваться слои нашей нейронной сети. Абстрактные классы нужны, чтобы можно было реализовывать не все методы, а только переопределить некоторые. Все в лучших традициях ООП.

Модуль — это такая чёрная коробка, которая
1. Умеет принимать какие-то входные данные $X$ и возращать какие-то выходные данные $Y$ (`forward`)
2. Возможно, имеет какие-то параметры, которые можно изменять, (`parameters`, `grad_parameters`)
3. Будучи встроенной в вычислительный граф, умеет по градиенту относительно своих выходных значений вычислять градиент относительно входных данных, а также собственных параметров (`backward`)
4. Умеет переключаться в режимы обучения и инференса, если они отличаются (`train`, `eval`)

Теперь поподробнее.

## Входные данные

Современные нейросети оптимизируют различными вариантами стохастического градиентного спуска, и мы тоже будем его использовать. Его отличие от обычного в том, что на каждом шаге мы не считаем градиент на всем датасете (это было бы слишком долго), а оцениваем его, усреднив градиенты на его случайно выбранной малой части, которую называют батчем (`batch`). Если батч формировать случайно, и если его размер достаточно большой, то мы можем быстро получить немного шумную, но приемлемую для нас оценку градиента, и не прогонять через сеть все миллионы примеров ради одного маленького шага. Эта интуиция ограничивает размер батча сверху.

Математик бы принял время прогона одного примера по всей сети за константу и пришел бы к выводу, что нужно считать по одному примеру и делать каждый раз один шаг, но маленький. Это верное заключение, но в реальности, если увеличить размер батча в $k$ раз, то он будет работать не в $k$ раз дольше, а намного меньше.

Самая долгая операция в большинстве нейросетей — это перемножение матриц. Начиная с каких-то размеров для их перемножения имеет использовать алгоритм Штрассена, который работает уже быстрее, чем линейно. Проведём небольшой вычислительный эксперимент:

```python
'''A = np.random.randn(256, 2000)
B = np.random.randn(2000, 800)

# помножить каждый вектор-строку на B и сконкатенировать
%time C = np.stack(np.dot(A[i].T, B) for i in range(256))

# это то же самое, что использовать одно-большое матричное умножение
%time C = np.dot(A, B)'''
```

    Wall time: 149 ms
    Wall time: 40 ms

    C:\Users\surta\anaconda3\lib\site-packages\IPython\core\magic.py:187: FutureWarning: arrays to stack must be passed as a "sequence" type such as list or tuple. Support for non-sequence iterables such as generators is deprecated as of NumPy 1.16 and will raise an error in the future.
      call = lambda f, *a, **k: f(*a, **k)

Такая чисто вычислительная причина ограничивает размер батча снизу. На практике, в большинстве случаев оптимальный размер батча — несколько сотен. В случае с CPU это несколько десятков, потому что выгода от распараллеливания вычислений не такая сильная.

Вообще, почти все наши слои будут работать с векторами независимо, но из-за вычилсительных причин мы будем объединять их в матрицы. Вообще, более сложные нейросети работают с тензорами. «Тензор» это, вообще говоря, сложный математический объект, но в DL этот термин используется просто в занчении «многомерный массив». Например, картинки — это четырехмерный тензор: `[batch, channel, x, y]`.

### Forward

Эта функция просто принимает тензор (`numpy.ndarray`) и возвращает какой-то другой, над которым применили соответствующие операции.

Важный нюанс: нам позже для реализации `backward` почти всегда будет нужно сохранять где-нибудь выход `forward` (это создает очень большую нагрузку на память при обучении; [в принципе это можно и не делать](https://arxiv.org/pdf/1604.06174.pdf), но так проще). Условимся сохранять его в `self.output`, сразу после того, как посчитали.

### Параметры

Параметр модели — это что-то, что можно поодгонять, чтобы функция потерь стала меньше. Он должен быть доступен оптимизатору, а оптимизатору не обязательно знать, как всё у слоя все внутри работает. Ему нужны просто градиенты — насколько ему нужно подвинуть параметры сети, чтобы стало лучше.

Общаться с ним мы будем посредством двух функций: `params` и `grad_params`. Обе возвращают списки из тензоров — значения параметров и их посчитанных градиентов (см. `backward`) соответственно. Питон делает shallow copy, поэтому у оптимизатора так есть доступ на их изменение.

### Backward

После того, как мы в вычислительном графе все последовательно посчитали и дошли до функции потерь, нам надо подогнать параметры так, чтобы на тех же данных при повторном прогоне она стала меньше — иными словами, нам надо сделать шаг против градиента функции потерь относительно параметров сети.

Посчитать эти градиенты — нетривиальная задача. Мы могли бы рассмотреть каждый параметр по отдельности и как-нибудь посчитать градиент для него. Но это очень долго — параметров в современных сетях бывает по несколько миллионов.

Вместо этого мы применим трюк, основанный на формуле для производной сложной функции:

$$ f(g(x))' = f'(g(x)) \cdot g'(x) $$

Представьте, что часть сети от параметра до выхода — это всего две последовательно выполненные функции: $g$ и $f$. Тогда, согласно формуле, нам для этого параметра достаточно посчитать и перемножить две величины — $g'(x)$ (производная текущего слоя) и f'(g(x)) (производная относительно выхода текущего слоя).

Какие-то другие параметры могли тоже зависеть от производной относительно выхода.  и мы получаем выигрыш за счет того, что считаем её только один раз и запоминаем. Можно сказать, что мы применяем таким образом динамическое программирование на вычислительном графе, чтобы посчитать градиенты относительно всех его параметров.

Обратный прогон (`backward`) определяется для каждого слоя и нужен как раз для подсчета градиентов, имея градиент относительно своих выходных значений (аналог $f'(g(x)))$.

Он должен делать две вещи:

1. Посчитать градиент относительно собственных параметров.
2. Посчитать и вернуть градиент относительно своих входных данных.

Для лучшего понимания рассмотрите пример с `Linear` и `ReLU`.

### train / eval

Некоторые слои ведут себя по-разному во время обучечния и предсказания (`inference`). Обычно, это связано с разного вида регуляризацией — например, так ведут себя `BatchNorm` и `Dropout`.

По сути, для таких слоев нужно просто написать два разных `forward`-а для обучения и инференса.

```python
class Module():
    def __init__(self):
        self._train = True

    def forward(self, input):
        raise NotImplementedError

    def backward(self,input, grad_output):
        raise NotImplementedError

    def parameters(self):
        'Возвращает список собственных параметров.'
        return []

    def grad_parameters(self):
        'Возвращает список тензоров-градиентов для своих параметров.'
        return []

    def train(self):
        self._train = True

    def eval(self):
        self._train = False
```

Это **абстрактный класс** — от него наследуются другие слои, в которых эти функции будут реализованы.

# Sequential

**Sequential** будет оборачивать список модулей и выполнять их последовательно.

Это своего рода контейнер, внутри которого есть какой-то пайплайн.

Можно даже засовывать один Sequential внутри другого.

Многие не знают, но в питоне почти всегда для итерирования используется не **deep copy**, а **shallow copy**. Это делается для экономии памяти.

```python
class Sequential(Module):
    def __init__ (self, *layers):
        super().__init__()
        self.layers = layers

    def forward(self, input):
        """
        Прогоните данные последовательно по всем слоям:

            y[0] = layers[0].forward(input)
            y[1] = layers[1].forward(y_0)
            ...
            output = module[n-1].forward(y[n-2])

        Это должен быть просто небольшой цикл: for layer in layers...

        Хранить выводы ещё раз не надо: они сохраняются внутри слоев после forward.
        """

        for layer in self.layers:
            input = layer.forward(input)

        self.output = input
        return self.output

    def backward(self, input, grad_output):
        """
        Backward -- это как forward, только наоборот. (с)

        Предназначение backward:
        1. посчитать посчитать градиенты для собственных параметров
        2. передать градиент относительно своего входа

        О своих параметрах модули сами позаботятся. Нам же нужно позаботиться о передачи градиента.

            g[n-1] = layers[n-1].backward(y[n-2], grad_output)
            g[n-2] = layers[n-2].backward(y[n-3], g[n-1])
            ...
            g[1] = layers[1].backward(y[0], g[2])
            grad_input = layers[0].backward(input, g[1])

        Тут цикл будет уже чуть посложнее.
        """

        for i in range(len(self.layers)-1, 0, -1):
            grad_output = self.layers[i].backward(self.layers[i-1].output, grad_output)

        grad_input = self.layers[0].backward(input, grad_output)

        return grad_input

    def parameters(self):
        'Можно просто сконкатенировать все параметры в один список.'
        res = []
        for l in self.layers:
            res += l.parameters()
        return res

    def grad_parameters(self):
        'Можно просто сконкатенировать все градиенты в один список.'
        res = []
        for l in self.layers:
            res += l.grad_parameters()
        return res

    def train(self):
        for layer in self.layers:
            layer.train()

    def eval(self):
        for layer in self.layers:
            layer.eval()
```

# Слои

Приступим к реализации содержательной части — самих слоев.

На вход всех слоев будет подаваться матрица размера `batch_size` $\times$ `n_features` (см. описание `forward`).

Начнем с основного: линейный слой, он же fully-conected.

$$ Y = X W + b $$

Правильнее его называть афинным: после матричного умножения добавляется вектор $b$.

`forward` у него трививальный, а `backward` уже сложнее: нужно посчитать градиенты относительно трёх вещей:
1. Входных данных. Автор добродушен и спалит вам ответ, а вам нужно его доказать: $\nabla X = W^T (\nabla Y)$.
2. Матрица весов $W$. Тут нужно подумать, как каждый вес влияет на каждое выходное значение, и выразить ваши мысли линейной алгеброй.
3. Вектор $b$. С ним всё будет просто.

Не забудьте, что `grad_params` должен иметь такие же размерности, как и соответствующие параметры.

```python
class Linear(Module):
    def __init__(self, dim_in, dim_out):
        super().__init__()

        # Xavier initialization: инциализируем так,
        # что если на вход идет N(0, 1)
        # то и на выходе должно идти N(0, 1)
        stdv = 1./np.sqrt(dim_in)
        self.W = np.random.uniform(-stdv, stdv, size=(dim_in, dim_out))
        self.b = np.random.uniform(-stdv, stdv, size=dim_out)

    def forward(self, input):
        self.output = np.dot(input, self.W) + self.b
        return self.output

    def backward(self, input, grad_output):
        self.grad_b = np.mean(grad_output, axis=0)

        #     in_dim x batch_size
        self.grad_W = np.dot(input.T, grad_output)
        #                 batch_size x out_dim

        grad_input = np.dot(grad_output, self.W.T)

        return grad_input

    def parameters(self):
        return [self.W, self.b]

    def grad_parameters(self):
        return [self.grad_W, self.grad_b]
```

## Функции активации

**ReLU** — одна из самых простых функций активации:

$$
ReLU(x)=
\begin{cases}
x, & x > 0\\
0, & x \leq 0\\
\end{cases}
$$

`ReLU` это очень простой слой, поэтому автору не жалко её реализовать его за вас:

```python
class ReLU(Module):
    def __init__(self):
         super().__init__()

    def forward(self, input):
        self.output = np.maximum(input, 0)
        return self.output

    def backward(self, input, grad_output):
        grad_input = np.multiply(grad_output, input > 0)
        return grad_input
```

У ReLU есть проблема — у него бесполезная нулевая производная при $x < 0$.

[**Leaky Rectified Linear Unit**](http://en.wikipedia.org/wiki%2FRectifier_%28neural_networks%29%23Leaky_ReLUs) — это его модифицированная версия, имеющая в отрицательных координатах не нулевой градиент, а просто помноженный на маленькую константу `slope`.

$$
LeakyReLU_k(x)=
\begin{cases}
x, & x > 0\\
kx, & x \leq 0\\
\end{cases}
$$

При `slope` = 0 он превращается в обычный `ReLU`.

```python
class LeakyReLU(Module):
    def __init__(self, slope=0.03):
        super().__init__()

        self.slope = slope

    def forward(self, input):
        self.output = np.multiply(input, input > 0) + np.multiply(self.slope*input, input < 0)
        return self.output

    def backward(self, input, grad_output):
        grad_input = np.multiply(grad_output, input > 0) + np.multiply(self.slope*grad_output, input < 0)
        return grad_input
```

**Сигмоида** определяется формулой $\sigma(x) = \frac{1}{1+e^{-x}}$.

Когда-то она была самой часто используемой функции активации, потому что имела логичную вероятностную интерпретацию (вероятность наличия какой-то фичи), но потом перестали, потому что на очень больших или маленьких значениях её производные почти нулевые (см. проблема затухающего градиента).

Также используют [гипреболический тангенс](https://ru.wikipedia.org/wiki/%D0%93%D0%B8%D0%BF%D0%B5%D1%80%D0%B1%D0%BE%D0%BB%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%B8%D0%B5_%D1%84%D1%83%D0%BD%D0%BA%D1%86%D0%B8%D0%B8), который на самом деле просто сигмоида, отнормированная так, чтобы значения были в $[-1, 1]$: $tanh(x) = 2 \sigma(x) - 1$. Мы его отдельно реализовывать не будем.

Давайте посчитаем её производную:

$$
\begin{align}
\sigma'(x) &= (\frac{1}{1+e^{-x}})'
\\         &= \frac{e^{-x}}{(1+e^{-x})^2}
\\         &= \frac{1+e^{-x}-1}{(1+e^{-x})^2}
\\         &= \frac{1+e^{-x}}{(1+e^{-x})^2} - \frac{1}{(1+e^{-x})^2}
\\         &= \frac{1}{1+e^{-x}} - \frac{1}{(1+e^{-x})^2}
\\         &= \sigma(x) - \sigma(x)^2
\\         &= \sigma(x)(1 - \sigma(x))
\end{align}
$$

```python
class Sigmoid(Module):
    def __init__(self, slope=0.03):
        super().__init__()

    def forward(self, input):
        self.output = 1/(1+np.exp(-input))
        return self.output

    def backward(self, input, grad_output):
        grad_input = np.multiply(self.output*(1-self.output), grad_output)
        return grad_input
```

**Софтмакс** определяется так:

$$ \sigma(x)_k = \frac{e^{x_k}}{\sum_{i=1}^n e^{x_i} }$$

Можно заметить, что сигмоида — это частный случай софтмакса. Его можно интерпретировать как вероятностное распределение: его выходы положительны и суммируются в единицу. Поэтому его используют как последний слой для классификации.

Софтмакс — самый сложный с точки зрения написания `backward`. Как и все остальное, оно считается в 5 строчек кода, но [вывести их трудно](https://deepnotes.io/softmax-crossentropy).

```python
class SoftMax(Module):
    def __init__(self):
         super().__init__()

    def forward(self, input: np.ndarray) -> np.ndarray:
        # N.B. If inputs are large, then
        # exponents would be much larger
        self.output = np.subtract(input, input.max(axis=1, keepdims=True))
        self.output = np.exp(self.output)
        self.output = np.divide(self.output, self.output.sum(axis=1, keepdims=True))

        return self.output

    def backward(self, input, grad_output):
        jacobian = np.einsum('ij,jk->ijk', self.output, np.eye(self.output.shape[-1])) \
        - np.einsum('ij,ik->ijk', self.output, self.output)
        grad_input = np.einsum('ij, ijk->ik', grad_output, jacobian)
        return grad_input
```

## Регуляризация

Самый популярный регуляризатор в нейросетях — [**дропаут**](https://www.cs.toronto.edu/~hinton/absps/JMLRdropout.pdf). Идея простая: просто помножим поэлементно входные данные на случайную бинарную маску того же размера, как и сами данные. Сгенерировать маску можно через `np.random.binomial`.

Дропаута обычно хватает как единственного регуляризатора. Если вы заметите, что сеть оверфитится — просто добавьте его побольше.

**У дропаута разное поведение в режимах `train` и `eval`**. При `eval` он не должен делать ничего, а в `train` вместо применения маски нужно ещё домножить вход на $p$, чтобы скомпенсировать дропаут при обучении (так математическое ожидание значений будет такое же, как на трейне).

```python
class Dropout(Module):
    def __init__(self, p=0.5):
        super().__init__()

        self.p = p
        self.mask = None

    def forward(self, input):
        if self._train:
            self.mask = (np.random.random(input.size) < (1 - self.p)).reshape(input.shape)
            self.output = input * self.mask / (1 - self.p)
        else:
            self.output = input
        return self.output

    def backward(self, input, grad_output):
        if self._train:
            grad_input = np.multiply(grad_output, self.mask / (1 - self.p))
        else:
            grad_input = grad_output
        return grad_input
```

`BatchNorm` -- относительно современный слой, сильно улучшающий сходимость. Всё, что он делает -- это нормирует свои входные значения так, что на выходе получаются значения со средним 0 и дисперсией 1.

Почитать про вывод градиента для него можно тут: https://wiseodd.github.io/techblog/2016/07/04/batchnorm/

BatchNorm тоже по-разному ведёт себя при обучении и инференсе. Во время инференса он использует в качестве оценки среднего и дисперсии свои экспоненциально усреднённые исторические значения. Это связано с тем, что батч может быть маленьким, и оценки среднего и дисперсии будут неточными (при батче размера 1 дисперсия вообще будет нулевая, и нам в алгоритме нужно будет делить на ноль).

```python
class BatchNorm(Module):
    def __init__(self, num_features, eps=1e-05, momentum=0.1):
        super().__init__()
        self.eps = eps
        self.momentum = momentum
        self.gamma = np.ones(num_features)
        self.beta = np.zeros(num_features)
        self.moving_mean = 0
        self.moving_var = 0

    def forward(self, input):
        if self._train:
            sample_mean = input.mean(axis=0)
            sample_var = input.var(axis=0)

            self.moving_mean = sample_mean * self.momentum + self.moving_mean * (1 - self.momentum)
            self.moving_var = sample_var * self.momentum + self.moving_var * (1 - self.momentum)

            self.std = np.sqrt(sample_var + self.eps)
            self.centered = input - sample_mean
            self.normalized = self.centered / self.std
            self.output = self.gamma * self.normalized + self.beta

        else:
            input_norm = (input - self.moving_mean) / np.sqrt(self.moving_var + self.eps)
            self.output = self.gamma * input_norm + self.beta

        return self.output

    def backward(self, input, grad_output):
        if self._train:
            N = grad_output.shape[0]

            self.grad_gamma = (grad_output * self.normalized).sum(axis=0)
            self.grad_beta = grad_output.sum(axis=0)

            grad_input_norm = grad_output * self.gamma
            grad_input = 1/N / self.std * (N * grad_input_norm - grad_input_norm.sum(axis=0) -
            self.normalized * (grad_input_norm * grad_input_norm).sum(axis=0))

        else:
            grad_input = grad_output

        return grad_input

    def parameters(self):
        return [self.gamma, self.beta]

    def grad_parameters(self):
        return [self.grad_gamma, self.grad_beta]

```

    '\nclass BatchNorm(Module):\n    def __init__(self, num_features, gamma):\n        super().__init__()\n        self.gamma = gamma\n        self.mu = # ...\n        self.sigma = # ...\n    \n    def forward(self, input):\n        if self._train:\n            # ...\n        else:\n            # ...\n        return self.output\n    \n    def backward(self, input, grad_output):\n        if self._train:\n            # ...\n        else:\n            # ...\n        return grad_input\n'

## Критерии

Критерии — это специальные функции, которые меряют качество, имея реальные данные и предсказанные. Все критерии возвращают скаляр — одно число, усреднённое значение метрики по всему батчу.

По сути это тоже модули, но мы всё равно создадим для них отдельный класс, потому что у них нет `train` / `eval`, а `backward` не требует `grad_output` — эта вершина и так конечная в вычислительном графе. Также нам не понадобится сохранять для них `output`.

```python
class Criterion():
    def forward(self, input, target):
        raise NotImplementedError

    def backward(self, input, target):
        raise NotImplementedError
```

В качестве примера реализуем среднюю квадратичную ошибку (`MSE`).

Обратите внимание, что в критериях мы делим итоговое число на размер батча — мы не хотим, чтобы функция потерь зависела от количества примеров.

```python
class MSE(Criterion):
    def forward(self, input, target):
        batch_size = input.shape[0]
        self.output = np.sum(np.power(input - target, 2)) / input.shape[0]
        return self.output

    def backward(self, input, target):
        self.grad_output  = (input - target) * 2 / input.shape[0]
        return self.grad_output
```

Ваша задача посложнее: вам нужно реализовать кроссэнтропию — это стандартная функция потерь для классификации. Тут можно почитать про вывод её градиентов, а также софтмакса: https://deepnotes.io/softmax-crossentropy

Напоминаем интуицию за принципом максимального правдоподобия: мы максимизируем произведение предсказанных вероятностей реально случившихся событий $ L = \prod p_i $.

Произведение оптимизировать не очень удобно, и поэтому мы возьмём логарифм (любой, ведь все логарифмы отличаются в константу раз) и будем вместо него максимизировать сумму:

$$ \log L = \log \prod p_i = \sum \log p_i $$

Эту штуку называют кроссэнтропией. Такое название пошло из теории информации, но нам пока знать это не надо.

Для удобноства вместо чисел — от 0 до 9 — будем использовать вектора размера 10, где будет стоять единица в нужном месте (такое кодирование называется one-hot).

```python
class CrossEntropy(Criterion):
    def __init__(self):
        super().__init__()

    def forward(self, input, target):
        # To avoid log(0):
        eps = 1e-9
        input_clamp = np.clip(input, eps, 1 - eps)
        self.output = -1 * np.sum(np.log(input_clamp) * target) / input.shape[0]

        return self.output

    def backward(self, input, target):
        eps = 1e-9
        input_clamp = np.clip(input, eps, 1 - eps)

        grad_input = -1 * target / input_clamp
        return grad_input
```

max_stars_repo_name: Rexol/TinkoffDLcourse

max_issues_repo_name: Rexol/TinkoffDLcourse

max_forks_repo_name: Rexol/TinkoffDLcourse


text:
```python
! pip install qeds
import pandas as pd
import numpy as np
import qeds
```

    Requirement already satisfied: qeds in c:\users\asus\anaconda3\lib\site-packages (0.6.2)
    Requirement already satisfied: openpyxl in c:\users\asus\anaconda3\lib\site-packages (from qeds) (3.0.0)
    Requirement already satisfied: quandl in c:\users\asus\anaconda3\lib\site-packages (from qeds) (3.5.0)
    Requirement already satisfied: statsmodels in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.10.1)
    Requirement already satisfied: pandas in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.25.1)
    Requirement already satisfied: scikit-learn in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.21.3)
    Requirement already satisfied: numpy in c:\users\asus\anaconda3\lib\site-packages (from qeds) (1.16.5)
    Requirement already satisfied: pandas-datareader in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.8.1)
    Requirement already satisfied: requests in c:\users\asus\anaconda3\lib\site-packages (from qeds) (2.22.0)
    Requirement already satisfied: matplotlib in c:\users\asus\anaconda3\lib\site-packages (from qeds) (3.1.1)
    Requirement already satisfied: plotly in c:\users\asus\anaconda3\lib\site-packages (from qeds) (4.5.4)
    Requirement already satisfied: pyarrow in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.16.0)
    Requirement already satisfied: seaborn in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.9.0)
    Requirement already satisfied: scipy in c:\users\asus\anaconda3\lib\site-packages (from qeds) (1.3.1)
    Requirement already satisfied: quantecon in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.4.6)
    Requirement already satisfied: jdcal in c:\users\asus\anaconda3\lib\site-packages (from openpyxl->qeds) (1.4.1)
    Requirement already satisfied: et-xmlfile in c:\users\asus\anaconda3\lib\site-packages (from openpyxl->qeds) (1.0.1)
    Requirement already satisfied: more-itertools in c:\users\asus\anaconda3\lib\site-packages (from quandl->qeds) (7.2.0)
    Requirement already satisfied: six in c:\users\asus\anaconda3\lib\site-packages (from quandl->qeds) (1.12.0)
    Requirement already satisfied: python-dateutil in c:\users\asus\anaconda3\lib\site-packages (from quandl->qeds) (2.8.0)
    Requirement already satisfied: inflection>=0.3.1 in c:\users\asus\anaconda3\lib\site-packages (from quandl->qeds) (0.3.1)
    Requirement already satisfied: patsy>=0.4.0 in c:\users\asus\anaconda3\lib\site-packages (from statsmodels->qeds) (0.5.1)
    Requirement already satisfied: pytz>=2017.2 in c:\users\asus\anaconda3\lib\site-packages (from pandas->qeds) (2019.3)
    Requirement already satisfied: joblib>=0.11 in c:\users\asus\anaconda3\lib\site-packages (from scikit-learn->qeds) (0.13.2)
    Requirement already satisfied: lxml in c:\users\asus\anaconda3\lib\site-packages (from pandas-datareader->qeds) (4.4.1)
    Requirement already satisfied: chardet<3.1.0,>=3.0.2 in c:\users\asus\anaconda3\lib\site-packages (from requests->qeds) (3.0.4)
    Requirement already satisfied: urllib3!=1.25.0,!=1.25.1,<1.26,>=1.21.1 in c:\users\asus\anaconda3\lib\site-packages (from requests->qeds) (1.24.2)
    Requirement already satisfied: idna<2.9,>=2.5 in c:\users\asus\anaconda3\lib\site-packages (from requests->qeds) (2.8)
    Requirement already satisfied: certifi>=2017.4.17 in c:\users\asus\anaconda3\lib\site-packages (from requests->qeds) (2019.9.11)
    Requirement already satisfied: cycler>=0.10 in c:\users\asus\anaconda3\lib\site-packages (from matplotlib->qeds) (0.10.0)
    Requirement already satisfied: kiwisolver>=1.0.1 in c:\users\asus\anaconda3\lib\site-packages (from matplotlib->qeds) (1.1.0)
    Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in c:\users\asus\anaconda3\lib\site-packages (from matplotlib->qeds) (2.4.2)
    Requirement already satisfied: retrying>=1.3.3 in c:\users\asus\anaconda3\lib\site-packages (from plotly->qeds) (1.3.3)
    Requirement already satisfied: sympy in c:\users\asus\anaconda3\lib\site-packages (from quantecon->qeds) (1.4)
    Requirement already satisfied: numba>=0.38 in c:\users\asus\anaconda3\lib\site-packages (from quantecon->qeds) (0.45.1)
    Requirement already satisfied: setuptools in c:\users\asus\anaconda3\lib\site-packages (from kiwisolver>=1.0.1->matplotlib->qeds) (41.4.0)
    Requirement already satisfied: mpmath>=0.19 in c:\users\asus\anaconda3\lib\site-packages (from sympy->quantecon->qeds) (1.1.0)
    Requirement already satisfied: llvmlite>=0.29.0dev0 in c:\users\asus\anaconda3\lib\site-packages (from numba>=0.38->quantecon->qeds) (0.29.0)

# Cleaning Data - Exercises

## Exercise 1

**Convert the string below into a number.**

```python
c2n = "#39"

new_c2n = c2n.replace("#", "")

print(float(new_c2n))
```

    39.0

## Exercises 2-3

**For these exercises, we create the following DataFrame:**

```python
df = pd.DataFrame({"numbers": ["#23", "#24", "#18", "#14", "#12", "#10", "#35"],
                   "nums": ["23", "24", "18", "14", np.nan, "XYZ", "35"],
                   "colors": ["green", "red", "yellow", "orange", "purple", "blue", "pink"],
                   "other_column": [0, 1, 0, 2, 1, 0, 2]})
df.head()
```

      numbers
      nums
      colors
      other_column

      0
      #23
      23
      green
      0

      1
      #24
      24
      red
      1

      2
      #18
      18
      yellow
      0

      3
      #14
      14
      orange
      2

      4
      #12
      NaN
      purple
      1

## Exercise 2

   **Make a new column called `colors_upper` that contains the elements of `colors` with all uppercase letters.**

```python
df['colors_upper']=df['colors'].str.upper()
df.head()

#if we want only the first letter as capital letter
#df['colors_upper'] = df["colors"].str.capitalize()
```

      numbers
      nums
      colors
      other_column
      colors_upper

      0
      #23
      23
      green
      0
      GREEN

      1
      #24
      24
      red
      1
      RED

      2
      #18
      18
      yellow
      0
      YELLOW

      3
      #14
      14
      orange
      2
      ORANGE

      4
      #12
      NaN
      purple
      1
      PURPLE

## Exercise 3

1. **Convert the column `"nums"` to a numeric type using `pd.to_numeric` and save it to the DataFrame as `"nums_tonumeric"`.**

    - Notice that there is a missing value, and a value that is not a number. Look at the documentation for `pd.to_numeric` and think about how to overcome this.

```python
df["nums_tonumeric"] = pd.to_numeric(df["nums"], errors='coerce')
df
# errors "coerce" means that we force the conversion of string to NaN
```

      numbers
      nums
      colors
      other_column
      colors_upper
      nums_tonumeric

      0
      #23
      23
      green
      0
      GREEN
      23.0

      1
      #24
      24
      red
      1
      RED
      24.0

      2
      #18
      18
      yellow
      0
      YELLOW
      18.0

      3
      #14
      14
      orange
      2
      ORANGE
      14.0

      4
      #12
      NaN
      purple
      1
      PURPLE
      NaN

      5
      #10
      XYZ
      blue
      0
      BLUE
      NaN

      6
      #35
      35
      pink
      2
      PINK
      35.0

2. **Think about why this could be a bad idea of used without knowing what your data looks like. (Think about what happens when you apply it to the `"numbers"` column before replacing the `"#"`.)**

```python
#that would set all values to NaN
```

## Exercise 4

**For this exercise, we use data from an article written by The Upshot at the NYTimes, which has order information from almost 2,000 Chipotle orders and includes information on what was ordered and how much it cost.**

```python
chipotle = qeds.data.load("chipotle_raw")
chipotle.head()
```

      order_id
      quantity
      item_name
      choice_description
      item_price

      0
      1
      1
      Chips and Fresh Tomato Salsa
      NaN
      $2.39

      1
      1
      1
      Izze
      [Clementine]
      $3.39

      2
      1
      1
      Nantucket Nectar
      [Apple]
      $3.39

      3
      1
      1
      Chips and Tomatillo-Green Chili Salsa
      NaN
      $2.39

      4
      2
      2
      Chicken Bowl
      [Tomatillo-Red Chili Salsa (Hot), [Black Beans...
      $16.98

1. **Make sure the `item_price` column has a numeric `dtype` (probably float).**

```python
chipotle["item_price_numeric"] = pd.to_numeric(chipotle["item_price"].str.replace("$", ""))
chipotle.dtypes
chipotle.head()
```

      order_id
      quantity
      item_name
      choice_description
      item_price
      item_price_numeric

      0
      1
      1
      Chips and Fresh Tomato Salsa
      NaN
      $2.39
      2.39

      1
      1
      1
      Izze
      [Clementine]
      $3.39
      3.39

      2
      1
      1
      Nantucket Nectar
      [Apple]
      $3.39
      3.39

      3
      1
      1
      Chips and Tomatillo-Green Chili Salsa
      NaN
      $2.39
      2.39

      4
      2
      2
      Chicken Bowl
      [Tomatillo-Red Chili Salsa (Hot), [Black Beans...
      $16.98
      16.98

2. **What is the average price of an item with chicken?**

```python
chipotle.loc[chipotle['item_name'].str.match('Chicken'),"item_price_numeric"].mean()
```

    10.133724358974309

3. **What is the average price of an item with steak?**

```python
chipotle.loc[chipotle['item_name'].str.match('Steak'),"item_price_numeric"].mean()
```

    10.518888888888851

4. **Did chicken or steak produce more revenue (total)?**

```python
chipotle.loc[chipotle['item_name'].str.match('Chicken'),"item_price_numeric"].sum()>chipotle.loc[chipotle['item_name'].str.match('Steak'),"item_price_numeric"].sum()

#we create a boolean which confirms that chicken produces more total revenue compared to steak
```

    True

```python
#Alternative solution:
chipotle['revenue']=chipotle['item_price_numeric']*chipotle['quantity']

steak= chipotle[chipotle['item_name'].str.contains(pat='Steak')].sum()
print('Product containing steak generated a revenue of', steak.loc['revenue'], 'dollars')
steak= chipotle[chipotle['item_name'].str.contains(pat='Chicken')].sum()
print('Product containing chicken generated a revenue of', steak.loc['revenue'], 'dollars')
```

    Product containing steak generated a revenue of 8072.619999999973 dollars
    Product containing chicken generated a revenue of 17742.14999999992 dollars

5. **How many missing items are there in this dataset? How many missing items in each column?**

```python
#in the dataset
chipotle.isnull().sum().sum()
```

    1246

```python
#in each column
chipotle.isnull().sum()
```

    order_id                 0
    quantity                 0
    item_name                0
    choice_description    1246
    item_price               0
    item_price_numeric       0
    revenue                  0
    dtype: int64

max_stars_repo_name: remi-sudo/Classes

max_issues_repo_name: remi-sudo/Classes

max_forks_repo_name: remi-sudo/Classes


text:
```python
import numpy as np
import numba
import matplotlib.pyplot as plt
import sympy as sym
plt.style.use('presentation')
%matplotlib notebook

def d2np(d):

    names = []
    numbers = ()
    dtypes = []
    for item in d:
        names += item
        if type(d[item]) == float:
            numbers += (d[item],)
            dtypes += [(item,float)]
        if type(d[item]) == int:
            numbers += (d[item],)
            dtypes += [(item,int)]
        if type(d[item]) == np.ndarray:
            numbers += (d[item],)
            dtypes += [(item,np.float64,d[item].shape)]
    return np.array([numbers],dtype=dtypes)
```

```python
psi_ds,psi_qs,psi_dr,psi_qr = sym.symbols('psi_ds,psi_qs,psi_dr,psi_qr')
i_ds,i_qs,i_dr,i_qr = sym.symbols('i_ds,i_qs,i_dr,i_qr')
L_s,L_r,L_m = sym.symbols('L_s,L_r,L_m')

eq_ds = (L_s+L_m)*i_ds + L_m*i_dr - psi_ds
eq_qs = (L_s+L_m)*i_qs + L_m*i_qr - psi_qs
eq_dr = (L_r+L_m)*i_dr + L_m*i_ds - psi_dr
eq_qr = (L_r+L_m)*i_qr + L_m*i_qs - psi_qr

s = sym.solve([eq_ds,eq_qs,eq_dr,eq_qr],[i_ds,i_qs,i_dr,i_qr])

for item in s:
    print(item, '=', s[item])
```

    i_qs = (-L_m*psi_qr + L_m*psi_qs + L_r*psi_qs)/(L_m*L_r + L_m*L_s + L_r*L_s)
    i_dr = (L_m*psi_dr - L_m*psi_ds + L_s*psi_dr)/(L_m*L_r + L_m*L_s + L_r*L_s)
    i_ds = (-L_m*psi_dr + L_m*psi_ds + L_r*psi_ds)/(L_m*L_r + L_m*L_s + L_r*L_s)
    i_qr = (L_m*psi_qr - L_m*psi_qs + L_s*psi_qr)/(L_m*L_r + L_m*L_s + L_r*L_s)

```python
# [1] T. Demiray, F. Milano, and G. Andersson,
# “Dynamic phasor modeling of the doubly-fed induction generator under unbalanced conditions,” 2007 IEEE Lausanne POWERTECH, Proc., no. 2, pp. 1049–1054, 2007.

@numba.jit(nopython=True, cache=True)
def dfim(struct,i):

    x_idx = struct[i]['dfim_idx']

    L_m = struct[i]['L_m']
    L_r = struct[i]['L_r']
    L_s = struct[i]['L_s']
    R_r = struct[i]['R_r']
    R_s = struct[i]['R_s']

    Dt  = struct[i]['Dt']

    psi_ds = float(struct[i]['x'][x_idx+0,0])
    psi_qs = float(struct[i]['x'][x_idx+1,0])
    psi_dr = float(struct[i]['x'][x_idx+2,0])
    psi_qr = float(struct[i]['x'][x_idx+3,0])

    v_ds = struct[i]['v_ds']
    v_qs = struct[i]['v_qs']
    v_dr = struct[i]['v_dr']
    v_qr = struct[i]['v_qr']

    omega_r = struct[i]['omega_r']
    omega_s = struct[i]['omega_s']

    sigma = (omega_s - omega_r)/omega_s

    i_qr = (L_m*psi_qr - L_m*psi_qs + L_s*psi_qr)/(L_m*L_r + L_m*L_s + L_r*L_s)
    i_qs = (-L_m*psi_qr + L_m*psi_qs + L_r*psi_qs)/(L_m*L_r + L_m*L_s + L_r*L_s)
    i_dr = (L_m*psi_dr - L_m*psi_ds + L_s*psi_dr)/(L_m*L_r + L_m*L_s + L_r*L_s)
    i_ds = (-L_m*psi_dr + L_m*psi_ds + L_r*psi_ds)/(L_m*L_r + L_m*L_s + L_r*L_s)

    dpsi_ds = v_ds - R_s*i_ds + omega_s*psi_qs
    dpsi_qs = v_qs - R_s*i_qs - omega_s*psi_ds
    dpsi_dr = v_dr - R_r*i_dr + sigma*omega_s*psi_qr
    dpsi_qr = v_qr - R_r*i_qr - sigma*omega_s*psi_dr

    tau_e  = psi_qr*i_dr - psi_dr*i_qr

    struct[i]['f'][x_idx+0,0] = float(dpsi_ds)
    struct[i]['f'][x_idx+1,0] = float(dpsi_qs)
    struct[i]['f'][x_idx+2,0] = float(dpsi_dr)
    struct[i]['f'][x_idx+3,0] = float(dpsi_qr)

    struct[i]['i_ds'] = i_ds
    struct[i]['i_qs'] = i_qs
    struct[i]['i_dr'] = i_dr
    struct[i]['i_qr'] = i_qr

    struct[i]['psi_ds'] = psi_ds
    struct[i]['psi_qs'] = psi_qs
    struct[i]['psi_dr'] = psi_dr
    struct[i]['psi_qr'] = psi_qr

    struct[i]['tau_e'] = tau_e
    struct[i]['sigma'] = sigma

    return tau_e

@numba.jit(nopython=True, cache=True)
def wecs_mech_1(struct,i):

    x_idx = struct[i]['mech_idx']
    omega_t = struct[i]['x'][x_idx,0]  # rad/s
    tau_t   = struct[i]['tau_t']
    tau_r   = struct[i]['tau_r']

    J_t  = struct[i]['J_t']
    N_tr = struct[i]['N_tr']
    Dt   = struct[i]['Dt']

    domega_t = 1.0/J_t*(tau_t - N_tr*tau_r)

    omega_r = N_tr*omega_t

    struct[i]['f'][x_idx,0] = domega_t

    struct[i]['omega_r'] = omega_r
    struct[i]['omega_t'] = omega_t

    return omega_t

```

```python
Omega_b = 2.0*np.pi*50.0
S_b = 1.0e6
U_b = 690.0
Z_b = U_b**2/S_b
#nu_w =np.linspace(0.1,15,N)
H = 0.001

# H = 0.5*J*Omega_t_n**2/S_b
S_b = 2.0e6
Omega_t_n = 1.5
J_t = 2*H*S_b/Omega_t_n**2

#Z_b = 1.0
#Omega_b = 1.0
d =dict(R_r = 0.01*Z_b,
        R_s = 0.01*Z_b,
        L_r = 0.08*Z_b/Omega_b,
        L_s = 0.1*Z_b/Omega_b,
        L_m = 3.0*Z_b/Omega_b,
        psi_ds = 0.0,
        psi_qs = 0.0,
        psi_dr = 0.0,
        psi_qr = 0.0,
        i_ds   = 0.0,
        i_qs   = 0.0,
        i_dr   = 0.0,
        i_qr   = 0.0,
        v_ds   = 0.0,
        v_qs   = 0.0,
        v_dr   = 0.0,
        v_qr   = 0.0,
        omega_r   = Omega_b*0.99,
        omega_s   = Omega_b,
        sigma   = 0.0,
        tau_e   = 0.0,
        x = np.zeros((6,1)),
        f = np.zeros((6,1)),
        Dt = 0.0,
        J_t = J_t,
        omega_t = 0.0,
        tau_t = 0.0,
        tau_r = 0.0,
        N_tr = 20.0,
        dfim_idx = 0,
        mech_idx = 4,
        )

struct = d2np(d)
dfim(struct,0)
wecs_mech_1(struct,0)

```

    0.0

```python
struct = d2np(d)
struct['v_ds'] = 325.0
struct['v_qs'] = 0.0
struct['omega_r'] = Omega_b*1.0
Dt = 1.0e-4
struct[0]['x']= np.zeros((6,1))
struct[0]['x'][4,0] =   Omega_b*1.0
Tau_e = []
T =[]
N_steps = 100000
X = np.zeros((N_steps,6))

def f_eval(struct):
    dfim(struct,0)
    wecs_mech_1(struct,0)
    return struct[0]['f']

for it in range(N_steps):
    t = Dt*it
    if t>1.0:
        struct[0]['tau_t'] = 1.0e4
    dfim(struct,0)
    f1 = np.copy(f_eval(struct))
    x1 = np.copy(struct[0]['x'])
    struct[0]['x'] = np.copy(x1 + Dt*f1)
    dfim(struct,0)
    struct[0]['x'] = np.copy(x1 + 0.5*Dt*(f1 + f_eval(struct)))
    struct[0]['tau_r'] =  -struct[0]['tau_e']
    Tau_e += [float(struct['tau_e'])]
    X[it,:] = np.copy(struct[0]['x'][:]).T
    T +=[t]

```

```python
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(8, 5), sharex = True)
axes.plot(T,X[:,0])
axes.plot(T,X[:,1])
axes.plot(T,X[:,2])
axes.plot(T,X[:,3])

fig.savefig('dfim_tau_e.svg', bbox_inches='tight')

```

    <IPython.core.display.Javascript object>

```python
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(8, 5), sharex = True)
axes.plot(T,Tau_e)

fig.savefig('dfim_tau_e.svg', bbox_inches='tight')

```

    <IPython.core.display.Javascript object>

```python
fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(8, 5), sharex = True)
axes.plot(T,X[:,4])
```

    <IPython.core.display.Javascript object>

    [<matplotlib.lines.Line2D at 0x7fa1a7d0d7b8>]

```python
X[:,4]
```

    array([ 314.15926536,  314.15926536,  314.15926535, ...,   15.82717562,
             15.82717562,   15.82717562])

```python
20*15.68
```

    313.6

```python

```

```python

```

```python

```

max_stars_repo_name: pydgrid/pydgrid

max_issues_repo_name: pydgrid/pydgrid

max_forks_repo_name: pydgrid/pydgrid


text:
# Timeseries prediction of the Mackey-Glass Equation with ESNs

## Introduction

The Mackey-Glass system is essentially the differential equation, where we set the parameters to $\alpha = 0.2$, $\beta = 10$, $\gamma = 0.1$ and the time delay $\tau = 17$ in  order to have a mildly chaotic attractor.

\begin{align}
\label{eq:MackeyGlass}
\dot{y}(t) = \alpha y(t-\tau) / (1 + y(t - \tau)^{\beta}) - \gamma y(t)
\end{align}

```python
import numpy as np
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error, make_scorer
from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import TimeSeriesSplit, GridSearchCV

from matplotlib import pyplot as plt
plt.rcParams["font.family"] = "Times New Roman"
plt.rcParams["font.size"] = 10

from pyrcn.echo_state_network import ESNRegressor
from pyrcn.extreme_learning_machine import ELMRegressor
from pyrcn.linear_model import IncrementalRegression
from pyrcn.model_selection import SequentialSearchCV
from pyrcn.datasets import mackey_glass
```

Load the dataset and rescale it to a range of [-1, 1]

```python
# Load the dataset
X, y = mackey_glass(n_timesteps=20000)
scaler = MinMaxScaler(feature_range=(-1, 1)).fit(X=X)
X = scaler.transform(X=X)
y = scaler.transform(y)
```

Define Train/Test lengths

```python
trainLen = 1900 # number of time steps during which we train the network
testLen = 2000 # number of time steps during which we test/run the network

X_train = X[:trainLen]
y_train = y[:trainLen]
X_test = X[trainLen:trainLen+testLen]
y_test = y[trainLen:trainLen+testLen]
```

Visualization

```python
fig = plt.figure()
plt.plot(X_train, label="Training input")
plt.plot(y_train, label="Training target")
plt.xlabel("n")
plt.xlim([0, 200])
plt.ylabel("u[n]")
plt.grid()
plt.legend()
fig.set_size_inches(4, 2.5)
```

Training and Prediction using vanilla ESNs and ELMs

```python
# initialize an ESNRegressor
esn = ESNRegressor()  # IncrementalRegression()

# initialize an ELMRegressor
elm = ELMRegressor(regressor=Ridge())  # Ridge()

# train a model
esn.fit(X=X_train, y=y_train)
elm.fit(X=X_train, y=y_train)

# evaluate the models
y_test_pred = esn.predict(X=X_test)
print(mean_squared_error(y_test, y_test_pred))
y_test_pred = elm.predict(X=X_test)
print(mean_squared_error(y_test, y_test_pred))
```

Hyperparameter optimization ESN

```python
# Echo State Network sequential hyperparameter tuning
initially_fixed_esn_params = {'hidden_layer_size': 100,
                              'input_activation': 'identity',
                              'bias_scaling': 0.0,
                              'random_state': 42,
                              'k_in': 1,
                              'reservoir_activation': 'tanh',
                              'leakage': 1.0,
                              'bi_directional': False,
                              'k_rec': 10,
                              'wash_out': 0,
                              'continuation': False,
                              'alpha': 1e-5 }

step1_esn_params = {'input_scaling': np.linspace(0.1, 5.0, 50),
                    'spectral_radius': np.linspace(0.0, 1.5, 16)}
step2_esn_params = {'leakage': np.linspace(0.1, 1.0, 10)}
step3_esn_params = {'bias_scaling': np.linspace(0.0, 1.5, 16)}

scorer = make_scorer(score_func=mean_squared_error, greater_is_better=False)

kwargs = {'verbose': 5,
          'scoring': scorer,
          'n_jobs': -1,
          'cv': TimeSeriesSplit()}

esn = ESNRegressor(regressor=Ridge(), **initially_fixed_esn_params)

searches = [('step1', GridSearchCV, step1_esn_params, kwargs),
            ('step2', GridSearchCV, step2_esn_params, kwargs),
            ('step3', GridSearchCV, step3_esn_params, kwargs)]

sequential_search_esn = SequentialSearchCV(esn, searches=searches).fit(X_train, y_train)
```

Hyperparameter optimization ELM

```python
# Extreme Learning Machine sequential hyperparameter tuning
initially_fixed_elm_params = {'hidden_layer_size': 100,
                              'activation': 'tanh',
                              'k_in': 1,
                              'alpha': 1e-5,
                              'random_state': 42 }

step1_elm_params = {'input_scaling': np.linspace(0.1, 5.0, 50)}
step2_elm_params = {'bias_scaling': np.linspace(0.0, 1.5, 16)}

scorer = make_scorer(score_func=mean_squared_error, greater_is_better=False)

kwargs = {'verbose': 5,
          'scoring': scorer,
          'n_jobs': -1,
          'cv': TimeSeriesSplit()}

elm = ELMRegressor(regressor=Ridge(), **initially_fixed_elm_params)

searches = [('step1', GridSearchCV, step1_elm_params, kwargs),
            ('step2', GridSearchCV, step2_elm_params, kwargs)]

sequential_search_elm = SequentialSearchCV(elm, searches=searches).fit(X_train, y_train)
```

Final prediction and visualization

```python
sequential_search_esn.all_best_score_
```

```python
sequential_search_elm.all_best_score_
```

```python
esn = sequential_search_esn.best_estimator_
elm = sequential_search_elm.best_estimator_

y_train_pred_esn = esn.predict(X=X_train)
y_train_pred_elm = elm.predict(X=X_train)
y_test_pred_esn = esn.predict(X=X_test)
y_test_pred_elm = elm.predict(X=X_test)

test_err_esn = mean_squared_error(y_true=y_test, y_pred=y_test_pred_esn)
test_err_elm = mean_squared_error(y_true=y_test, y_pred=y_test_pred_elm)

print("Test MSE ESN:\t{0}".format(test_err_esn))
print("Test MSE ELM:\t{0}".format(test_err_elm))

# Prediction of the test set.
fig = plt.figure()
plt.plot(y_test_pred_esn, label="ESN prediction", linewidth=1)
plt.plot(y_test_pred_elm, label="ELM prediction", linewidth=1)
plt.plot(y_test, label="Test target", linewidth=.5, color="black")
plt.xlabel("n")
plt.xlim([0, 100])
plt.ylabel("u[n]")
plt.grid()
plt.legend()
fig.set_size_inches(4, 2.5)
plt.savefig("test_data.pdf")
```

```python
fig = plt.figure()
im = plt.imshow(np.abs(elm.hidden_layer_state[:100, 1:].T),vmin=0, vmax=1.0)
plt.xlim([0, 100])
plt.ylim([0, elm.hidden_layer_state.shape[1] - 1])
plt.xlabel('n')
plt.ylabel('R\'[n]')
plt.colorbar(im)
plt.grid()
fig.set_size_inches(2, 1.25)
plt.savefig("elm_reservoir_state.pdf", bbox_inches = 'tight', pad_inches = 0)
```

```python
fig = plt.figure()
im = plt.imshow(np.abs(esn.hidden_layer_state[:100, 1:].T),vmin=0, vmax=1.0)
plt.xlim([0, 100])
plt.ylim([0, esn.hidden_layer_state.shape[1] - 1])
plt.xlabel('n')
plt.ylabel('R[n]')
plt.colorbar(im)
plt.grid()
fig.set_size_inches(2, 1.25)
plt.savefig("esn_reservoir_state.pdf", bbox_inches = 'tight', pad_inches = 0)
```

```python

```

max_stars_repo_name: TUD-STKS/PyRCN

max_issues_repo_name: TUD-STKS/PyRCN

max_forks_repo_name: TUD-STKS/PyRCN