text:
# Batch Normalization
One way to make deep networks easier to train is to use more sophisticated optimization procedures such as SGD+momentum, RMSProp, or Adam. Another strategy is to change the architecture of the network to make it easier to train.
One idea along these lines is batch normalization which was proposed by [1] in 2015.

The idea is relatively straightforward. Machine learning methods tend to work better when their input data consists of uncorrelated features with zero mean and unit variance. When training a neural network, we can preprocess the data before feeding it to the network to explicitly decorrelate its features; this will ensure that the first layer of the network sees data that follows a nice distribution. However, even if we preprocess the input data, the activations at deeper layers of the network will likely no longer be decorrelated and will no longer have zero mean or unit variance since they are output from earlier layers in the network. Even worse, during the training process the distribution of features at each layer of the network will shift as the weights of each layer are updated.

The authors of [1] hypothesize that the shifting distribution of features inside deep neural networks may make training deep networks more difficult. To overcome this problem, [1] proposes to insert batch normalization layers into the network. At training time, a batch normalization layer uses a minibatch of data to estimate the mean and standard deviation of each feature. These estimated means and standard deviations are then used to center and normalize the features of the minibatch. A running average of these means and standard deviations is kept during training, and at test time these running averages are used to center and normalize features.

It is possible that this normalization strategy could reduce the representational power of the network, since it may sometimes be optimal for certain layers to have features that are not zero-mean or unit variance. To this end, the batch normalization layer includes learnable shift and scale parameters for each feature dimension.

[1] [Sergey Ioffe and Christian Szegedy, "Batch Normalization: Accelerating Deep Network Training by Reducing
Internal Covariate Shift", ICML 2015.](https://arxiv.org/abs/1502.03167)

```python
# As usual, a bit of setup
import time
import numpy as np
import matplotlib.pyplot as plt
from cs231n.classifiers.fc_net import *
from cs231n.data_utils import get_CIFAR10_data
from cs231n.gradient_check import eval_numerical_gradient, eval_numerical_gradient_array
from cs231n.solver import Solver

%matplotlib inline
plt.rcParams['figure.figsize'] = (10.0, 8.0) # set default size of plots
plt.rcParams['image.interpolation'] = 'nearest'
plt.rcParams['image.cmap'] = 'gray'

# for auto-reloading external modules
# see http://stackoverflow.com/questions/1907993/autoreload-of-modules-in-ipython
%load_ext autoreload
%autoreload 2

def rel_error(x, y):
    """ returns relative error """
    return np.max(np.abs(x - y) / (np.maximum(1e-8, np.abs(x) + np.abs(y))))

def print_mean_std(x,axis=0):
    print('  means: ', x.mean(axis=axis))
    print('  stds:  ', x.std(axis=axis))
    print()
```

    run the following from the cs231n directory and try again:
    python setup.py build_ext --inplace
    You may also need to restart your iPython kernel

```python
# Load the (preprocessed) CIFAR10 data.
data = get_CIFAR10_data()
for k, v in data.items():
  print('%s: ' % k, v.shape)
```

    X_train:  (49000, 3, 32, 32)
    y_train:  (49000,)
    X_val:  (1000, 3, 32, 32)
    y_val:  (1000,)
    X_test:  (1000, 3, 32, 32)
    y_test:  (1000,)

## Batch normalization: forward
In the file `cs231n/layers.py`, implement the batch normalization forward pass in the function `batchnorm_forward`. Once you have done so, run the following to test your implementation.

Referencing the paper linked to above in [1] may be helpful!

```python
# Check the training-time forward pass by checking means and variances
# of features both before and after batch normalization

# Simulate the forward pass for a two-layer network
np.random.seed(231)
N, D1, D2, D3 = 200, 50, 60, 3
X = np.random.randn(N, D1)
W1 = np.random.randn(D1, D2)
W2 = np.random.randn(D2, D3)
a = np.maximum(0, X.dot(W1)).dot(W2)

print('Before batch normalization:')
print_mean_std(a,axis=0)

gamma = np.ones((D3,))
beta = np.zeros((D3,))
# Means should be close to zero and stds close to one
print('After batch normalization (gamma=1, beta=0)')
a_norm, _ = batchnorm_forward(a, gamma, beta, {'mode': 'train'})
print_mean_std(a_norm,axis=0)

gamma = np.asarray([1.0, 2.0, 3.0])
beta = np.asarray([11.0, 12.0, 13.0])
# Now means should be close to beta and stds close to gamma
print('After batch normalization (gamma=', gamma, ', beta=', beta, ')')
a_norm, _ = batchnorm_forward(a, gamma, beta, {'mode': 'train'})
print_mean_std(a_norm,axis=0)
```

    Before batch normalization:
      means:  [ -2.3814598  -13.18038246   1.91780462]
      stds:   [27.18502186 34.21455511 37.68611762]

    After batch normalization (gamma=1, beta=0)
      means:  [5.99520433e-17 6.93889390e-17 8.32667268e-19]
      stds:   [0.99999999 1.         1.        ]

    After batch normalization (gamma= [1. 2. 3.] , beta= [11. 12. 13.] )
      means:  [11. 12. 13.]
      stds:   [0.99999999 1.99999999 2.99999999]

```python
# Check the test-time forward pass by running the training-time
# forward pass many times to warm up the running averages, and then
# checking the means and variances of activations after a test-time
# forward pass.

np.random.seed(231)
N, D1, D2, D3 = 200, 50, 60, 3
W1 = np.random.randn(D1, D2)
W2 = np.random.randn(D2, D3)

bn_param = {'mode': 'train'}
gamma = np.ones(D3)
beta = np.zeros(D3)

for t in range(50):
  X = np.random.randn(N, D1)
  a = np.maximum(0, X.dot(W1)).dot(W2)
  batchnorm_forward(a, gamma, beta, bn_param)

bn_param['mode'] = 'test'
X = np.random.randn(N, D1)
a = np.maximum(0, X.dot(W1)).dot(W2)
a_norm, _ = batchnorm_forward(a, gamma, beta, bn_param)

# Means should be close to zero and stds close to one, but will be
# noisier than training-time forward passes.
print('After batch normalization (test-time):')
print_mean_std(a_norm,axis=0)
```

    After batch normalization (test-time):
      means:  [-0.03927354 -0.04349152 -0.10452688]
      stds:   [1.01531428 1.01238373 0.97819988]

## Batch normalization: backward
Now implement the backward pass for batch normalization in the function `batchnorm_backward`.

To derive the backward pass you should write out the computation graph for batch normalization and backprop through each of the intermediate nodes. Some intermediates may have multiple outgoing branches; make sure to sum gradients across these branches in the backward pass.

Once you have finished, run the following to numerically check your backward pass.

```python
# Gradient check batchnorm backward pass
np.random.seed(231)
N, D = 4, 5
x = 5 * np.random.randn(N, D) + 12
gamma = np.random.randn(D)
beta = np.random.randn(D)
dout = np.random.randn(N, D)

bn_param = {'mode': 'train'}
fx = lambda x: batchnorm_forward(x, gamma, beta, bn_param)[0]
fg = lambda a: batchnorm_forward(x, a, beta, bn_param)[0]
fb = lambda b: batchnorm_forward(x, gamma, b, bn_param)[0]

dx_num = eval_numerical_gradient_array(fx, x, dout)
da_num = eval_numerical_gradient_array(fg, gamma.copy(), dout)
db_num = eval_numerical_gradient_array(fb, beta.copy(), dout)

_, cache = batchnorm_forward(x, gamma, beta, bn_param)
dx, dgamma, dbeta = batchnorm_backward(dout, cache)
#You should expect to see relative errors between 1e-13 and 1e-8
print('dx error: ', rel_error(dx_num, dx))
print('dgamma error: ', rel_error(da_num, dgamma))
print('dbeta error: ', rel_error(db_num, dbeta))
```

    dx error:  1.6674604875341426e-09
    dgamma error:  7.417225040694815e-13
    dbeta error:  2.379446949959628e-12

## Batch normalization: alternative backward
In class we talked about two different implementations for the sigmoid backward pass. One strategy is to write out a computation graph composed of simple operations and backprop through all intermediate values. Another strategy is to work out the derivatives on paper. For example, you can derive a very simple formula for the sigmoid function's backward pass by simplifying gradients on paper.

Surprisingly, it turns out that you can do a similar simplification for the batch normalization backward pass too!

In the forward pass, given a set of inputs $X=\begin{bmatrix}x_1\\x_2\\...\\x_N\end{bmatrix}$,

we first calculate the mean $\mu$ and variance $v$.
With $\mu$ and $v$ calculated, we can calculate the standard deviation $\sigma$  and normalized data $Y$.
The equations and graph illustration below describe the computation ($y_i$ is the i-th element of the vector $Y$).

\begin{align}
& \mu=\frac{1}{N}\sum_{k=1}^N x_k  &  v=\frac{1}{N}\sum_{k=1}^N (x_k-\mu)^2 \\
& \sigma=\sqrt{v+\epsilon}         &  y_i=\frac{x_i-\mu}{\sigma}
\end{align}

The meat of our problem during backpropagation is to compute $\frac{\partial L}{\partial X}$, given the upstream gradient we receive, $\frac{\partial L}{\partial Y}.$ To do this, recall the chain rule in calculus gives us $\frac{\partial L}{\partial X} = \frac{\partial L}{\partial Y} \cdot \frac{\partial Y}{\partial X}$.

The unknown/hart part is $\frac{\partial Y}{\partial X}$. We can find this by first deriving step-by-step our local gradients at
$\frac{\partial v}{\partial X}$, $\frac{\partial \mu}{\partial X}$,
$\frac{\partial \sigma}{\partial v}$,
$\frac{\partial Y}{\partial \sigma}$, and $\frac{\partial Y}{\partial \mu}$,
and then use the chain rule to compose these gradients (which appear in the form of vectors!) appropriately to compute $\frac{\partial Y}{\partial X}$.

If it's challenging to directly reason about the gradients over $X$ and $Y$ which require matrix multiplication, try reasoning about the gradients in terms of individual elements $x_i$ and $y_i$ first: in that case, you will need to come up with the derivations for $\frac{\partial L}{\partial x_i}$, by relying on the Chain Rule to first calculate the intermediate $\frac{\partial \mu}{\partial x_i}, \frac{\partial v}{\partial x_i}, \frac{\partial \sigma}{\partial x_i},$ then assemble these pieces to calculate $\frac{\partial y_i}{\partial x_i}$.

You should make sure each of the intermediary gradient derivations are all as simplified as possible, for ease of implementation.

After doing so, implement the simplified batch normalization backward pass in the function `batchnorm_backward_alt` and compare the two implementations by running the following. Your two implementations should compute nearly identical results, but the alternative implementation should be a bit faster.

```python
np.random.seed(231)
N, D = 100, 500
x = 5 * np.random.randn(N, D) + 12
gamma = np.random.randn(D)
beta = np.random.randn(D)
dout = np.random.randn(N, D)

bn_param = {'mode': 'train'}
out, cache = batchnorm_forward(x, gamma, beta, bn_param)

t1 = time.time()
dx1, dgamma1, dbeta1 = batchnorm_backward(dout, cache)
t2 = time.time()
dx2, dgamma2, dbeta2 = batchnorm_backward_alt(dout, cache)
t3 = time.time()

print('dx difference: ', rel_error(dx1, dx2))
print('dgamma difference: ', rel_error(dgamma1, dgamma2))
print('dbeta difference: ', rel_error(dbeta1, dbeta2))
try:
    print('speedup: %.2fx' % ((t2 - t1) / (t3 - t2)))
except ZeroDivisionError:
    pass
```

    dx difference:  7.666703093523215e-13
    dgamma difference:  0.0
    dbeta difference:  0.0
    speedup: 1.92x

## Fully Connected Nets with Batch Normalization
Now that you have a working implementation for batch normalization, go back to your `FullyConnectedNet` in the file `cs231n/classifiers/fc_net.py`. Modify your implementation to add batch normalization.

Concretely, when the `normalization` flag is set to `"batchnorm"` in the constructor, you should insert a batch normalization layer before each ReLU nonlinearity. The outputs from the last layer of the network should not be normalized. Once you are done, run the following to gradient-check your implementation.

HINT: You might find it useful to define an additional helper layer similar to those in the file `cs231n/layer_utils.py`. If you decide to do so, do it in the file `cs231n/classifiers/fc_net.py`.

```python
np.random.seed(231)
N, D, H1, H2, C = 2, 15, 20, 30, 10
X = np.random.randn(N, D)
y = np.random.randint(C, size=(N,))

# You should expect losses between 1e-4~1e-10 for W,
# losses between 1e-08~1e-10 for b,
# and losses between 1e-08~1e-09 for beta and gammas.
for reg in [0, 3.14]:
  print('Running check with reg = ', reg)
  model = FullyConnectedNet([H1, H2], input_dim=D, num_classes=C,
                            reg=reg, weight_scale=5e-2, dtype=np.float64,
                            normalization='batchnorm')

  loss, grads = model.loss(X, y)
  print('Initial loss: ', loss)

  for name in sorted(grads):
    f = lambda _: model.loss(X, y)[0]
    grad_num = eval_numerical_gradient(f, model.params[name], verbose=False, h=1e-5)
    print('%s relative error: %.2e' % (name, rel_error(grad_num, grads[name])))
  if reg == 0: print()
```

    Running check with reg =  0
    Initial loss:  2.2611955101340957
    W1 relative error: 1.10e-04
    W2 relative error: 3.11e-06
    W3 relative error: 4.05e-10
    b1 relative error: 5.55e-09
    b2 relative error: 6.66e-08
    b3 relative error: 1.01e-10
    beta1 relative error: 7.33e-09
    beta2 relative error: 1.89e-09
    gamma1 relative error: 6.96e-09
    gamma2 relative error: 2.41e-09

    Running check with reg =  3.14
    Initial loss:  6.996533220108303
    W1 relative error: 1.98e-06
    W2 relative error: 2.29e-06
    W3 relative error: 2.79e-08
    b1 relative error: 5.55e-09
    b2 relative error: 4.44e-08
    b3 relative error: 2.10e-10
    beta1 relative error: 6.65e-09
    beta2 relative error: 3.39e-09
    gamma1 relative error: 6.27e-09
    gamma2 relative error: 5.28e-09

# Batchnorm for deep networks
Run the following to train a six-layer network on a subset of 1000 training examples both with and without batch normalization.

```python
np.random.seed(231)
# Try training a very deep net with batchnorm
hidden_dims = [100, 100, 100, 100, 100]

num_train = 1000
small_data = {
  'X_train': data['X_train'][:num_train],
  'y_train': data['y_train'][:num_train],
  'X_val': data['X_val'],
  'y_val': data['y_val'],
}

weight_scale = 2e-2
bn_model = FullyConnectedNet(hidden_dims, weight_scale=weight_scale, normalization='batchnorm')
model = FullyConnectedNet(hidden_dims, weight_scale=weight_scale, normalization=None)

print('Solver with batch norm:')
bn_solver = Solver(bn_model, small_data,
                num_epochs=10, batch_size=50,
                update_rule='adam',
                optim_config={
                  'learning_rate': 1e-3,
                },
                verbose=True,print_every=20)
bn_solver.train()

print('\nSolver without batch norm:')
solver = Solver(model, small_data,
                num_epochs=10, batch_size=50,
                update_rule='adam',
                optim_config={
                  'learning_rate': 1e-3,
                },
                verbose=True, print_every=20)
solver.train()
```

    Solver with batch norm:
    (Iteration 1 / 200) loss: 2.297074
    (Epoch 0 / 10) train acc: 0.102000; val_acc: 0.107000
    (Epoch 1 / 10) train acc: 0.354000; val_acc: 0.279000
    (Iteration 21 / 200) loss: 1.979541
    (Epoch 2 / 10) train acc: 0.434000; val_acc: 0.301000
    (Iteration 41 / 200) loss: 1.776136
    (Epoch 3 / 10) train acc: 0.489000; val_acc: 0.320000
    (Iteration 61 / 200) loss: 1.638009
    (Epoch 4 / 10) train acc: 0.494000; val_acc: 0.294000
    (Iteration 81 / 200) loss: 1.378472
    (Epoch 5 / 10) train acc: 0.582000; val_acc: 0.313000
    (Iteration 101 / 200) loss: 1.476406
    (Epoch 6 / 10) train acc: 0.627000; val_acc: 0.341000
    (Iteration 121 / 200) loss: 1.053954
    (Epoch 7 / 10) train acc: 0.696000; val_acc: 0.324000
    (Iteration 141 / 200) loss: 1.200633
    (Epoch 8 / 10) train acc: 0.729000; val_acc: 0.330000
    (Iteration 161 / 200) loss: 0.884257
    (Epoch 9 / 10) train acc: 0.766000; val_acc: 0.321000
    (Iteration 181 / 200) loss: 0.750285
    (Epoch 10 / 10) train acc: 0.805000; val_acc: 0.299000

    Solver without batch norm:
    (Iteration 1 / 200) loss: 2.302600
    (Epoch 0 / 10) train acc: 0.124000; val_acc: 0.106000
    (Epoch 1 / 10) train acc: 0.185000; val_acc: 0.168000
    (Iteration 21 / 200) loss: 2.248002
    (Epoch 2 / 10) train acc: 0.220000; val_acc: 0.226000
    (Iteration 41 / 200) loss: 1.977625
    (Epoch 3 / 10) train acc: 0.257000; val_acc: 0.218000
    (Iteration 61 / 200) loss: 2.037953
    (Epoch 4 / 10) train acc: 0.244000; val_acc: 0.227000
    (Iteration 81 / 200) loss: 2.000312
    (Epoch 5 / 10) train acc: 0.299000; val_acc: 0.237000
    (Iteration 101 / 200) loss: 1.748184
    (Epoch 6 / 10) train acc: 0.362000; val_acc: 0.279000
    (Iteration 121 / 200) loss: 1.724046
    (Epoch 7 / 10) train acc: 0.342000; val_acc: 0.258000
    (Iteration 141 / 200) loss: 1.696783
    (Epoch 8 / 10) train acc: 0.366000; val_acc: 0.252000
    (Iteration 161 / 200) loss: 1.539274
    (Epoch 9 / 10) train acc: 0.426000; val_acc: 0.256000
    (Iteration 181 / 200) loss: 1.462405
    (Epoch 10 / 10) train acc: 0.404000; val_acc: 0.273000

Run the following to visualize the results from two networks trained above. You should find that using batch normalization helps the network to converge much faster.

```python
def plot_training_history(title, label, baseline, bn_solvers, plot_fn, bl_marker='.', bn_marker='.', labels=None):
    """utility function for plotting training history"""
    plt.title(title)
    plt.xlabel(label)
    bn_plots = [plot_fn(bn_solver) for bn_solver in bn_solvers]
    bl_plot = plot_fn(baseline)
    num_bn = len(bn_plots)
    for i in range(num_bn):
        label='with_norm'
        if labels is not None:
            label += str(labels[i])
        plt.plot(bn_plots[i], bn_marker, label=label)
    label='baseline'
    if labels is not None:
        label += str(labels[0])
    plt.plot(bl_plot, bl_marker, label=label)
    plt.legend(loc='lower center', ncol=num_bn+1)

plt.subplot(3, 1, 1)
plot_training_history('Training loss','Iteration', solver, [bn_solver], \
                      lambda x: x.loss_history, bl_marker='o', bn_marker='o')
plt.subplot(3, 1, 2)
plot_training_history('Training accuracy','Epoch', solver, [bn_solver], \
                      lambda x: x.train_acc_history, bl_marker='-o', bn_marker='-o')
plt.subplot(3, 1, 3)
plot_training_history('Validation accuracy','Epoch', solver, [bn_solver], \
                      lambda x: x.val_acc_history, bl_marker='-o', bn_marker='-o')

plt.gcf().set_size_inches(15, 15)
plt.show()
```

# Batch normalization and initialization
We will now run a small experiment to study the interaction of batch normalization and weight initialization.

The first cell will train 8-layer networks both with and without batch normalization using different scales for weight initialization. The second layer will plot training accuracy, validation set accuracy, and training loss as a function of the weight initialization scale.

```python
np.random.seed(231)
# Try training a very deep net with batchnorm
hidden_dims = [50, 50, 50, 50, 50, 50, 50]
num_train = 1000
small_data = {
  'X_train': data['X_train'][:num_train],
  'y_train': data['y_train'][:num_train],
  'X_val': data['X_val'],
  'y_val': data['y_val'],
}

bn_solvers_ws = {}
solvers_ws = {}
weight_scales = np.logspace(-4, 0, num=20)
for i, weight_scale in enumerate(weight_scales):
  print('Running weight scale %d / %d' % (i + 1, len(weight_scales)))
  bn_model = FullyConnectedNet(hidden_dims, weight_scale=weight_scale, normalization='batchnorm')
  model = FullyConnectedNet(hidden_dims, weight_scale=weight_scale, normalization=None)

  bn_solver = Solver(bn_model, small_data,
                  num_epochs=10, batch_size=50,
                  update_rule='adam',
                  optim_config={
                    'learning_rate': 1e-3,
                  },
                  verbose=False, print_every=200)
  bn_solver.train()
  bn_solvers_ws[weight_scale] = bn_solver

  solver = Solver(model, small_data,
                  num_epochs=10, batch_size=50,
                  update_rule='adam',
                  optim_config={
                    'learning_rate': 1e-3,
                  },
                  verbose=False, print_every=200)
  solver.train()
  solvers_ws[weight_scale] = solver
```

    Running weight scale 1 / 20
    Running weight scale 2 / 20
    Running weight scale 3 / 20
    Running weight scale 4 / 20
    Running weight scale 5 / 20
    Running weight scale 6 / 20
    Running weight scale 7 / 20
    Running weight scale 8 / 20
    Running weight scale 9 / 20
    Running weight scale 10 / 20
    Running weight scale 11 / 20
    Running weight scale 12 / 20
    Running weight scale 13 / 20
    Running weight scale 14 / 20
    Running weight scale 15 / 20
    Running weight scale 16 / 20
    Running weight scale 17 / 20
    Running weight scale 18 / 20
    Running weight scale 19 / 20
    Running weight scale 20 / 20

```python
# Plot results of weight scale experiment
best_train_accs, bn_best_train_accs = [], []
best_val_accs, bn_best_val_accs = [], []
final_train_loss, bn_final_train_loss = [], []

for ws in weight_scales:
  best_train_accs.append(max(solvers_ws[ws].train_acc_history))
  bn_best_train_accs.append(max(bn_solvers_ws[ws].train_acc_history))

  best_val_accs.append(max(solvers_ws[ws].val_acc_history))
  bn_best_val_accs.append(max(bn_solvers_ws[ws].val_acc_history))

  final_train_loss.append(np.mean(solvers_ws[ws].loss_history[-100:]))
  bn_final_train_loss.append(np.mean(bn_solvers_ws[ws].loss_history[-100:]))

plt.subplot(3, 1, 1)
plt.title('Best val accuracy vs weight initialization scale')
plt.xlabel('Weight initialization scale')
plt.ylabel('Best val accuracy')
plt.semilogx(weight_scales, best_val_accs, '-o', label='baseline')
plt.semilogx(weight_scales, bn_best_val_accs, '-o', label='batchnorm')
plt.legend(ncol=2, loc='lower right')

plt.subplot(3, 1, 2)
plt.title('Best train accuracy vs weight initialization scale')
plt.xlabel('Weight initialization scale')
plt.ylabel('Best training accuracy')
plt.semilogx(weight_scales, best_train_accs, '-o', label='baseline')
plt.semilogx(weight_scales, bn_best_train_accs, '-o', label='batchnorm')
plt.legend()

plt.subplot(3, 1, 3)
plt.title('Final training loss vs weight initialization scale')
plt.xlabel('Weight initialization scale')
plt.ylabel('Final training loss')
plt.semilogx(weight_scales, final_train_loss, '-o', label='baseline')
plt.semilogx(weight_scales, bn_final_train_loss, '-o', label='batchnorm')
plt.legend()
plt.gca().set_ylim(1.0, 3.5)

plt.gcf().set_size_inches(15, 15)
plt.show()
```

## Inline Question 1:
Describe the results of this experiment. How does the scale of weight initialization affect models with/without batch normalization differently, and why?

## Answer:
The models without BN effectively don't learn over a large range of weight initializations.  The BN models fairly consistently outperforms ther non BN models.  This is again due to vanishing/exploding gradients which can arise due to poor initialization causing the gradient update to go to zero or inf.  By enforcing zero mean and unit variance (with a learned shift and scaling factor) the model at each layer is able to compensate and not cause earlier layers to stop learning effectively.

# Batch normalization and batch size
We will now run a small experiment to study the interaction of batch normalization and batch size.

The first cell will train 6-layer networks both with and without batch normalization using different batch sizes. The second layer will plot training accuracy and validation set accuracy over time.

```python
def run_batchsize_experiments(normalization_mode):
    np.random.seed(231)
    # Try training a very deep net with batchnorm
    hidden_dims = [100, 100, 100, 100, 100]
    num_train = 1000
    small_data = {
      'X_train': data['X_train'][:num_train],
      'y_train': data['y_train'][:num_train],
      'X_val': data['X_val'],
      'y_val': data['y_val'],
    }
    n_epochs=10
    weight_scale = 2e-2
    batch_sizes = [5,10,50]
    lr = 10**(-3.5)
    solver_bsize = batch_sizes[0]

    print('No normalization: batch size = ',solver_bsize)
    model = FullyConnectedNet(hidden_dims, weight_scale=weight_scale, normalization=None)
    solver = Solver(model, small_data,
                    num_epochs=n_epochs, batch_size=solver_bsize,
                    update_rule='adam',
                    optim_config={
                      'learning_rate': lr,
                    },
                    verbose=False)
    solver.train()

    bn_solvers = []
    for i in range(len(batch_sizes)):
        b_size=batch_sizes[i]
        print('Normalization: batch size = ',b_size)
        bn_model = FullyConnectedNet(hidden_dims, weight_scale=weight_scale, normalization=normalization_mode)
        bn_solver = Solver(bn_model, small_data,
                        num_epochs=n_epochs, batch_size=b_size,
                        update_rule='adam',
                        optim_config={
                          'learning_rate': lr,
                        },
                        verbose=False)
        bn_solver.train()
        bn_solvers.append(bn_solver)

    return bn_solvers, solver, batch_sizes

batch_sizes = [5,10,50]
bn_solvers_bsize, solver_bsize, batch_sizes = run_batchsize_experiments('batchnorm')
```

    No normalization: batch size =  5
    Normalization: batch size =  5
    Normalization: batch size =  10
    Normalization: batch size =  50

```python
plt.subplot(2, 1, 1)
plot_training_history('Training accuracy (Batch Normalization)','Epoch', solver_bsize, bn_solvers_bsize, \
                      lambda x: x.train_acc_history, bl_marker='-^', bn_marker='-o', labels=batch_sizes)
plt.subplot(2, 1, 2)
plot_training_history('Validation accuracy (Batch Normalization)','Epoch', solver_bsize, bn_solvers_bsize, \
                      lambda x: x.val_acc_history, bl_marker='-^', bn_marker='-o', labels=batch_sizes)

plt.gcf().set_size_inches(15, 10)
plt.show()
```

## Inline Question 2:
Describe the results of this experiment. What does this imply about the relationship between batch normalization and batch size? Why is this relationship observed?

## Answer:
Interestingly the baseline outperforms the batch norm trials.  Regardless, the lower the batch size, the less reliable the mean and variance estimations are (potentially acting as a harsher regularizer) which will slow down training.  This slower learning may indicate that the BN models should be trained for more epochs.  The smaller the batch size, the lower the training accuracy for each epoch.

# Layer Normalization
Batch normalization has proved to be effective in making networks easier to train, but the dependency on batch size makes it less useful in complex networks which have a cap on the input batch size due to hardware limitations.

Several alternatives to batch normalization have been proposed to mitigate this problem; one such technique is Layer Normalization [2]. Instead of normalizing over the batch, we normalize over the features. In other words, when using Layer Normalization, each feature vector corresponding to a single datapoint is normalized based on the sum of all terms within that feature vector.

[2] [Ba, Jimmy Lei, Jamie Ryan Kiros, and Geoffrey E. Hinton. "Layer Normalization." stat 1050 (2016): 21.](https://arxiv.org/pdf/1607.06450.pdf)

## Inline Question 3:
Which of these data preprocessing steps is analogous to batch normalization, and which is analogous to layer normalization?

1. Scaling each image in the dataset, so that the RGB channels for each row of pixels within an image sums up to 1.
2. Scaling each image in the dataset, so that the RGB channels for all pixels within an image sums up to 1.
3. Subtracting the mean image of the dataset from each image in the dataset.
4. Setting all RGB values to either 0 or 1 depending on a given threshold.

## Answer:
1) Layer norm
2) Layer norm
3) Batch norm
4) Layer norm

# Layer Normalization: Implementation

Now you'll implement layer normalization. This step should be relatively straightforward, as conceptually the implementation is almost identical to that of batch normalization. One significant difference though is that for layer normalization, we do not keep track of the moving moments, and the testing phase is identical to the training phase, where the mean and variance are directly calculated per datapoint.

Here's what you need to do:

* In `cs231n/layers.py`, implement the forward pass for layer normalization in the function `layernorm_backward`.

Run the cell below to check your results.
* In `cs231n/layers.py`, implement the backward pass for layer normalization in the function `layernorm_backward`.

Run the second cell below to check your results.
* Modify `cs231n/classifiers/fc_net.py` to add layer normalization to the `FullyConnectedNet`. When the `normalization` flag is set to `"layernorm"` in the constructor, you should insert a layer normalization layer before each ReLU nonlinearity.

Run the third cell below to run the batch size experiment on layer normalization.

```python
# Check the training-time forward pass by checking means and variances
# of features both before and after layer normalization

# Simulate the forward pass for a two-layer network
np.random.seed(231)
N, D1, D2, D3 =4, 50, 60, 3
X = np.random.randn(N, D1)
W1 = np.random.randn(D1, D2)
W2 = np.random.randn(D2, D3)
a = np.maximum(0, X.dot(W1)).dot(W2)

print('Before layer normalization:')
print_mean_std(a,axis=1)

gamma = np.ones(D3)
beta = np.zeros(D3)
# Means should be close to zero and stds close to one
print('After layer normalization (gamma=1, beta=0)')
a_norm, _ = layernorm_forward(a, gamma, beta, {'mode': 'train'})
print_mean_std(a_norm,axis=1)

gamma = np.asarray([3.0,3.0,3.0])
beta = np.asarray([5.0,5.0,5.0])
# Now means should be close to beta and stds close to gamma
print('After layer normalization (gamma=', gamma, ', beta=', beta, ')')
a_norm, _ = layernorm_forward(a, gamma, beta, {'mode': 'train'})
print_mean_std(a_norm,axis=1)
```

    Before layer normalization:
      means:  [-59.06673243 -47.60782686 -43.31137368 -26.40991744]
      stds:   [10.07429373 28.39478981 35.28360729  4.01831507]

    After layer normalization (gamma=1, beta=0)
      means:  [ 4.81096644e-16 -7.40148683e-17  2.22044605e-16 -5.92118946e-16]
      stds:   [0.99999995 0.99999999 1.         0.99999969]

    After layer normalization (gamma= [3. 3. 3.] , beta= [5. 5. 5.] )
      means:  [5. 5. 5. 5.]
      stds:   [2.99999985 2.99999998 2.99999999 2.99999907]

```python
# Gradient check batchnorm backward pass
np.random.seed(231)
N, D = 4, 5
x = 5 * np.random.randn(N, D) + 12
gamma = np.random.randn(D)
beta = np.random.randn(D)
dout = np.random.randn(N, D)

ln_param = {}
fx = lambda x: layernorm_forward(x, gamma, beta, ln_param)[0]
fg = lambda a: layernorm_forward(x, a, beta, ln_param)[0]
fb = lambda b: layernorm_forward(x, gamma, b, ln_param)[0]

dx_num = eval_numerical_gradient_array(fx, x, dout)
da_num = eval_numerical_gradient_array(fg, gamma.copy(), dout)
db_num = eval_numerical_gradient_array(fb, beta.copy(), dout)

_, cache = layernorm_forward(x, gamma, beta, ln_param)
dx, dgamma, dbeta = layernorm_backward(dout, cache)

#You should expect to see relative errors between 1e-12 and 1e-8
print('dx error: ', rel_error(dx_num, dx))
print('dgamma error: ', rel_error(da_num, dgamma))
print('dbeta error: ', rel_error(db_num, dbeta))
```

    dx error:  2.107277492956569e-09
    dgamma error:  4.519489546032799e-12
    dbeta error:  2.5842537629899423e-12

# Layer Normalization and batch size

We will now run the previous batch size experiment with layer normalization instead of batch normalization. Compared to the previous experiment, you should see a markedly smaller influence of batch size on the training history!

```python
ln_solvers_bsize, solver_bsize, batch_sizes = run_batchsize_experiments('layernorm')

plt.subplot(2, 1, 1)
plot_training_history('Training accuracy (Layer Normalization)','Epoch', solver_bsize, ln_solvers_bsize, \
                      lambda x: x.train_acc_history, bl_marker='-^', bn_marker='-o', labels=batch_sizes)
plt.subplot(2, 1, 2)
plot_training_history('Validation accuracy (Layer Normalization)','Epoch', solver_bsize, ln_solvers_bsize, \
                      lambda x: x.val_acc_history, bl_marker='-^', bn_marker='-o', labels=batch_sizes)

plt.gcf().set_size_inches(15, 10)
plt.show()
```

## Inline Question 4:
When is layer normalization likely to not work well, and why?

1. Using it in a very deep network
2. Having a very small dimension of features
3. Having a high regularization term

## Answer:
1, likely to do fine, experimentally layer norm sees use in RNN which can be viewed as deep networks when unrolled through time.
2, similar to the batch size and batch normalization argument.
3, similar to layer regularization in that it's invariant with batch/batch size.

```python

```


text:
# The Index

**Prerequisites**

- [Introduction to pandas](https://datascience.quantecon.org/intro.html)

**Outcomes**

- Understand how the index is used to align data
- Know how to set and reset the index
- Understand how to select subsets of data by slicing on index and columns
- Understand that for DataFrames, the column names also align data

## Outline

- [The Index](#The-Index)
  - [So What is this Index?](#So-What-is-this-Index?)
  - [Setting the Index](#Setting-the-Index)
  - [Re-setting the Index](#Re-setting-the-Index)
  - [Choose the Index Carefully](#Choose-the-Index-Carefully)
  - [Exercises](#Exercises)

```python
# Uncomment following line to install on colab
! pip install qeds
```

    Requirement already satisfied: qeds in c:\users\asus\anaconda3\lib\site-packages (0.6.2)
    Requirement already satisfied: quandl in c:\users\asus\anaconda3\lib\site-packages (from qeds) (3.5.0)
    Requirement already satisfied: pandas-datareader in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.8.1)
    Requirement already satisfied: pandas in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.25.1)
    Requirement already satisfied: statsmodels in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.10.1)
    Requirement already satisfied: plotly in c:\users\asus\anaconda3\lib\site-packages (from qeds) (4.5.4)
    Requirement already satisfied: numpy in c:\users\asus\anaconda3\lib\site-packages (from qeds) (1.16.5)
    Requirement already satisfied: quantecon in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.4.6)
    Requirement already satisfied: openpyxl in c:\users\asus\anaconda3\lib\site-packages (from qeds) (3.0.0)
    Requirement already satisfied: pyarrow in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.16.0)
    Requirement already satisfied: scikit-learn in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.21.3)
    Requirement already satisfied: scipy in c:\users\asus\anaconda3\lib\site-packages (from qeds) (1.3.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: seaborn in c:\users\asus\anaconda3\lib\site-packages (from qeds) (0.9.0)
    Requirement already satisfied: six in c:\users\asus\anaconda3\lib\site-packages (from quandl->qeds) (1.12.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: python-dateutil in c:\users\asus\anaconda3\lib\site-packages (from quandl->qeds) (2.8.0)
    Requirement already satisfied: more-itertools in c:\users\asus\anaconda3\lib\site-packages (from quandl->qeds) (7.2.0)
    Requirement already satisfied: lxml in c:\users\asus\anaconda3\lib\site-packages (from pandas-datareader->qeds) (4.4.1)
    Requirement already satisfied: pytz>=2017.2 in c:\users\asus\anaconda3\lib\site-packages (from pandas->qeds) (2019.3)
    Requirement already satisfied: patsy>=0.4.0 in c:\users\asus\anaconda3\lib\site-packages (from statsmodels->qeds) (0.5.1)
    Requirement already satisfied: retrying>=1.3.3 in c:\users\asus\anaconda3\lib\site-packages (from plotly->qeds) (1.3.3)
    Requirement already satisfied: numba>=0.38 in c:\users\asus\anaconda3\lib\site-packages (from quantecon->qeds) (0.45.1)
    Requirement already satisfied: sympy in c:\users\asus\anaconda3\lib\site-packages (from quantecon->qeds) (1.4)
    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: joblib>=0.11 in c:\users\asus\anaconda3\lib\site-packages (from scikit-learn->qeds) (0.13.2)
    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: idna<2.9,>=2.5 in c:\users\asus\anaconda3\lib\site-packages (from requests->qeds) (2.8)
    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: 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: llvmlite>=0.29.0dev0 in c:\users\asus\anaconda3\lib\site-packages (from numba>=0.38->quantecon->qeds) (0.29.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: setuptools in c:\users\asus\anaconda3\lib\site-packages (from kiwisolver>=1.0.1->matplotlib->qeds) (41.4.0)

```python
import pandas as pd
import numpy as np
```

## So What is this Index?

Every Series or DataFrame has an index.

We told you that the index was the “row labels” for the data.

This is true, but an index in pandas does much more than label the rows.

The purpose of this lecture is to understand the importance of the index.

The [pandas
documentation](https://pandas.pydata.org/pandas-docs/stable/dsintro.html)
says

> Data alignment is intrinsic. The link between labels and data will
not be broken unless done so explicitly by you.

In practice, the index and column names are used to make sure the data is
properly aligned when operating on multiple DataFrames.

This is a somewhat abstract concept that is best understood by
example…

Let’s begin by loading some data on GDP components that we collected from
the World Bank’s World Development Indicators Dataset.

```python
url = "https://datascience.quantecon.org/assets/data/wdi_data.csv"
df = pd.read_csv(url)
df.info()

df.head()
```

    <class 'pandas.core.frame.DataFrame'>
    RangeIndex: 72 entries, 0 to 71
    Data columns (total 7 columns):
    country        72 non-null object
    year           72 non-null int64
    GovExpend      72 non-null float64
    Consumption    72 non-null float64
    Exports        72 non-null float64
    Imports        72 non-null float64
    GDP            72 non-null float64
    dtypes: float64(5), int64(1), object(1)
    memory usage: 4.1+ KB

      country
      year
      GovExpend
      Consumption
      Exports
      Imports
      GDP

      0
      Canada
      2017
      0.372665
      1.095475
      0.582831
      0.600031
      1.868164

      1
      Canada
      2016
      0.364899
      1.058426
      0.576394
      0.575775
      1.814016

      2
      Canada
      2015
      0.358303
      1.035208
      0.568859
      0.575793
      1.794270

      3
      Canada
      2014
      0.353485
      1.011988
      0.550323
      0.572344
      1.782252

      4
      Canada
      2013
      0.351541
      0.986400
      0.518040
      0.558636
      1.732714

We’ll also extract a couple smaller DataFrames we can use in examples.

```python
df_small = df.head(5)
df_small
```

      country
      year
      GovExpend
      Consumption
      Exports
      Imports
      GDP

      0
      Canada
      2017
      0.372665
      1.095475
      0.582831
      0.600031
      1.868164

      1
      Canada
      2016
      0.364899
      1.058426
      0.576394
      0.575775
      1.814016

      2
      Canada
      2015
      0.358303
      1.035208
      0.568859
      0.575793
      1.794270

      3
      Canada
      2014
      0.353485
      1.011988
      0.550323
      0.572344
      1.782252

      4
      Canada
      2013
      0.351541
      0.986400
      0.518040
      0.558636
      1.732714

```python
df_tiny = df.iloc[[0, 3, 2, 4], :]
df_tiny
```

      country
      year
      GovExpend
      Consumption
      Exports
      Imports
      GDP

      0
      Canada
      2017
      0.372665
      1.095475
      0.582831
      0.600031
      1.868164

      3
      Canada
      2014
      0.353485
      1.011988
      0.550323
      0.572344
      1.782252

      2
      Canada
      2015
      0.358303
      1.035208
      0.568859
      0.575793
      1.794270

      4
      Canada
      2013
      0.351541
      0.986400
      0.518040
      0.558636
      1.732714

```python
im_ex = df_small[["Imports", "Exports"]]
im_ex_copy = im_ex.copy()
im_ex_copy
```

      Imports
      Exports

      0
      0.600031
      0.582831

      1
      0.575775
      0.576394

      2
      0.575793
      0.568859

      3
      0.572344
      0.550323

      4
      0.558636
      0.518040

Observe what happens when we evaluate `im_ex + im_ex_copy`.

```python
im_ex + im_ex_copy
```

      Imports
      Exports

      0
      1.200063
      1.165661

      1
      1.151550
      1.152787

      2
      1.151585
      1.137718

      3
      1.144688
      1.100646

      4
      1.117272
      1.036081

Notice that this operated *elementwise*, meaning that the `+`
operation was applied to each element of `im_ex` and the corresponding
element of `im_ex_copy`.

Let’s take a closer look at `df_tiny`:

```python
df_tiny
```

      country
      year
      GovExpend
      Consumption
      Exports
      Imports
      GDP

      0
      Canada
      2017
      0.372665
      1.095475
      0.582831
      0.600031
      1.868164

      3
      Canada
      2014
      0.353485
      1.011988
      0.550323
      0.572344
      1.782252

      2
      Canada
      2015
      0.358303
      1.035208
      0.568859
      0.575793
      1.794270

      4
      Canada
      2013
      0.351541
      0.986400
      0.518040
      0.558636
      1.732714

Relative to `im_ex` notice a few things:

- The row labeled `1` appears in `im_ex` but not `df_tiny`.
- For row labels that appear in both, they are not in the same position
  within each DataFrame.
- Certain columns appear only in `df_tiny`.
- The `Imports` and `Exports` columns are the 6th and 5th columns of
  `df_tiny` and the 1st and 2nd of `im_ex`, respectively.

Now, let’s see what happens when we try `df_tiny + im_ex`.

```python
im_ex_tiny = df_tiny + im_ex
im_ex_tiny
```

      Consumption
      Exports
      GDP
      GovExpend
      Imports
      country
      year

      0
      NaN
      1.165661
      NaN
      NaN
      1.200063
      NaN
      NaN

      1
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN

      2
      NaN
      1.137718
      NaN
      NaN
      1.151585
      NaN
      NaN

      3
      NaN
      1.100646
      NaN
      NaN
      1.144688
      NaN
      NaN

      4
      NaN
      1.036081
      NaN
      NaN
      1.117272
      NaN
      NaN

Whoa, a lot happened! Let’s break it down.

### Automatic Alignment

For all (row, column) combinations that appear in both DataFrames (e.g.
rows `[1, 3]` and columns `[Imports, Exports]`), the value of `im_ex_tiny`
is equal to `df_tiny.loc[row, col] + im_ex.loc[row, col]`.

This happened even though the rows and columns were not in the same
order.

We refer to this as pandas *aligning* the data for us.

To see how awesome this is, think about how to do something similar in
Excel:

- `df_tiny` and `im_ex` would be in different sheets.
- The index and column names would be the first column and row in each
  sheet.
- We would have a third sheet to hold the sum.
- For each label in the first row and column of *either* the `df_tiny`
  sheet or the `im_ex` sheet we would have to do a `IFELSE` to check
  if the label exists in the other sheet and then a `VLOOKUP` to
  extract the value.

In pandas, this happens automatically, behind the scenes, and *very
quickly*.

### Handling Missing Data

For all elements in row `1` or columns
`["country", "year", "GovExpend", "Consumption", "GDP"]`,
the value in `im_ex_tiny` is `NaN`.

This is how pandas represents *missing data*.

So, when pandas was trying to look up the values in `df_tiny` and `im_ex`, it could
only find a value in one DataFrame: the other value was missing.

When pandas tries to add a number to something that is missing, it says
that the result is missing (spelled `NaN`).

## Setting the Index

For a DataFrame `df`, the `df.set_index` method allows us to use one
(or more) of the DataFrame’s columns as the index.

Here’s an example.

```python
# first, create the DataFrame
df_year = df.set_index(["year"])
df_year.head()
```

      country
      GovExpend
      Consumption
      Exports
      Imports
      GDP

      year

      2017
      Canada
      0.372665
      1.095475
      0.582831
      0.600031
      1.868164

      2016
      Canada
      0.364899
      1.058426
      0.576394
      0.575775
      1.814016

      2015
      Canada
      0.358303
      1.035208
      0.568859
      0.575793
      1.794270

      2014
      Canada
      0.353485
      1.011988
      0.550323
      0.572344
      1.782252

      2013
      Canada
      0.351541
      0.986400
      0.518040
      0.558636
      1.732714

Now that the year is on the index, we can use `.loc` to extract all the
data for a specific year.

```python
df_year.loc[2010]
```

      country
      GovExpend
      Consumption
      Exports
      Imports
      GDP

      year

      2010
      Canada
      0.347332
      0.921952
      0.469949
      0.500341
      1.613543

      2010
      Germany
      0.653386
      1.915481
      1.443735
      1.266126
      3.417095

      2010
      United Kingdom
      0.521146
      1.598563
      0.690824
      0.745065
      2.452900

      2010
      United States
      2.510143
      10.185836
      1.846280
      2.360183
      14.992053

This would be helpful, for example, if we wanted to compute the difference
in the average of all our variables from one year to the next.

```python
df_year.loc[2009].mean() - df_year.loc[2008].mean()
```

    GovExpend      0.033317
    Consumption   -0.042998
    Exports       -0.121425
    Imports       -0.140042
    GDP           -0.182610
    dtype: float64

Notice that pandas did a few things for us.

- After computing `.mean()`, the row labels (index) were the former column names.
- These column names were used to align data when we wanted asked pandas to
  compute the difference.

Suppose that someone asked you, “What was the GDP in the US in 2010?”

To compute that using `df_year` you might do something like this:

```python
df_year.loc[df_year["country"] == "United States", "GDP"].loc[2010]
```

    14.992052727

That was a lot of work!

Now, suppose that after seeing you extract that data, your friend asks you
“What about GDP in Germany and the UK in 2010?”

To answer that question, you might write.

```python
df_year.loc[df_year["country"].isin(["United Kingdom", "Germany"]), "GDP"].loc[2010]
```

    year
    2010    3.417095
    2010    2.452900
    Name: GDP, dtype: float64

Notice that this code is similar to the code above, but now provides a result
that is ambiguous.

The two elements in the series both have with label 2010.

How do we know which is which?

We might think that the first value corresponds to the United Kingdom because
that is what we listed first in the call to `isin`, but we would be wrong!

Let’s check.

```python
df_year.loc[2010]
```

      country
      GovExpend
      Consumption
      Exports
      Imports
      GDP

      year

      2010
      Canada
      0.347332
      0.921952
      0.469949
      0.500341
      1.613543

      2010
      Germany
      0.653386
      1.915481
      1.443735
      1.266126
      3.417095

      2010
      United Kingdom
      0.521146
      1.598563
      0.690824
      0.745065
      2.452900

      2010
      United States
      2.510143
      10.185836
      1.846280
      2.360183
      14.992053

Setting just the year as index has one more potential issue: we will
get data alignment only on the year, which may not be sufficient.

To demonstrate this point, suppose now you are asked to use our WDI dataset
to compute an approximation for net exports and investment in in 2009.

As a seasoned economist, you would remember the expenditure formula for GDP is
written

$$
GDP = Consumption + Investment + GovExpend + Net Exports
$$

which we can rearrange to compute investment as a function of the variables in
our DataFrame…

$$
Investment = GDP - Consumption - GovExpend - Net Exports
$$

Note that we can compute NetExports as `Exports - Imports`.

```python
nx = df_year["Exports"] - df_year["Imports"]
nx.head(19)
```

    year
    2017   -0.017201
    2016    0.000619
    2015   -0.006934
    2014   -0.022021
    2013   -0.040596
    2012   -0.041787
    2011   -0.035878
    2010   -0.030393
    2009    0.000896
    2008    0.004068
    2007    0.032451
    2006    0.053530
    2005    0.072729
    2004    0.091902
    2003    0.097794
    2002    0.121850
    2001    0.122673
    2000    0.118702
    2017    0.264214
    dtype: float64

Now, suppose that we accidentally had a bug in our code that swapped
the data for Canada and Germany’s net exports in 2017.

>**Note**
>
>This example is contrived, but if you were getting unclean data from
some resource or doing more complicated operations, this type of mistake
becomes increasingly likely.

```python
ca17 = nx.iloc[[0]]
g17 = nx.iloc[[18]]
nx.iloc[[0]] = g17
nx.iloc[[18]] = ca17

nx.head(19)
```

    year
    2017    0.264214
    2016    0.000619
    2015   -0.006934
    2014   -0.022021
    2013   -0.040596
    2012   -0.041787
    2011   -0.035878
    2010   -0.030393
    2009    0.000896
    2008    0.004068
    2007    0.032451
    2006    0.053530
    2005    0.072729
    2004    0.091902
    2003    0.097794
    2002    0.121850
    2001    0.122673
    2000    0.118702
    2017   -0.017201
    dtype: float64

Notice that if we now add `nx` to the DataFrame and compute investment
pandas doesn’t complain.

```python
df_year["NetExports"] = nx
df_year["Investment"] = df_year.eval("GDP - Consumption - GovExpend - NetExports")
df_year.head(19)
```

      country
      GovExpend
      Consumption
      Exports
      Imports
      GDP
      NetExports
      Investment

      year

      2017
      Canada
      0.372665
      1.095475
      0.582831
      0.600031
      1.868164
      0.264214
      0.135811

      2016
      Canada
      0.364899
      1.058426
      0.576394
      0.575775
      1.814016
      0.000619
      0.390072

      2015
      Canada
      0.358303
      1.035208
      0.568859
      0.575793
      1.794270
      -0.006934
      0.407692

      2014
      Canada
      0.353485
      1.011988
      0.550323
      0.572344
      1.782252
      -0.022021
      0.438800

      2013
      Canada
      0.351541
      0.986400
      0.518040
      0.558636
      1.732714
      -0.040596
      0.435369

      2012
      Canada
      0.354342
      0.961226
      0.505969
      0.547756
      1.693428
      -0.041787
      0.419647

      2011
      Canada
      0.351887
      0.943145
      0.492349
      0.528227
      1.664240
      -0.035878
      0.405086

      2010
      Canada
      0.347332
      0.921952
      0.469949
      0.500341
      1.613543
      -0.030393
      0.374652

      2009
      Canada
      0.339686
      0.890078
      0.440692
      0.439796
      1.565291
      0.000896
      0.334631

      2008
      Canada
      0.330766
      0.889602
      0.506350
      0.502281
      1.612862
      0.004068
      0.388425

      2007
      Canada
      0.318777
      0.864012
      0.530453
      0.498002
      1.596876
      0.032451
      0.381636

      2006
      Canada
      0.311382
      0.827643
      0.524461
      0.470931
      1.564608
      0.053530
      0.372053

      2005
      Canada
      0.303043
      0.794390
      0.519950
      0.447222
      1.524608
      0.072729
      0.354447

      2004
      Canada
      0.299854
      0.764357
      0.508657
      0.416754
      1.477317
      0.091902
      0.321203

      2003
      Canada
      0.294335
      0.741796
      0.481993
      0.384199
      1.433089
      0.097794
      0.299164

      2002
      Canada
      0.286094
      0.721974
      0.490465
      0.368615
      1.407725
      0.121850
      0.277806

      2001
      Canada
      0.279767
      0.694230
      0.484696
      0.362023
      1.366590
      0.122673
      0.269921

      2000
      Canada
      0.270553
      0.677713
      0.499526
      0.380823
      1.342805
      0.118702
      0.275837

      2017
      Germany
      0.745579
      2.112009
      1.930563
      1.666348
      3.883870
      -0.017201
      1.043482

Because we didn’t also have data alignment on the country, we would have overstated Canada’s investment by 281 billion USD and understated Germany’s by the
same amount.

To make these types operation easier, we need to include both the year
and country in the index…

### Setting a Hierarchical Index

Include multiple columns in the index is advantageous in some situations.

These situations might include:

- When we need more than one piece of information (column) to identify an
  observation (as in the Germany and UK GDP example above)
- When we need data-alignment by more than one column

To achieve multiple columns in the index, we pass a list of multiple column
names to `set_index`.

```python
wdi = df.set_index(["country", "year"])
wdi.head(20)
```

      GovExpend
      Consumption
      Exports
      Imports
      GDP

      country
      year

      Canada
      2017
      0.372665
      1.095475
      0.582831
      0.600031
      1.868164

      2016
      0.364899
      1.058426
      0.576394
      0.575775
      1.814016

      2015
      0.358303
      1.035208
      0.568859
      0.575793
      1.794270

      2014
      0.353485
      1.011988
      0.550323
      0.572344
      1.782252

      2013
      0.351541
      0.986400
      0.518040
      0.558636
      1.732714

      2012
      0.354342
      0.961226
      0.505969
      0.547756
      1.693428

      2011
      0.351887
      0.943145
      0.492349
      0.528227
      1.664240

      2010
      0.347332
      0.921952
      0.469949
      0.500341
      1.613543

      2009
      0.339686
      0.890078
      0.440692
      0.439796
      1.565291

      2008
      0.330766
      0.889602
      0.506350
      0.502281
      1.612862

      2007
      0.318777
      0.864012
      0.530453
      0.498002
      1.596876

      2006
      0.311382
      0.827643
      0.524461
      0.470931
      1.564608

      2005
      0.303043
      0.794390
      0.519950
      0.447222
      1.524608

      2004
      0.299854
      0.764357
      0.508657
      0.416754
      1.477317

      2003
      0.294335
      0.741796
      0.481993
      0.384199
      1.433089

      2002
      0.286094
      0.721974
      0.490465
      0.368615
      1.407725

      2001
      0.279767
      0.694230
      0.484696
      0.362023
      1.366590

      2000
      0.270553
      0.677713
      0.499526
      0.380823
      1.342805

      Germany
      2017
      0.745579
      2.112009
      1.930563
      1.666348
      3.883870

      2016
      0.734014
      2.075615
      1.844949
      1.589495
      3.801859

Notice that in the display above, the row labels seem to have two
*levels* now.

The *outer* (or left-most) level is named `country` and the *inner* (or
right-most) level is named `year`.

When a DataFrame’s index has multiple levels, we (and the pandas documentation)
refer to the DataFrame as having a hierarchical index.

### Slicing a Hierarchical Index

Now, we can answer our friend’s questions in a much more straightforward way.

```python
wdi.loc[("United States", 2010), "GDP"]
```

    14.992052727

```python
wdi.loc[(["United Kingdom", "Germany"], 2010), "GDP"]
```

    country         year
    Germany         2010    3.417095
    United Kingdom  2010    2.452900
    Name: GDP, dtype: float64

As shown above, we can use `wdi.loc` to extract different slices of our
national accounts data.

The rules for using `.loc` with a hierarchically-indexed DataFrame are
similar to the ones we’ve learned for standard DataFrames, but they are a bit
more elaborate as we now have more structure to our data.

**Slicing rules**

pandas slicing reacts differently to `list`s and `tuple`s.

It does this to provide more flexibility to select the
data you want.

`list` in row slicing will be an “or” operation, where it chooses rows
based on whether the index value corresponds to any element of the list.

`tuple` in row slicing will be used to denote a single hierarchical
index and must include a value for each level.

**Row slicing examples**

1. `wdi.loc["United States"]`: all rows where the *outer* most index value is
  equal to `United States`
1. `wdi.loc[("United States", 2010)]`: all rows where the *outer-most* index value
  is equal to `"United States` and the second level is equal to `2010`
1. `wdi.loc[["United States", "Canada"]]`: all rows where the *outer-most* index is
  either `"United States"` or `"Canada"`
1. `wdi.loc[(["United States", "Canada"], [2010, 2011]), :]`: all rows where the
  *outer-most* index is either `"United States` or `"Canada"` AND where the
  second level index is either `2010` or `2011`
1. `wdi.loc[[("United States", 2010), ("Canada", 2011)], :]`: all rows where the the
  two hierarchical indices are either `("United States", 2010)` or
  `("Canada", 2011)`

We can also restrict `.loc` to extract certain columns by doing:

1. `wdi.loc[rows, GDP]`: return the rows specified by rows (see rules
  above) and only column named `GDP` (returned object will be a
  Series)
1. `df.loc[rows, ["GDP", "Consumption"]]`: return the rows specified by rows
  (see rules above) and only columns `GDP` and `Consumption`

### Alignment with `MultiIndex`

The data alignment features we talked about above also apply to a
`MultiIndex` DataFrame.

### `pd.IndexSlice`

When we want to extract rows for a few values of the outer index and all
values for an inner index level, we can use the convenient
`df.loc[[id11, id22]]` shorthand.

We can use this notation to extract all the data for the United States and
Canada.

```python
wdi.loc[["United States", "Canada"]]
```

      GovExpend
      Consumption
      Exports
      Imports
      GDP

      country
      year

      Canada
      2017
      0.372665
      1.095475
      0.582831
      0.600031
      1.868164

      2016
      0.364899
      1.058426
      0.576394
      0.575775
      1.814016

      2015
      0.358303
      1.035208
      0.568859
      0.575793
      1.794270

      2014
      0.353485
      1.011988
      0.550323
      0.572344
      1.782252

      2013
      0.351541
      0.986400
      0.518040
      0.558636
      1.732714

      2012
      0.354342
      0.961226
      0.505969
      0.547756
      1.693428

      2011
      0.351887
      0.943145
      0.492349
      0.528227
      1.664240

      2010
      0.347332
      0.921952
      0.469949
      0.500341
      1.613543

      2009
      0.339686
      0.890078
      0.440692
      0.439796
      1.565291

      2008
      0.330766
      0.889602
      0.506350
      0.502281
      1.612862

      2007
      0.318777
      0.864012
      0.530453
      0.498002
      1.596876

      2006
      0.311382
      0.827643
      0.524461
      0.470931
      1.564608

      2005
      0.303043
      0.794390
      0.519950
      0.447222
      1.524608

      2004
      0.299854
      0.764357
      0.508657
      0.416754
      1.477317

      2003
      0.294335
      0.741796
      0.481993
      0.384199
      1.433089

      2002
      0.286094
      0.721974
      0.490465
      0.368615
      1.407725

      2001
      0.279767
      0.694230
      0.484696
      0.362023
      1.366590

      2000
      0.270553
      0.677713
      0.499526
      0.380823
      1.342805

      United States
      2017
      2.405743
      12.019266
      2.287071
      3.069954
      17.348627

      2016
      2.407981
      11.722133
      2.219937
      2.936004
      16.972348

      2015
      2.373130
      11.409800
      2.222228
      2.881337
      16.710459

      2014
      2.334071
      11.000619
      2.209555
      2.732228
      16.242526

      2013
      2.353381
      10.687214
      2.118639
      2.600198
      15.853796

      2012
      2.398873
      10.534042
      2.045509
      2.560677
      15.567038

      2011
      2.434378
      10.378060
      1.978083
      2.493194
      15.224555

      2010
      2.510143
      10.185836
      1.846280
      2.360183
      14.992053

      2009
      2.507390
      10.010687
      1.646432
      2.086299
      14.617299

      2008
      2.407771
      10.137847
      1.797347
      2.400349
      14.997756

      2007
      2.351987
      10.159387
      1.701096
      2.455016
      15.018268

      2006
      2.314957
      9.938503
      1.564920
      2.395189
      14.741688

      2005
      2.287022
      9.643098
      1.431205
      2.246246
      14.332500

      2004
      2.267999
      9.311431
      1.335978
      2.108585
      13.846058

      2003
      2.233519
      8.974708
      1.218199
      1.892825
      13.339312

      2002
      2.193188
      8.698306
      1.192180
      1.804105
      12.968263

      2001
      2.112038
      8.480461
      1.213253
      1.740797
      12.746262

      2000
      2.040500
      8.272097
      1.287739
      1.790995
      12.620268

However, suppose we wanted to extract the data for all countries, but only the
years 2005, 2007, and 2009.

We cannot do this using `wdi.loc` because the year is on the second level,
not outer-most level of our index.

To get around this limitation, we can use the `pd.IndexSlice` helper.

Here’s an example.

```python
wdi.loc[pd.IndexSlice[:, [2005, 2007, 2009]], :]
```

      GovExpend
      Consumption
      Exports
      Imports
      GDP

      country
      year

      Canada
      2009
      0.339686
      0.890078
      0.440692
      0.439796
      1.565291

      2007
      0.318777
      0.864012
      0.530453
      0.498002
      1.596876

      2005
      0.303043
      0.794390
      0.519950
      0.447222
      1.524608

      Germany
      2009
      0.645023
      1.908393
      1.260525
      1.121914
      3.283144

      2007
      0.605624
      1.894219
      1.442436
      1.213835
      3.441356

      2005
      0.591184
      1.866253
      1.175200
      1.028094
      3.213777

      United Kingdom
      2009
      0.519716
      1.587152
      0.653830
      0.689011
      2.411632

      2007
      0.504549
      1.644789
      0.710200
      0.767699
      2.527327

      2005
      0.490806
      1.578914
      0.640088
      0.715951
      2.403352

      United States
      2009
      2.507390
      10.010687
      1.646432
      2.086299
      14.617299

      2007
      2.351987
      10.159387
      1.701096
      2.455016
      15.018268

      2005
      2.287022
      9.643098
      1.431205
      2.246246
      14.332500

Notice that the `:` in the first part of `[:, ["A", "D"]]`
instructed pandas to give us rows for all values of the outer most index
level and that the `:` just before `]` said grab all the columns.

### Multi-index Columns

The functionality of `MultiIndex` also applies to the column names.

Let’s see how it works.

```python
wdiT = wdi.T  # .T means "transpose" or "swap rows and columns"
wdiT
```

      country
      Canada
      ...
      United States

      year
      2017
      2016
      2015
      2014
      2013
      2012
      2011
      2010
      2009
      2008
      ...
      2009
      2008
      2007
      2006
      2005
      2004
      2003
      2002
      2001
      2000

      GovExpend
      0.372665
      0.364899
      0.358303
      0.353485
      0.351541
      0.354342
      0.351887
      0.347332
      0.339686
      0.330766
      ...
      2.507390
      2.407771
      2.351987
      2.314957
      2.287022
      2.267999
      2.233519
      2.193188
      2.112038
      2.040500

      Consumption
      1.095475
      1.058426
      1.035208
      1.011988
      0.986400
      0.961226
      0.943145
      0.921952
      0.890078
      0.889602
      ...
      10.010687
      10.137847
      10.159387
      9.938503
      9.643098
      9.311431
      8.974708
      8.698306
      8.480461
      8.272097

      Exports
      0.582831
      0.576394
      0.568859
      0.550323
      0.518040
      0.505969
      0.492349
      0.469949
      0.440692
      0.506350
      ...
      1.646432
      1.797347
      1.701096
      1.564920
      1.431205
      1.335978
      1.218199
      1.192180
      1.213253
      1.287739

      Imports
      0.600031
      0.575775
      0.575793
      0.572344
      0.558636
      0.547756
      0.528227
      0.500341
      0.439796
      0.502281
      ...
      2.086299
      2.400349
      2.455016
      2.395189
      2.246246
      2.108585
      1.892825
      1.804105
      1.740797
      1.790995

      GDP
      1.868164
      1.814016
      1.794270
      1.782252
      1.732714
      1.693428
      1.664240
      1.613543
      1.565291
      1.612862
      ...
      14.617299
      14.997756
      15.018268
      14.741688
      14.332500
      13.846058
      13.339312
      12.968263
      12.746262
      12.620268

5 rows × 72 columns

Notice that `wdiT` seems to have two levels of names for the columns.

The same logic laid out in the above row slicing rules applies when we
have a hierarchical index for column names.

```python
wdiT.loc[:, "United States"]
```

      year
      2017
      2016
      2015
      2014
      2013
      2012
      2011
      2010
      2009
      2008
      2007
      2006
      2005
      2004
      2003
      2002
      2001
      2000

      GovExpend
      2.405743
      2.407981
      2.373130
      2.334071
      2.353381
      2.398873
      2.434378
      2.510143
      2.507390
      2.407771
      2.351987
      2.314957
      2.287022
      2.267999
      2.233519
      2.193188
      2.112038
      2.040500

      Consumption
      12.019266
      11.722133
      11.409800
      11.000619
      10.687214
      10.534042
      10.378060
      10.185836
      10.010687
      10.137847
      10.159387
      9.938503
      9.643098
      9.311431
      8.974708
      8.698306
      8.480461
      8.272097

      Exports
      2.287071
      2.219937
      2.222228
      2.209555
      2.118639
      2.045509
      1.978083
      1.846280
      1.646432
      1.797347
      1.701096
      1.564920
      1.431205
      1.335978
      1.218199
      1.192180
      1.213253
      1.287739

      Imports
      3.069954
      2.936004
      2.881337
      2.732228
      2.600198
      2.560677
      2.493194
      2.360183
      2.086299
      2.400349
      2.455016
      2.395189
      2.246246
      2.108585
      1.892825
      1.804105
      1.740797
      1.790995

      GDP
      17.348627
      16.972348
      16.710459
      16.242526
      15.853796
      15.567038
      15.224555
      14.992053
      14.617299
      14.997756
      15.018268
      14.741688
      14.332500
      13.846058
      13.339312
      12.968263
      12.746262
      12.620268

```python
wdiT.loc[:, ["United States", "Canada"]]
```

      country
      Canada
      ...
      United States

      year
      2017
      2016
      2015
      2014
      2013
      2012
      2011
      2010
      2009
      2008
      ...
      2009
      2008
      2007
      2006
      2005
      2004
      2003
      2002
      2001
      2000

      GovExpend
      0.372665
      0.364899
      0.358303
      0.353485
      0.351541
      0.354342
      0.351887
      0.347332
      0.339686
      0.330766
      ...
      2.507390
      2.407771
      2.351987
      2.314957
      2.287022
      2.267999
      2.233519
      2.193188
      2.112038
      2.040500

      Consumption
      1.095475
      1.058426
      1.035208
      1.011988
      0.986400
      0.961226
      0.943145
      0.921952
      0.890078
      0.889602
      ...
      10.010687
      10.137847
      10.159387
      9.938503
      9.643098
      9.311431
      8.974708
      8.698306
      8.480461
      8.272097

      Exports
      0.582831
      0.576394
      0.568859
      0.550323
      0.518040
      0.505969
      0.492349
      0.469949
      0.440692
      0.506350
      ...
      1.646432
      1.797347
      1.701096
      1.564920
      1.431205
      1.335978
      1.218199
      1.192180
      1.213253
      1.287739

      Imports
      0.600031
      0.575775
      0.575793
      0.572344
      0.558636
      0.547756
      0.528227
      0.500341
      0.439796
      0.502281
      ...
      2.086299
      2.400349
      2.455016
      2.395189
      2.246246
      2.108585
      1.892825
      1.804105
      1.740797
      1.790995

      GDP
      1.868164
      1.814016
      1.794270
      1.782252
      1.732714
      1.693428
      1.664240
      1.613543
      1.565291
      1.612862
      ...
      14.617299
      14.997756
      15.018268
      14.741688
      14.332500
      13.846058
      13.339312
      12.968263
      12.746262
      12.620268

5 rows × 36 columns

```python
wdiT.loc[:, (["United States", "Canada"], 2010)]
```

      country
      Canada
      United States

      year
      2010
      2010

      GovExpend
      0.347332
      2.510143

      Consumption
      0.921952
      10.185836

      Exports
      0.469949
      1.846280

      Imports
      0.500341
      2.360183

      GDP
      1.613543
      14.992053

## Re-setting the Index

The `df.reset_index` method will move one or more level of the index
back into the DataFrame as a normal column.

With no additional arguments, it moves all levels out of the index and
sets the index of the returned DataFrame to the default of
`range(df.shape[0])`.

```python
wdi.reset_index()
```

      country
      year
      GovExpend
      Consumption
      Exports
      Imports
      GDP

      0
      Canada
      2017
      0.372665
      1.095475
      0.582831
      0.600031
      1.868164

      1
      Canada
      2016
      0.364899
      1.058426
      0.576394
      0.575775
      1.814016

      2
      Canada
      2015
      0.358303
      1.035208
      0.568859
      0.575793
      1.794270

      3
      Canada
      2014
      0.353485
      1.011988
      0.550323
      0.572344
      1.782252

      4
      Canada
      2013
      0.351541
      0.986400
      0.518040
      0.558636
      1.732714

      ...
      ...
      ...
      ...
      ...
      ...
      ...
      ...

      67
      United States
      2004
      2.267999
      9.311431
      1.335978
      2.108585
      13.846058

      68
      United States
      2003
      2.233519
      8.974708
      1.218199
      1.892825
      13.339312

      69
      United States
      2002
      2.193188
      8.698306
      1.192180
      1.804105
      12.968263

      70
      United States
      2001
      2.112038
      8.480461
      1.213253
      1.740797
      12.746262

      71
      United States
      2000
      2.040500
      8.272097
      1.287739
      1.790995
      12.620268

72 rows × 7 columns

## Choose the Index Carefully

So, now that we know that we use index and column names for
aligning data, “how should we pick the index?” is a natural question to ask.

To guide us to the right answer, we will list the first two components
to [Hadley Wickham’s](http://hadley.nz/) description of [tidy
data](http://vita.had.co.nz/papers/tidy-data.html):

1. Each column should each have one variable.
1. Each row should each have one observation.

If we strive to have our data in a tidy form (we should), then when
choosing the index, we should set:

- the row labels (index) to be a unique identifier for an observation
  of data
- the column names to identify one variable

For example, suppose we are looking data on interest rates.

Each column might represent one bond or asset and each row might
represent the date.

Using hierarchical row and column indices allows us to store higher
dimensional data in our (inherently) two dimensional DataFrame.

### Know Your Goal

The correct column(s) to choose for the index often depends on the context of
your analysis.

For example, if I were studying how GDP and consumption evolved over time for
various countries, I would want time (year) and country name on the index

On the other hand, if I were trying to look at the differences across countries
and variables within a particular year, I may opt to put the country and
variable on the index and have years be columns.

Following the tidy data rules above and thinking about how you intend to *use*
the data – and a little practice – will enable you to consistently select the
correct index.

ext: ipynb

max_stars_repo_name: remi-sudo/Classes

max_issues_repo_name: remi-sudo/Classes

max_forks_repo_name: remi-sudo/Classes


text:
# Big Signal Analysis of Reaction Networks

This notebook explores how to analyze the qualitative characteristics of a reaction network, such as:
* Number of fixed points
* Types of fixed points
* Bifurcations

The core challenge is to find the number of fixed points and explore the conditions under which their characteristics change.
Finding the number of fixed points requires solving a quadratic system.
Characterizing the fixed points requires calculating characteristic equations.
The technical approach here is primarily using symbolic algebra.

# Preliminaries

## Imports

```python
import tellurium as te
import sympy
import matplotlib.pyplot as plt
import numpy as np
from common_python.sympy import sympyUtil as su
from common_python.ODEModel.ODEModel import ODEModel
import matplotlib.pyplot as plt
```

## Constants

```python
su.addSymbols("S0 S1 S2 S3 S4 S5 S6 S7 S8 S9 k0 k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 k17 k18 k19 k20")
su.addSymbols("J0 J1 J2 J3 J4 J5 J6 J7 J8 J9 J10 J11 J12 J13 J14 J15 J16 J`7 J18 J19 J20")
FLUXES = [J0, J1, J2, J3, J4, J5, J6, J7, J8, J9, J10, J11, J12, J13, J14, J15]
SPECIES = [S0, S1, S2, S3, S4, S5, S6, S7, S8]
```

```python
MODEL = """
J0: S0 -> S2; k0*S0
J1: S3 + S3 -> S0; k5*S3*S3
J2: S2 -> S3; k6*S2
J3: S3 ->; k9*S2*S3
J4: -> S3; k10*S0
k0 = 6+0.9011095014634776
k5 = 1.4823891153952284
k6 = -10+15.149868787476994
k9 = 91.19197034598812
k10 = 200
S0 = 1.0
S1 = 5.0
S2 = 9.0
S3 = 3.0
S4 = 10.0
"""
RR = te.loada(MODEL)
```

# Helper Functions

```python
# Create dictionary relating reaction names to species
def mkStoichiometryExpressionDct(rr):
    stoichiometryArr = rr.getFullStoichiometryMatrix()
    reactionNames = [n[1:] if n[0] == "_" else n for n in stoichiometryArr.colnames]
    stoichiometryArr.colnames = reactionNames
    speciesNames = stoichiometryArr.rownames
    dct = {}
    for idx, species in enumerate(speciesNames):
        sym = eval(species)
        entry = ["%d*%s" % (stoichiometryArr[idx, n], reactionNames[n]) for n in range(len(reactionNames))]
        expressionStr = " + ".join(entry)
        dct[sym] = eval(expressionStr)
        dct[sym] = sympy.simplify(dct[sym])
    return dct

# Tests
stoichiometryDct = mkStoichiometryExpressionDct(RR)
stoichiometryDct
```

    {S0: -J0 + J1, S2: J0 - J2, S3: -2*J1 + J2 - J3 + J4}

```python
# Do plot
def plotEigenInfo(rr, title="", k10Range=None, fixedPointIdx=1, **kwargs):
    """
    Plots information about the dominant eigenvalue and fixed point for the
    indicated fixed point.

    Parameters
    ----------
    rr: ExtendedRoadrunner
    title: str
    k10Range: range of k10
    fixedPointIdx: int
        index of the fixed point to study
    kwargs: dict
        fractional adjustment in value
    """
    if k10Range is None:
        k10Range = [150 + v for v in range(60)]
    # Construct the data
    subDct = {k0: rr.k0, k5: rr.k5, k6: rr.k6, k9: rr.k9, k10: rr.k10}
    for key, value in kwargs.items():
        if isinstance(key, str):
            newKey = [s for s in subDct.keys() if s.name == key][0]
        else:
            newKey = key
        subDct[newKey] = subDct[newKey] * value
    xvs = []
    reals = []
    imags = []
    fixedPointDcts = []
    for c10 in k10Range:
        subDct[k10] = c10
        dominantReal = None
        dominantImag = None
        model = ODEModel(REDUCED_STATE_DCT, subs=subDct, isEigenvecs=False)
        # Find the dominant eigenvalue for the fixed points
        if len(model.fixedPoints) > fixedPointIdx:
            fixedPointDcts.append(dict(model.fixedPoints[fixedPointIdx].valueDct))
            # Find the dominant eigenvalue
            for entry in model.fixedPoints[fixedPointIdx].eigenEntries:
                value = entry.value
                if isinstance(value, complex):
                    real, imag = su.asRealImag(value)
                else:
                    real = value
                    imag = 0
                if (dominantReal is None) or (real > dominantReal):
                    dominantReal = real
                    dominantImag = np.abs(imag)
        xvs.append(c10)
        reals.append(dominantReal)
        imags.append(dominantImag)
    # Plot the dominant eigenvalue
    _, ax = plt.subplots(1)
    ax.plot(xvs, reals, color="blue")
    ax.plot(xvs, imags, color="brown")
    ax.plot([xvs[0], xvs[-1]], [0, 0], linestyle="--", color="black")
    ax.legend(["real", "imag"])
    ax.set_title(title)
    ax.set_xlabel("k10")
    # Plot the indexed fixed point
    states = list(fixedPointDcts[0].keys())
    _, ax = plt.subplots(1)
    COLORS = ["red", "green", "brown"]
    for idx, state in enumerate(states):
        yvs = [f[state] for f in fixedPointDcts]
        ax.plot(xvs, yvs, color=COLORS[idx])
    ax.legend(states)
    ax.set_title("Fixed Points")
    ax.set_xlabel("k10")
    return fixedPointDcts

# Test
dcts = plotEigenInfo(RR, k10Range=[100 + 5*v for v in range(5)], k9=1, title="Dominant eigenvalue for 2nd fixed point.")
```

```python
def runSim(model=MODEL, endTime=100, startTime=0, **kwargs):
    def findIdx(arr, time):
        """Finds the index of the time in the simulation results array."""
        bestIdx = 0
        diff = np.abs(arr[0, 0] - time)
        for idx, value in enumerate(arr[:, 0]):
            if np.abs(value - time) < diff:
                diff = np.abs(value - time)
                bestIdx = idx
        return bestIdx
    rr = te.loada(MODEL)
    # Adjust the parameters
    for key, value in kwargs.items():
        if isinstance(key, sympy.core.symbol.Symbol):
            newKey = s.name
        else:
            newKey = key
        rr[newKey] = rr[newKey] * value
    #rr.plot(rr.simulate(startTime, endTime, 10*endTime))
    arr = rr.simulate(0, endTime, 10*endTime)
    _, ax = plt.subplots(1)
    startIdx = findIdx(arr, startTime)
    endIdx = findIdx(arr, endTime)
    for idx in range(len(arr.colnames[1:])):
        ax.plot(arr[startIdx:endIdx,0], arr[startIdx:endIdx, idx+1])
    ax.legend(arr.colnames[1:])

# Tests
runSim(k9=0.5, startTime=900, endTime=910)
```

# Damped Model 2

```python
MODEL2 = """
var S0
var S1
var S2
var S3
var S4
var S5
var S6
var S7
ext S8
J0: S4 -> S7+S5; k0*S4
J1: S2 -> S4+S4; k1*S2
J2: S4 -> S3+S3; k2*S4
J3: S4 -> S2+S3; k3*S4
J4: S0 -> S5; k4*S0
J5: S5 + S4 -> S5; k5*S5*S4
J6: S5 -> S3; k6*S5
J7: S8 + S3 -> S0; k7*S8*S3
J8: S3 -> S6+S5; k8*S3
J9: S6 + S5 -> S4; k9*S6*S5
J10: S7 + S5 -> S0 + S2; k10*S7*S5
J11: S3 -> S5+S6; k11*S3
J12: S6 + S1 -> S5; k12*S6*S1
J13: S5 -> S5; k13*S5
J14: S1 + S7 -> S1 + S1; k14*S1*S7
k0 = 2.5920480618068815
k1 = 422.2728070204247
k2 = 28.978192374985912
k3 = 29.723263589242986
k4 = 21.04114996098882
k5 = 1.5111236529181926
k6 = 14.363185343334044
k7 = 0.8231126169112812
k8 = 54.27226867691914
k9 = 58.17954213283633
k10 = 10.682986014127339
k11 = 194.08273474192015
k12 = 15.989508525207631
k13 = 13.186614071108659
k14 = 35.67582901156382
S0 = 1.0
S1 = 5.0
S2 = 9.0
S3 = 3.0
S4 = 10.0
S5 = 3.0
S6 = 7.0
S7 = 1.0
S8 = 6.0
"""

rr = te.loada(MODEL2)
rr.plot(rr.simulate())
```

```python
mat = sympy.Matrix(rr.getFullStoichiometryMatrix())
mat
```

```python
SPECIES_FLUX_DCT = mkStoichiometryExpressionDct(rr)
SPECIES_FLUX_DCT
```

```python
nullspace = mat.nullspace()
```

```python
# Kinetics dictionary
kineticDct = {
    J0:  k0*S4,
    J1:  k1*S2,
    J2:  k2*S4,
    J3:  k3*S4,
    J4:  k4*S0,
    J5:  k5*S5*S4,
    J6:  k6*S5,
    J7:  k7*S8*S3,
    J8:  k8*S3,
    J9:  k9*S6*S5,
    J10:  k10*S7*S5,
    J11:  k11*S3,
    J12:  k12*S6*S1,
    J13:  k13*S5,
    J14:  k14*S1*S7,
}
```

```python
STATE_DCT = {s: SPECIES_FLUX_DCT[s].subs(kineticDct) for s in SPECIES_FLUX_DCT.keys() }
STATE_DCT
```

```python
MODEL = ODEModel(STATE_DCT, isFixedPoints=False)
```

```python
# Need to find a linear combination of values in the null space
# such that the kinetic equations hold.
# Have N reactions, M species. So, N - M constants to find.
su.addSymbols("c c_0 c_1 c_2 c_3 c_4 c_5 c_6")
c = sympy.Matrix([c_0, c_1, c_2, c_3, c_4, c_5, c_6])
mat = sympy.Matrix(nullspace).reshape(15, 7)
mat * c
```

```python
105/7
```

```python
# Solve for log(S*)
exprs = [ j - b for j, b in zip(kineticDct.values(), mat*c) ]
exprs = [e.subs({S2: 0, S3: 0}) for e in exprs]
sympy.solve(exprs, [ S5])
```

**Approach**
1. $N$ =  stoichiometry matrix
1. $M$ = nullspace of $N$
1. Substitute 0 for any state variable that must be zero for all vectors $M \star c$.
1. Solve for log of $x_n$ (state variable) in terms of log of $J_n$ (fluxes)
1. We know that the fluxes for the fixed points must be in $M \star c$, where $c$ is a vector.
1. Substitute $J_n$ value from previous (5) into (4) to give an expression for $x_n$ in terms of $c_n$.

Issue: How do I find the $c_i$?

**Approach 2**
1. Solve for $x_n$ in terms of $J$

```python
exprs[1]
```

```python
list(kineticDct.values())[1]
```

```python
sympy.shape(mat)
```

```python
exprs
```

```python
sympy.solve(exprs, [S0, S1, S2, S3, S4, S5, S6, S7, S8])
```

```python
# Looks like I can manually solve for most species
SPECIES_FLUX_DCT = {
    S0: J4/k4,
    S1: (J12/k12) /((J9 / k9) / (J6 / k6)),
    S2: J1 / k1,
    S3: J8 / k8,
    S4: J0 / k0,
    S5: J6 / k6,
    S6: (J9 / k9) / (J6 / k6),
    S7: (J10 / k10) / (J6 / k6),
    S8: (J7 / k7) / (J8 / k8),

}
```

```python
dstateDct =  {s: SPECIES_FLUX_DCT[s].subs(kineticDct) for s in SPECIES_FLUX_DCT.keys()}
```

```python
dstateDct
```

```python
solnDct = sympy.solve(list(SPECIES_FLUX_DICT.values()), list(kineticDct.keys()))
solnDct
```

```python
exprs = [solnDct[j].subs(kineticDct) - kineticDct[j] for j in solnDct.keys()]
exprs
```

```python
#sympy.solve(exprs, list(dstateDct.keys()))
```

# Reduced Model

```python
su.addSymbols("S0 S1 S2 S3 S4 S5 S6")
su.addSymbols("k0 k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 k17 k18 k19 k20")
su.addSymbols("J0 J1 J2 J3 J4 J5 J6 J7 J8 J9 J10 J11 J12 J13 J14 J15 J16 J`7 J18 J19 J20")
REDUCED_FLUXES = [J0, J1, J2, J3, J4]
REDUCED_SPECIES = [S0, S1, S2, S3, S4]
```

```python
MODEL = """
J0: S0 -> S2; k0*S0
J1: S3 + S3 -> S0; k5*S3*S3
J2: S2 -> S3; k6*S2
J3: S3 ->; k9*S2*S3
J4: -> S3; k10*S0
k0 = 6+0.9011095014634776
k5 = 1.4823891153952284
k6 = -10+15.149868787476994
k9 = 91.19197034598812
k10 = 200
S0 = 1.0
S1 = 5.0
S2 = 9.0
S3 = 3.0
S4 = 10.0
"""
```

```python
MODEL = """
J0: S0 -> S2; k0*S0
J1: S3 + S3 -> S0; k5*S3*S3
J2: S2 -> S3; k6*S2
J3: S3 ->; k9*S2*S3
J4: -> S3; k10*S0
k0 = (6+0.9011095014634776)
k5 = 1.4823891153952284
k6 = (-10+15.149868787476994)
k9 = 91.19197034598812  # At k9 * 0.5, use k10 = k155
k10 = 200 # 164, 165 ins a transition from damped to stable oscillations
S0 = 1.0
S2 = 9.0
S3 = 3.0
"""

rr = te.loada(MODEL)
rr.plot(rr.simulate(0, 100, 100))
```

```python
runSim(k9=0.5, k10=140, endTime=1000)
```

```python
runSim(k9=0.5, k10=150, endTime=1000)
```

```python
plot1(k0=0.01, k5=0.1, k6=0.1, k9=0.1, base=100)
```

```python
REDUCED_SPECIES_FLUX_DCT = mkStoichiometryExpressionDct(rr)
REDUCED_SPECIES_FLUX_DCT
```

```python
kineticDct = {
    J0: k0 * S0,
    J1: k5 * S3 * S3,
    J2: k6 * S2,
    J3: k9 * S2 * S3,
    J4: k10*S0,  # Is this really mass action?

}
```

```python
# State equation is wrong for S2. Should be - S2*k6
REDUCED_STATE_DCT = {s: REDUCED_SPECIES_FLUX_DCT[s].subs(kineticDct) for s in REDUCED_SPECIES_FLUX_DCT.keys()}
REDUCED_STATE_DCT
```

```python
sympy.solve(list(REDUCED_STATE_DCT.values()), list(REDUCED_STATE_DCT.keys()))
```

```python
reducedModel = ODEModel(REDUCED_STATE_DCT)
```

```python
# Fixed points
[f.valueDct for f in reducedModel.fixedPoints]
```

```python
# Verify that these are fixed points
for fp in reducedModel.fixedPoints:
    print([sympy.simplify(e.subs(fp.valueDct)) for e in REDUCED_STATE_DCT.values()])
```

```python
# Look at the eigenvectors
if False:
    for entry in reducedModel.fixedPoints[1].eigenEntries:
        for vector in entry.vectors:
            print(vector)
```

```python
eigenvalues = [e.value for e in reducedModel.fixedPoints[1].eigenEntries]
```

**Approach**
1. Find the fixed points.
1. For non-zero fixed points:
    1. Find the eigenvalues in terms of each constant in turn, setting the other constants to 1.
    1. Search for values of constants that result in a positive but near zero real value and significant non-zero imaginary part

**Issue**
1. Eigenvalues have no relationship to the system behavior

## Finding Parameter Values

Given an ODEModel, find values of parameters that result in oscillations at different frequencies.

```python
c0 = rr.k0
c5 = rr.k5
c6 = rr.k6
c9 = rr.k9
c10 = rr.k10
for c10 in [150 + n for n in range(50)]:
    subDct = {k0: c0, k5: c5, k6: c6, k9: c9, k10: c10}
    model = ODEModel(REDUCED_STATE_DCT, subs=subDct)
    entries = model.fixedPoints[1].eigenEntries
    print((c10, [e.value for e in entries]))
```

## Plots

```python
dcts = plotEigenInfo(RR, k10Range=[100 + 5*v for v in range(25)], k9=1, title="Dominant eigenvalue for 2nd fixed point.")
```

```python
runSim(k10=100/200, startTime=0, endTime=10)
```

```python
runSim(k10=160/200, startTime=0, endTime=10)
```

```python
runSim(k9=1, k10=160/200, startTime=990, endTime=1000)
```

```python
runSim(k9=1, k10=170/200, startTime=990, endTime=1000)
```

```python
runSim(k9=1, k10=200/200, startTime=990, endTime=1000)
```

```python
4.5 / (2*np.pi)
```

```python
def plot1(base=150, **kwargs):
    k10Range=[base + 10*v for v in range(10)]
    title = ""
    for key, value in kwargs.items():
        title += " %s: %3.2f " % (str(key), value)
    plotComplexEigenvalue(rr, k10Range=k10Range, title=title, **kwargs)

plot1(k0=1, k5=1, k6=1, k9=1)
```

```python
plot1(k0=0.01, k5=0.1, k6=0.1, k9=0.1, base=100)
```

```python
plot1(k0=1, k5=1, k6=1, k9=0.5, base=100)

```

```python
runSim(k0=1, k5=1, k6=1, k9=0.5, k10=100, endTime=1000)
```

Am i excluding the dominant eigenvalue? Do the plots for all eigenvalues.

max_stars_repo_path: notebooks/Big Signal Analysis of Reaction Networks.ipynb

max_stars_repo_name: sys-bio/IntegralControl

max_issues_repo_name: sys-bio/IntegralControl

max_forks_repo_name: sys-bio/IntegralControl


text:
# Manufactured solution for BP3-like problem

The manufactured solution is going to exactly solve the following problem:

$$
\begin{aligned}
    -\frac{\partial\sigma_{ij}(\boldsymbol{u})}{\partial x_j} &= 0 & \text{ in } & \Omega \\
    \sigma_{ij} &= c_{ijkl} \frac{\partial u_k}{\partial x_l} & \text{ in } & \Omega \\
    u_i &= u_i^* & \text{ on } & \Gamma_D \\
    u_1^- - u_1^+ &= 0 & \text{ on } & \Gamma_F \\
    u_2^- - u_2^+ &= S_2 & \text{ on } & \Gamma_F \\
    u_3^- - u_3^+ &= -S_1 & \text{ on } & \Gamma_F \\
    -(\tau_i^0 + \tau_i) &= (\sigma_n^0-\sigma_n)f(|V|,\psi)\frac{V_i}{|V|} + \eta V_i & \text{ on } & \Gamma_F \\
    \frac{d\psi}{dt} &= g(|V|,\psi) + s(\vec x, t) & \text{ on } & \Gamma_F \\
    \frac{dS}{dt} &= V & \text{ on } & \Gamma_F
\end{aligned}
$$
where
$$
\begin{aligned}
    f(|V|, \psi) &= a \cdot \mathrm{asinh}\left(\frac{|V|}{2V_0} \exp\left(\frac{\psi}{a}\right)\right) \\
    g(|V|, \psi) &= \frac{b V_0}{L} \left(\exp\left(\frac{f_0-\psi}{b}\right) - \frac{|V|}{V_0}\right) \\
    s(\vec{x}, t) &= - g(|V^*|, \psi^*) + \frac{d\psi^*}{dt} \\
    \tau_1 &= -\sigma_{13} \\
    \tau_2 &= \sigma_{12} \\
    \sigma_n &= \sigma_{11}
\end{aligned}
$$

```python
from sympy import *
init_printing()
```

We first define symbols and general functions, which are going to be specified later.

```python
x, y, z, t = symbols('x y z t')
te, tw = symbols('t_e t_w')
V0, Vp = symbols('V_0 V_p')
a, sn, eta = symbols('a sigma_n eta')
K = Function('K')(t)
V1_star = Function('V_1^*')(z, t)
tau1_star = Function('tau_1^*')(z, t)
tau1_pre = Function('tau_1^0')(z)
sn_star = symbols('sigma_n^*')
sn_pre = Function('sigma_n^0')(z)
E, nu = symbols('E nu', real=True, positive=True)
lamda = E*nu/((1-2*nu)*(1+nu))
mu = E/(2*(1+nu))
```

$\psi^*$ is the state variable.

```python
psi_star = a * log((2*V0/V1_star) * sinh((tau1_pre + tau1_star + eta * V1_star) / (a*(sn_star - sn_pre))))
psi_star
```

We need the time derivative $\frac{d\psi^*}{dt}$

```python
psi_star.diff(t).simplify()
```

$u_i^*$ is solution for the displacement.

```python
k = symbols('k', real=True, positive=True)
f = K*x*exp(-k*x)
phi = cos(k*z)*f
phi
```

```python
sxx = phi.diff(z, 2)
szz = phi.diff(x, 2)
sxz = -phi.diff(x).diff(z)
exx = (sxx - nu/(1-nu)*szz)*(1-nu**2)/E
ezz = (szz - nu/(1-nu)*sxx)*(1-nu**2)/E
exz = sxz*(1+nu)/E
ux_star = integrate(exx, x)
uz_star = integrate(ezz, z)
ux_star.simplify()
```

```python
uz_star.simplify()
```

```python
simplify(ux_star.diff(z) + uz_star.diff(x) - 2*exz)
```

Let's introduce slip:

```python
uxm_star = ux_star.subs(k, -k)
uzm_star = uz_star.subs(k, -k)
uxp_star = ux_star
uzp_star = uz_star
```

```python
u = Array([uxm_star, 0, uzm_star])
I = Array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
Ju = permutedims(derive_by_array(u, [x, y, z]), (1, 0))
strain = 0.5 * simplify(Ju + Ju.transpose())
stress1 = simplify(lamda * tensorcontraction(strain, (0, 1)) * I + 2.0 * mu * strain)
eq1 = simplify(stress1[0,0].diff(x) + stress1[0,1].diff(y) +  + stress1[0,2].diff(z))
eq2 = simplify(stress1[1,0].diff(x) + stress1[1,1].diff(y) +  + stress1[1,2].diff(z))
eq3 = simplify(stress1[2,0].diff(x) + stress1[2,1].diff(y) +  + stress1[2,2].diff(z))
(eq1, eq2, eq3)
```

```python
u = Array([uxp_star, 0, uzp_star])
I = Array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
Ju = permutedims(derive_by_array(u, [x, y, z]), (1, 0))
strain = 0.5 * simplify(Ju + Ju.transpose())
stress2 = simplify(lamda * tensorcontraction(strain, (0, 1)) * I + 2.0 * mu * strain)
```

```python
simplify((stress2-stress1).subs(x, 0))
```

$\displaystyle \left[\begin{matrix}0 & 0 & 0\\0 & - 4.0 k \nu K{\left(t \right)} \cos{\left(k z \right)} & 0\\0 & 0 & - 4 k K{\left(t \right)} \cos{\left(k z \right)}\end{matrix}\right]$

$S^*$ is the on-fault displacement.

```python
S1_star_fun = -(uzm_star - uzp_star).subs(x, 0).simplify()
S1_star_fun
```

$\tau^*$ is the on fault shear traction.

```python
s13 = mu * (uxm_star.diff(z) + uzm_star.diff(x))
tau1_star_fun = -s13.subs(x, 0).simplify()
tau1_star_fun
```

Check continuity:

```python
s13 = mu * (uxp_star.diff(z) + uzp_star.diff(x))
tau1_star_fun = -s13.subs(x, 0).simplify()
tau1_star_fun
```

$\sigma_n^*$ is the on fault normal stress.

```python
s11 = (lamda + 2*mu) * uxm_star.diff(x) + lamda * uzm_star.diff(z)
sn_star_fun = s11.subs(x, 0).simplify()
sn_star_fun
```

Check continuity:

```python
s11 = (lamda + 2*mu) * uxp_star.diff(x) + lamda * uzp_star.diff(z)
sn_star_fun = s11.subs(x, 0).simplify()
sn_star_fun
```

Slip-rate $V^*$

```python
V1_star_fun = S1_star_fun.diff(t)
V1_star_fun
```

Acceleration $\frac{dV^*}{dt}$

```python
V1_star_fun.diff(t)
```

The following lines define the functions $K$ and necessary time derivatives.

```python
K_fun = 2.0 * (atan((t-te)/tw) + pi/2) / pi
K_fun
```

```python
K_fun.subs(t, 0)
```

```python
K_fun.subs(t, te)
```

```python
K_fun.diff(t)
```

```python
K_fun.diff(t).diff(t).simplify()
```

```python

```

max_stars_repo_name: NicoSchlw/tandem

max_issues_repo_name: NicoSchlw/tandem

max_forks_repo_name: NicoSchlw/tandem


text:
```python
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sympy import diff, Symbol
import numpy as np

"""
Basic version of Newton-Raphson Algorithm with symbolic math.

Just define the number of variables, the function f(z1,z2,...,zn) and the algorithm
does the calculation.

Issues: Convergence is not always granted. With more than 2 variables the plot of the
third cell is not available.
"""

n = 3 # numero de variables
z = [Symbol('z'+str(i+1)) for i in range(n)] # definicion de variables z1,z2,...,zn

f = 0.6*z[0]**4+2*z[1]**4+3*z[0]**3- 2*z[1]**2 +30*z[0] +3 #definicion de funcion f(z)
f = (z[0]-2)**2 + (z[1]-1)**2 -z[2]*(z[0]+z[1]-1)#definicion de funcion f(z)

#f = z[0]**2 + 2*z[1]**2 #definicion de funcion f(z)

```

```python
np.set_printoptions(formatter={'float': lambda x: "{0:0.3f}".format(x)})

gradient=[f.diff(zi) for zi in z] # calculo del gradiente f'(z)
hessian=[[f.diff(zi).diff(zj) for zi in z] for zj in z]

x=[np.array([0.5,0.5,-1])] #iteraciones del metodo del gradiente, valor inicial x[0]=(2,1)

#t=0.1 #valor del paso fijo
maxiter=1000 #numero maximo de iteraciones

print("Initial Point: "+str(x[0]))

for k in range(maxiter):
    nablaf_xk=np.array([nf.evalf(subs = dict(zip(z,x[k]))) for nf in gradient],dtype="float32") #evaluar x[k] en gradiente
    hessianf_xk=np.matrix([[hij.evalf(subs = dict(zip(z,x[k]))) for hij in hi] for hi in hessian],dtype="float32") #evaluar x[k] en matriz hessiana

    x.append(x[k]-np.tensordot(np.linalg.inv(hessianf_xk),nablaf_xk,1))

    if np.linalg.norm(nablaf_xk)<0.00001: # termino del algoritmo
        print("Optimo: "+str(x[k+1]))
        break
    else:
        print("Iteration "+str(k+1)+", gradient: "+str(nablaf_xk)+", x: "+str(x[k+1]))
```

    Initial Point: [0.500 0.500 -1.000]
    Iteration 1, gradient: [-2.000 0.000 0.000], x: [1.000 0.000 -2.000]
    Optimo: [1.000 0.000 -2.000]

```python
import plotly.graph_objects as go

xd=np.linspace(-10,5,101)
yd=np.linspace(-5,5,101)
X,Y=np.meshgrid(xd,yd)
Z=np.array([[f.evalf(subs = dict(zip(z,[x,y]))) for x in xd] for y in yd]).astype(float)
fig = go.Figure()
fig.add_trace(go.Surface(z=Z, x=X, y=Y, opacity=0.3))

xp=np.round([p[0] for p in x])
yp=np.round([p[1] for p in x])
zp=np.array([f.evalf(subs = dict(zip(z,[p[0],p[1]]))) for p in x]).astype(float)

fig.add_trace(go.Scatter3d(x=xp, y=yp, z=zp,marker=dict(size=4,color="red"),
                                line=dict(color='darkblue',width=2)))

fig.update_layout(title='Newton-Raphson Optimization', autosize=True,
                  width=1000, height=1000,
                  margin=dict(l=65, r=50, b=65, t=90))

fig.show()

```

<meta charset="utf-8" />

max_stars_repo_name: maxigaarp/teaching

max_issues_repo_name: maxigaarp/teaching

max_forks_repo_name: maxigaarp/teaching