text:
# License
https://raw.githubusercontent.com/computational-sediment-hyd/1DBedVariationMixed/master/LICENSE

# import module

```python
import numpy as np
import pandas as pd
import numba
import sys
import os
import json
import glob as glob
```

```python
screenclass = np.array( [ 2.0**0, 2.0**1, 2.0**2, 2.0**3, 2.0**4, 2.0**5, 2.0**6] )
dsize = np.array( [ 10**(0.5*(np.log10(screenclass[i]) + np.log10(screenclass[i+1]))) for i in range(len(screenclass)-1) ] )

h=1.0
ib = 1.0/200.0
hini = 1.0
n = 0.03
v = ib**0.5 * h**(2/3) / n
fr = v/np.sqrt(9.8*h)
print('Froude Number:' + str(fr))

columns = ['diameter[mm]' ]
df = pd.DataFrame(dsize, columns=columns)
df['tau_star'] = h*ib/dsize/1.65*1000
df
```

    Froude Number:0.7529232524210427

      diameter[mm]
      tau_star

      0
      1.414214
      2.142748

      1
      2.828427
      1.071374

      2
      5.656854
      0.535687

      3
      11.313708
      0.267843

      4
      22.627417
      0.133922

      5
      45.254834
      0.066961

# Governing Equation of River flow

\begin{align}
    \frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x} = 0 \\
\end{align}

\begin{align}
    \frac{\partial Q}{\partial t}+\frac{\partial}{\partial x}\left(\frac{ Q^2}{ A}\right) + gA\frac{\partial H}{\partial x}+gAi_e = 0
\end{align}

```python
@numba.jit(nopython=True, parallel=False)
def flowmodel(A, Q, Adown, Qup, dAb, zb, B, dx, dt, g, manning):
    imax = len(A)
    Anew, Qnew = np.zeros(imax), np.zeros(imax)

# continuous equation
    for i in numba.prange(1, imax-1) : Anew[i] = A[i] - dt * ( Q[i] - Q[i-1] ) / dx

    Anew[0], Anew[-1] = Anew[1], Adown

# moumentum equation
    for i in numba.prange(1,imax-1):
        ip, ic, im = (i+1, i, i-1)
        Cr1 = 0.5*( Q[ic]/A[ic] + Q[ip]/A[ip] )*dt/dx
        Cr2 = 0.5*( Q[ic]/A[ic] + Q[im]/A[im] )*dt/dx
        dHdx1 = ( Anew[ip]/B[ip] + zb[ip] + dAb[ip]/B[ip] - Anew[ic]/B[ic] - zb[ic] - dAb[ic]/B[ic] ) / dx
        dHdx2 = ( Anew[ic]/B[ic] + zb[ic] + dAb[ic]/B[ic] - Anew[im]/B[im] - zb[im] - dAb[im]/B[im] ) / dx
        dHdx = (1.0 - Cr1) * dHdx1 + Cr2 * dHdx2

        Qnew[ic] = Q[ic] - dt * ( Q[ic]**2/A[ic] - Q[im]**2/A[im] ) / dx  \
                   - dt * g * Anew[ic] * dHdx \
                   - dt * g * A[ic] * manning**2 * Q[ic]**2 / B[ic]**2 / ( A[ic]/B[ic] )**(10.0/3.0)

    Qnew[0], Qnew[-1] = Qup, Qnew[-2]

    return Anew, Qnew
```

# Governing Equation of River Bed : Hirano model

\begin{align}
    (1-\lambda)\frac{\partial A_{b}}{\partial t}+\frac{\partial }{\partial x} \sum_{i=1}^n ( Q_{bi}P_i)= 0 \\
\end{align}

\begin{align}
    \frac{\partial P_i}{\partial t} &= - \frac{1}{(1-\lambda)E_dB} \frac{\partial (Q_{bi}P_i)}{\partial x} - \frac{P_{si}}{E_d B}\frac{\partial A_b}{\partial t} \\
    &= - \frac{1}{E_d B}\left(\frac{\partial A_{bi}}{\partial t} + P_{si}\frac{\partial A_b}{\partial t}\right)
\end{align}

\begin{align}
&P_{si} =
\left\{ \begin{array}{ll}
    P_i \mbox{ in exchange layer} & \left(\dfrac{\partial A_b}{\partial t} > 0 \right) \\
    P_i \mbox{ under exchange layer} & \left(\dfrac{\partial A_b}{\partial t} < 0 \right) \\
\end{array} \right. \\
\end{align}

$Q_{bi}$ : Ashida-Michiue Eq. and  Egiazaroff Eq.

$E_d$ : thickness of exchange layer

\begin{align}
    \frac{\tau_{*c i}}{\tau_{*cm}} &= 0.85 \frac{D_m}{D_i} \qquad & \left(\frac{D_i}{D_m} < 0.4 \right) \\
    \dfrac{\tau_{*ci}}{\tau_{*cm}} &= \left( \dfrac{ \displaystyle \log_e 19 }{ \displaystyle \log_e \left(19\dfrac{D_i}{D_m} \right) } \right)^2 \qquad & \left(\frac{D_i}{D_m} \geq 0.4 \right)
\end{align}

```python
@numba.jit(nopython=True, parallel=False)
def sedimentmodel(dAb, dratio, A, Q, B, dx, dt, g, manning, dsize, dratioStandard, hExlayer, Qbup):
    rhosw = 1.65 # grain specific gravity in water
    porosity = 0.4
    tscAve = 0.05 # critical tractive force of average grain size

    dAbnew = np.zeros_like( dAb )
    drationew = np.zeros_like( dratio )
    Qbsub = np.zeros_like( dratio )
    dAbsub = np.zeros_like( dratio )

    imax, lmax = len(dratio), len(dratio[0])

    for i in numba.prange(imax):
        Ap, Qp, dr, Bp = A[i], Q[i], dratio[i], B[i]
        dAve = np.sum(dsize * dr)
        us = np.sqrt(g * Ap/Bp) * Qp/Ap * manning / (Ap/Bp)**(2/3)
        tsAve = us**2.0/rhosw/g/dAve
        use = Qp/Ap / ( 6.0 + 2.5 * np.log( Ap/Bp/dAve/( 1.0+2.0*tsAve) ) )
        Kc=1.0

        for l in range(lmax):
            dri, dsi = dr[l], dsize[l]
            ts  = us**2.0 /rhosw/g/dsi
            tse = use**2.0/rhosw/g/dsi

# Egiazaroff Eq.
            x = dsi/dAve
            tscbytscm = (np.log(19)/np.log(19*x))**2 if x > 0.4 else 0.85/x
            tsc = Kc * tscbytscm * tscAve

            if ts > tsc :
# Ashida-Michiue Eq.
                Qbsub[i][l] = np.sqrt(rhosw * g * dsi**3.0) \
                        * 17.0 * tse**1.5 * ( 1.0 - tsc / ts) * ( 1.0 - np.sqrt(tsc / ts) ) \
                        * dri * B[i]
            else:
                Qbsub[i][l] = 0.0

    for i in numba.prange(imax):
        Qbs = Qbsub[i]
        if i == 0:
            if np.min(Qbup) < 0.0 : # equilibrium sand supply
                dAbsub[i][:] = 0.0
            else:
                for l in range(lmax):
                    dAbsub[i][l] = - dt / (1.0 - porosity) * ( Qbsub[i][l]-Qbup[l] )/dx
                    if np.abs(dAbsub[i][l]) > hExlayer :
                        print(i) ; print('dzsub-error')
        else:
            for l in range(lmax):
                dAbsub[i][l] = - dt / (1.0 - porosity) * ( Qbsub[i][l] - Qbsub[i-1][l] )/dx
                if np.abs(dAbsub[i][l]) > hExlayer :
                    print(i) ; print('dzsub-error')

    for i in numba.prange(imax):
        Qbs = Qbsub[i]
# update dAb
        dAball = np.sum( dAbsub[i] )
        if np.abs(dAball) > hExlayer:
            print(i) ; print('dz-error')

        dAbnew[i] = dAb[i] + dAball

# update dratio
        dratioIn = dratioStandard[i][:] if dAball < 0.0 else dratio[i][:]

        for l in range(lmax):
            drationew[i][l] = dratio[i][l] + ( dAbsub[i][l] - dAball * dratioIn[l] ) /hExlayer/B[i]
            if drationew[i][l] < 0.0 : drationew[i][l] = 0.0

# correct ration so that sum of ration become 100%
        sumdr = np.sum( drationew[i] )
        drationew[i] /= sumdr

    return dAbnew, drationew, Qbsub
```

# main

```python
def bedvariation(
dx,dt,manning,totalTime,outTimeStep,RunUpTime
,dsize ,dratioStandard ,dratio ,hExlayer ,A ,Q ,B ,zb ,dAb ,Qup ,Adown
,outputfilename
):

    g = 9.8
    Qbup = np.full( ( len(dsize) ), -9999.0 ) # when equilibrium sand supply, set qbup to minus value

# run-up calculation
    for i in range(int(RunUpTime/dt)):
        ib = ( ( dAb[-2]/B[-2] + zb[-2] ) - ( dAb[-1]/B[-1] + zb[-1] ) )/dx
        Adownp = Adown(0.0, Q[-1], dAb[-1]/B[-1], ib)
        Qupp = Qup(0.0)
        A, Q = flowmodel(A, Q, Adownp, Qupp, dAb, zb, B, dx, dt, g, manning)

# main calculation
    for i in range( int(totalTime/dt) ):

# cal bed variation
        dAb, dratio, Qbsub = sedimentmodel(dAb, dratio, A, Q, B, dx, dt, g, manning, dsize, dratioStandard, hExlayer, Qbup)

# cal water profile
        ib = ( ( dAb[-2]/B[-2] + zb[-2] ) - ( dAb[-1]/B[-1] + zb[-1] ) )/dx
        Adownp = Adown(i*dt, Q[-1], dAb[-1]/B[-1], ib)
        Qupp = Qup(i*dt)
        A, Q = flowmodel(A, Q, Adownp, Qupp, dAb, zb, B, dx, dt, g, manning)

# output
        if( int(i*dt) % int(outTimeStep) ) == 0 :
            print( str( int(i*dt) ) + ' second')
            profile = []
            for ii, (ap, qp, z, r, q) in enumerate(zip(A, Q, dAb, dratio, Qbsub)) :
                profile.append({'distance' : int(ii*dx), 'A':ap, 'Q':qp, 'dAb':z, 'ratio':list(r), 'Qb':list(q)})
            json.dump( {'time':int(i*dt), 'profile' : profile}, open('%010d' % int(i*dt) + 'sec.json', 'w') )

# join json files
    d = [ json.load( open(f, 'r') ) for f in glob.glob('*sec.json') ]
# delete json files
    for f in glob.glob('*sec.json') : os.remove(f)

    d.sort(key=lambda x: x['time'])

    cond = {'width':list(B), 'elevation':list(zb), 'screenclass':list(screenclass), 'dsize':list(dsize) }
    json.dump( {'condition':cond, 'output':d}, open(outputfilename, 'w') )
```

# Sample
## case1

```python
%%time
length = 10000.0
dx = 100.0
imax = int(length/dx) + 1
dt = 10.0
totalTime = 1000.1*3600.0
outTimeStep = 100.0*3600.0
RunUpTime = 3.0 * 3600.0
hini = 1.0
manning  = 0.03
ib = 1.0/200.0
ib2 = 1.0/1000.0
outputfilename = 'case01.json'

# grain diameter classification
screenclass = np.array( [ 2.0**0, 2.0**1, 2.0**2, 2.0**3, 2.0**4, 2.0**5, 2.0**6], dtype=float )/1000
dsize = np.array( [ 10**(0.5*(np.log10(screenclass[i]) + np.log10(screenclass[i+1]))) for i in range(len(screenclass)-1) ], dtype=float )

# percentage of grain size under exchange layer
dratioStandard1 = np.full_like(dsize, 1.0/len(dsize), dtype=float)
dratioStandard = np.full( (imax, len(dsize) ), dratioStandard1, dtype=float)

# initial percentage of grain size in exchange layer
dratio = np.copy(dratioStandard)

# thickness of exchange layer
hExlayer = dsize[-1]

# Initial & Boundary condition
B = np.full(imax, 1.0, dtype=float)
A = hini*B
Q = ib**0.5*(hini)**(5.0/3.0)/manning*B #normal flow
zb = np.zeros(imax)
for i in range(1,imax):
    zb[i] = zb[i-1] + ib2*dx if i < 50 else zb[i-1] + ib*dx
zb = zb[::-1]

dAb = np.zeros(imax)

Qup = Q[0]

def Adown(time, Q, dzb, ib):
    return ( manning**2*Q**2/ib/B[-1]**2 )**0.3 * B[-1]

def Qup(time):
    return Q[0]

bedvariation(
dx,dt,manning,totalTime,outTimeStep,RunUpTime
,dsize ,dratioStandard ,dratio ,hExlayer ,A ,Q ,B ,zb ,dAb ,Qup ,Adown
,outputfilename
    )
```

    0 second
    360000 second
    720000 second
    1080000 second
    1440000 second
    1800000 second
    2160000 second
    2520000 second
    2880000 second
    3240000 second
    3600000 second
    Wall time: 37.2 s

## case2

```python
%%time
length = 10000.0
dx = 100.0
imax = int(length/dx) + 1
dt = 10.0
totalTime = 3000.1*3600.0
outTimeStep = 300.0*3600.0
RunUpTime = 3.0 * 3600.0
hini = 1.0
manning  = 0.03
ib = 1.0/200.0
ib2 = 1.0/1000.0
outputfilename = 'case02.json'

# grain diameter classification
screenclass = np.array( [ 2.0**0, 2.0**1, 2.0**2, 2.0**3, 2.0**4, 2.0**5, 2.0**6], dtype=float )/1000
dsize = np.array( [ 10**(0.5*(np.log10(screenclass[i]) + np.log10(screenclass[i+1]))) for i in range(len(screenclass)-1) ], dtype=float )

# percentage of grain size under exchange layer
dratioStandard1 = np.full_like(dsize, 1.0/len(dsize), dtype=float)
dratioStandard2 = np.zeros(len(dsize))
dratioStandard2[-1] = 1.0

dratioStandard = np.full( (imax, len(dsize) ), dratioStandard1, dtype=float)
dratioStandard[:50] = dratioStandard2

# initial percentage of grain size in exchange layer
#     dratio = np.full( (imax, len(dsize) ), dratioStandard1, dtype=float)
dratio = np.copy(dratioStandard)

# thickness of exchange layer
hExlayer = dsize[-1]

# Initial & Boundary condition
B = np.full(imax, 1.0, dtype=float)
A = hini*B
Q = ib**0.5*(hini)**(5.0/3.0)/manning*B #normal flow
zb = np.zeros(imax)
for i in range(1,imax):
    zb[i] = zb[i-1] + ib2*dx if i < 50 else zb[i-1] + ib*dx
zb = zb[::-1]

dAb = np.zeros(imax)

Qup = Q[0]

def Adown(time, Q, dzb, ib):
    return ( manning**2*Q**2/ib/B[-1]**2 )**0.3 * B[-1]

def Qup(time):
    return Q[0]

bedvariation(
dx,dt,manning,totalTime,outTimeStep,RunUpTime
,dsize ,dratioStandard ,dratio ,hExlayer ,A ,Q ,B ,zb ,dAb ,Qup ,Adown
,outputfilename
    )
```

    0 second
    1080000 second
    2160000 second
    3240000 second
    4320000 second
    5400000 second
    6480000 second
    7560000 second
    8640000 second
    9720000 second
    10800000 second
    Wall time: 2min 31s

## case3

```python
%%time
length = 10000.0
dx = 100.0
imax = int(length/dx) + 1
dt = 10.0
totalTime = 3000.1*3600.0
outTimeStep = 300.0*3600.0
RunUpTime = 3.0 * 3600.0
hini = 1.0
manning  = 0.03
ib = 1.0/1000.0
outputfilename = 'case03.json'

# grain diameter classification
screenclass = np.array( [ 2.0**0, 2.0**1, 2.0**2, 2.0**3, 2.0**4, 2.0**5, 2.0**6], dtype=float )/1000
dsize = np.array( [ 10**(0.5*(np.log10(screenclass[i]) + np.log10(screenclass[i+1]))) for i in range(len(screenclass)-1) ], dtype=float )

# percentage of grain size under exchange layer
dratioStandard1 = np.full_like(dsize, 1.0/len(dsize), dtype=float)
dratioStandard = np.full( (imax, len(dsize) ), dratioStandard1, dtype=float)

# initial percentage of grain size in exchange layer
dratio = np.copy(dratioStandard)

# thickness of exchange layer
hExlayer = dsize[-1]

# Initial & Boundary condition
B = np.full(imax, 100.0, dtype=float)
for i in range(40, 50):
    B[i] = B[i-1] - 5.0
for i in range(50, 60):
    B[i] = B[i-1] + 5.0

A = hini*B
Q = ib**0.5*(hini)**(5.0/3.0)/manning*B #normal flow
zb = np.zeros(imax)
for i in range(1,imax):
    zb[i] = zb[i-1] + ib*dx
zb = zb[::-1]

dAb = np.zeros(imax)

Qup = Q[0]

def Adown(time, Q, dzb, ib):
    return ( manning**2*Q**2/ib/B[-1]**2 )**0.3 * B[-1]

def Qup(time):
    return Q[0]

bedvariation(
dx,dt,manning,totalTime,outTimeStep,RunUpTime
,dsize ,dratioStandard ,dratio ,hExlayer ,A ,Q ,B ,zb ,dAb ,Qup ,Adown
,outputfilename
    )
```

    0 second
    1080000 second
    2160000 second
    3240000 second
    4320000 second
    5400000 second
    6480000 second
    7560000 second
    8640000 second
    9720000 second
    10800000 second
    Wall time: 1min 53s

## case4

```python
%%time
length = 10000.0
dx = 100.0
imax = int(length/dx) + 1
dt = 10.0
totalTime = 5000.1*3600.0
outTimeStep = 500.0*3600.0
RunUpTime = 3.0 * 3600.0
hini = 1.0
manning  = 0.03
ib = 1.0/500.0
outputfilename = 'case04.json'

# grain diameter classification
screenclass = np.array( [ 2.0**0, 2.0**1, 2.0**2, 2.0**3, 2.0**4, 2.0**5, 2.0**6], dtype=float )/1000
dsize = np.array( [ 10**(0.5*(np.log10(screenclass[i]) + np.log10(screenclass[i+1]))) for i in range(len(screenclass)-1) ], dtype=float )

# percentage of grain size under exchange layer
dratioStandard1 = np.full_like(dsize, 1.0/len(dsize), dtype=float)
dratioStandard = np.full( (imax, len(dsize) ), dratioStandard1, dtype=float)

# initial percentage of grain size in exchange layer
dratio = np.copy(dratioStandard)

# thickness of exchange layer
hExlayer = dsize[-1]

# Initial & Boundary condition
B = np.full(imax, 1.0, dtype=float)
A = hini*B
Q = ib**0.5*(hini)**(5.0/3.0)/manning*B #normal flow
zb = np.zeros(imax)
for i in range(1,imax):
    zb[i] = zb[i-1] + ib*dx
zb = zb[::-1]

WLdown = 5.0
for i, (Ap, zbp) in enumerate(zip(A, zb)):
    if( Ap/B[i] + zbp < WLdown) : A[i] = (WLdown - zbp)*B[i]

dAb = np.zeros(imax)

Qup = Q[0]

def Adown(time, Q, dzb, ib):
    return (WLdown - (dzb + zb[-1]))*B[-1]
#     return ( manning**2*Q**2/ib/B[0]**2 )**0.3 * B[0]

def Qup(time):
    return Q[0]

bedvariation(
dx,dt,manning,totalTime,outTimeStep,RunUpTime
,dsize ,dratioStandard ,dratio ,hExlayer ,A ,Q ,B ,zb ,dAb ,Qup ,Adown
,outputfilename
    )
```

    0 second
    1800000 second
    3600000 second
    5400000 second
    7200000 second
    9000000 second
    10800000 second
    12600000 second
    14400000 second
    16200000 second
    18000000 second
    Wall time: 3min 46s


text:
# Weight initialization

Based on lectures from Stanford's CS231n Training Neural Networks, Part 1.

**Question:** *What happens when we start with inital weghts set to 0, i.e., `W = 0`?*

**Answer:** If every neuron in the network computes the same output, then they will also all compute the same gradients during backpropagation and undergo the same parameter updates. In other words, there is no source of asymmetry between neurons if their weights are initialized to be the same value.

Therefore, one should initialize `W` with random weights to promote difference between neuron outputs. However, given a deep neural network, we learned that having really high weigts can lead to the *exploding gradient* problem and having really low weights can read to the *valishing gradient* problem.

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

def init_weight_plot(init_func, act_name):
    # assume some unit Guassian 10-D input data
    D = np.random.randn(1000, 500)
    hidden_layer_sizes = [500] * 10
    nonlinearities = [act_name] * len(hidden_layer_sizes)

    act = {'relu' : lambda x : np.maximum(0, x), 'tanh' : lambda x : np.tanh(x)}
    Hs = {}
    for i in range(len(hidden_layer_sizes)):
        X = D if i == 0 else Hs[i - 1] # input at this layer
        fan_in = X.shape[1]
        fan_out = hidden_layer_sizes[i]
        W = init_func(fan_in, fan_out)

        H = np.dot(X, W) # matrix multiply
        H = act[nonlinearities[i]](H) # nonlinearity
        Hs[i] = H

    # look at distributions at each layer
    print('input layer had mean %f and std %f' % (np.mean(D), np.std(D)))
    layer_means = [np.mean(H) for i, H in Hs.items()]
    layer_stds = [np.std(H) for i, H in Hs.items()]
    for i, H in Hs.items():
        print('hidden layer %d had mean %f and std %f' % (i + 1, layer_means[i], layer_stds[i]))

    # plot the means and standard devations
    plt.figure()
    plt.subplot(121)
    plt.plot(Hs.keys(), layer_means, 'ob-')
    plt.title('layer mean')
    plt.subplot(122)
    plt.plot(Hs.keys(), layer_stds, 'or-')
    plt.title('layer std')

    # plot the raw distributions
    fig, axes = plt.subplots(nrows = 1, ncols = 10, figsize=(15, 5))
    for i, ax in enumerate(axes.flat, start = 1):
        ax.hist(Hs[i - 1].ravel(), 30, range = (-1, 1))

    fig.tight_layout()
    plt.show()
```

```python
func = lambda x, y: np.ones((x, y)) * 0.02
init_weight_plot(func, 'tanh')
```

## Option 1: Initialize using small random weights

One common practice is to initialize the weights of the neurons to small numbers, which is also referred to as *breaking symmetry*. The idea is that the neurons will start as random and unique, and will thus compute distinct updates and integrate themeselves as diverse parts of the full network. A default choice is to use Gausian Normal values scaled by 1e-2.

$W = 0.01 * np.random.randn(n_{l-1},n_{l})$

where $n_{l-1}$ is the number of inputs and $n_{l}$ is the number of hidden units for the $l$-th layer.

This tends to work well for small networks, but it is problematic still for deep networks, potentially leading to *valishing gradients* and dead nodes.

Let's see an example. First, we will initialize with small random weights.

```python
func = lambda x, y: np.random.randn(x, y) * 0.01
init_weight_plot(func, 'tanh')
```

**What happened?** Weights in the last layer are all 0 and we cannot make any useful prediction.

**Why did this happen?** Think about the backprop step. What do the gradients look like for one of the layers?

**Answer** Small weights will lead to small gradients (since the gradient is proportional to the value of the weights). This greatly diminshes the 'gradient signal' flowing backward through the network, and makes it hard to have clear predictions.

**What if we just increase the weights?**

```python
func = lambda x, y: np.random.randn(x, y) * 1.0
init_weight_plot(func, 'tanh')
```

This makes the neurons become saturated from the start, simply propagating -1 or 1 through the network layers.

## Option 2: Xavier Initialization [Glorot et al., 2010]

Xavier Glorot and Yoshua Bengio showed that, in the case of linear activations, the distribution of the outputs in a neuron with randomly initialized weights has a variance that grows with the number of inputs. The higher the number of inputs, the higher the variance. They thus suggest normalizing the variance of each neuron's output to 1 by scaling its weight vector by the square root of its fan-in (i.e., its number of inputs), as follows.

```python
func = lambda x, y: np.random.randn(x, y) / np.sqrt(x)
init_weight_plot(func, 'tanh')
```

**Why does it work?** Consider the dot product between the weight and input vector, $w^Tx$, which is the input to the non-linear activation function. For the case of simplicity, we'll ignore the bias term. Let us then examine the variance of the ouput values.

$
\begin{align}
Var(w^Tx) &= Var(\sum_i^n w_ix_i)\\
          &= \sum_i^n Var(w_ix_i)\\
          &= \sum_i^n [E(w_i)]^2 Var(x_i) + [E(x_i)]^2 Var(w_i) + Var(x_i)Var(w_i)\\
          &= \sum_i^n Var(x_i)Var(w_i)\\
          &= (nVar(w))Var(x)
\end{align}
$

The first 3 steps in the derivation are due to the properties of variance, while the 4th assumes that both the data and weights are mean centered (*note that will not be the case with outputs from ReLU activations*). The final step assumes an identical distribution of the $x_i$ and $w_i$ values.

If we wish to prevent exploding or vanishing gradients, we should strive for the variance of the activation inputs to be similar to the variance of the layer inputs. The consequence of the derivation above is that, if we want to maintain the variance of the inputs after the dot-product operation, the weights $w$ should be scaled in such a way that the $Var(w) = 1/n$, where $n$ is the number of inputs for the layer. Given that, for any scalar $a$ and random variable $X$, $Var(aX)=a^2Var(X)$, this implies that the initial weights should be unit Gaussian normal values scaled by $\sqrt{1/n}$, leading to weight vectors with variance $1/n$, just as we do in the Xavier Initialization.

## Xavier Initialization and ReLU activations

Unlike tanh gates, the output of the ReLU (and other activations) is non-negative, leading to non-zero positive mean, which breaks the assumptions above. And, as we will see, Xavier Initialization does not seem to work with ReLU nonlinearity.

```python
func = lambda x, y: np.random.randn(x, y) / np.sqrt(x)
init_weight_plot(func, 'relu')
```

## He Initialization for ReLU activations [He et al., 2015]

He et al. suggest that the variance of the neurons should be $2/n$ for ReLU activated nodes, rather than $1/n$, leading to initial weights drawn from a unit Gaussian normal and scaled by $\sqrt{n/2}$.

```python
# note the sqrt(x/2) instead of sqrt(x)
func = lambda x, y: np.random.randn(x, y) / np.sqrt(x/2)
init_weight_plot(func, 'relu')
```

ext: ipynb


text:
# Gaussian Processes

Based on code from: http://www.cs.ubc.ca/~nando/540-2013/lectures/gp.py

Also referenced:

* http://bridg.land/posts/gaussian-processes-1
* https://www.youtube.com/watch?v=4vGiHC35j9s
* "Gaussian Processes for Machine Learning" book

Key idea: as mentioned in these [two](https://www.youtube.com/watch?v=50Vgw11qn0o) [lectures](https://www.video.ethz.ch/conferences/2014/rnls/05_thursday/0f8953ea-b1fe-4960-8fa1-ef71bb642445.html), we use Gaussians to represent things because they make the math easy, and the assumptions we have to make to model real stuff as Gaussians are reasonable enough (same reason we assume linear dynamical systems). Gaussians turn marginalization and conditioning into linear algebra, which is why "Gaussian Processes are the linear algebra of inference". Simple matrix multiplications and inverses give us probabilistic inference of joint probability distributions.

Model $p(\mathbf{y}\vert\mathbf{x})$ using a multivariate normal:

$$p(\mathbf{y}\vert\mathbf{x})=\mathcal{N}(\mathbf{y}\vert m(\mathbf{x}),\mathbf{K})$$

where $\mathbf{K}=\kappa(\mathbf{x}, \mathbf{x})$ and $m(\mathbf{x})=\mathbf{0}$. When we assume zero-mean, the GP is completely defined by the covariance function $\mathbf{K}$.

This is a prior distribution representing the kind of outputs $\mathbf{y}$ that we expect to see over some inputs $\mathbf{x}$ _before_ we observe any data.

So we have some training data with inputs $\mathbf{x}$, and outputs $\mathbf{y}=f(\mathbf{x})$. Now lets say we have some new points $\mathbf{x}_*$ where we want to predict $\mathbf{y}_*=f(\mathbf{x}_*)$

GPs are non-parametric; as in, they don't specify a set of parameters for the function (e.g. linear regression $y=w^Tx$, neural networks/others $y=f(\theta, x)$). The way I think of it is: in GPs, every data point from the training set becomes a basis function. See Section 2.2 in Rasmussen book for illustration of this

Wikipedia explains it well [here](https://en.wikipedia.org/wiki/Nonparametric_statistics): "The second meaning of non-parametric covers techniques that do not assume that the structure of a model is fixed. Typically, the model grows in size to accommodate the complexity of the data", "The term non-parametric is not meant to imply that such models completely lack parameters but that the number and nature of the parameters are flexible and not fixed in advance." This means that unlike parametric methods like neural networks, GPs can grow in flexibility as more data is added to effectively "memorize" all the training data.

Note: GPs are a specific instantiation of kernel method, where the Gaussian form allows closed-form solutions for most operations.

Another interpretation that made sense to me from Wikipedia: "A machine-learning algorithm that involves a Gaussian process uses lazy learning and a measure of the similarity between points (the kernel function) to predict the value for an unseen point from training data." "Ultimately Gaussian processes translate as taking priors on functions and the smoothness of these priors can be induced by the covariance function. If we expect that for "near-by" input points x and x' their corresponding output points y and y' to be "near-by" also, then the assumption of continuity is present. If we wish to allow for significant displacement then we might choose a rougher covariance function. Extreme examples of the behaviour is the Ornstein–Uhlenbeck covariance function and the squared exponential where the former is never differentiable and the latter infinitely differentiable."

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

```python
# Real function to approximate
# xs = np.linspace(-7, 7, 100)
# plt_range = [-7, 7, -3, 3]
# f = lambda x: np.sin(0.9*x).flatten()

# Another function to try: 5th order polynomial
# Note that you'll have to change bounds on the plots, as well as samples for training
# coefs[i] is the coefficient of x^i
xs = np.linspace(-5.0, 3.5, 100)
plt_range = [-5, 3.5, -10, 15]
coefs = np.array([6, -2.5, -2.4, -0.1, 0.2, 0.03])
def f(x):
    total = 0
    for exp, coef in enumerate(coefs):
        total += coef * (x ** exp)
    return total.flatten()

ys = f(xs)

fig, ax = plt.subplots(figsize=(12, 5))
plt.plot(xs, ys, '*')
plt.show()
```

### GP fit

The kernel function acts as a prior on our function. Here, we assume smoothness and encode this in the kernel:

$$ \kappa(x, x')=\sigma^2\exp\left(-~\frac{(x-x')^2}{2l}\right)$$

This function is 1 when $x=x'$ and tends to zero as its arguments drift apart. It encourages the gaussian process to have similar valued outputs for similar valued inputs. Note that this function is infinitely differentiable.

Note that this is where design choices and hyperparameter tuning fall in GPs; kernels encode assumptions/background knowledge about the data. The hyperparameters can also be optimized with maximization of the marginal likelihood (some form of gradient descent). See slide 23 [here](https://people.cs.umass.edu/~wallach/talks/gp_intro.pdf)

```python
"""
GP squared exponential kernel
"""
def se_kernel(a, b):
    kernelParameter = 0.5 # horizontal scale
    sigma = 1 # vertical variation
    sqdist = np.sum(a**2,1).reshape(-1,1) + np.sum(b**2,1) - 2*np.dot(a, b.T)
    return sigma**2 * np.exp(-0.5 * (1/kernelParameter) * sqdist)
```

```python
n_train = 5         # number of training points.
n = 50         # number of test points.
s = 5e-5    # noise variance.
```

```python
"Construct function from random sample/training points along reference function"
# Sample some input points and noisy version of the function evaluated at these points.
X = np.random.uniform(-5, 5, size=(n_train,1))
y = f(X) + s*np.random.randn(n_train)

K = se_kernel(X, X)
L = np.linalg.cholesky(K + s*np.eye(n_train))
```

We represent the covariance matrix with its cholesky decomposition because it makes later operations (sampling from distribution, inverting covariance matrix) more efficient.

The following math is from Rasmussen pg 16.

##### Computing the mean at test points

Want to do: $\mu = K_* K^{-1} y$, where $K=k(X_\text{train}, X_\text{train})$ and $K_*=k(X_\text{train}, X_\text{test})$

Trying to find an intuition for this:

* $y$ $n_\text{train} \times n_\text{dims}$ (training data)
* $K$ $n_\text{train} \times n_\text{train}$ (covariance matrix of training data - inverse is information matrix)
    * $K^{-1}y$ : extract salient features from y?
* $K_*$ $n_\text{test} \times n_\text{train}$ (covariance between test and train data)
    * project what we know about the training data onto the test set

$K = L L^T \rightarrow K^{-1} = L^{-T} L^{-1}$

$L_k = L^{-1} K_*$

$L_k^T = K^{T}_* L^{-T}$

$\mu = K_* L^{-T} L^{-1} y = K_* K^{-1} y$

##### Computing the variance at test points

TODO clarify this

$K\_ = k(X_\text{test}, X_\text{test})$

$diag(K\_)$ is vector of variances

$L_k^2$ is element-wise square

$L_k = L^{-1} K_*$

$L_k^2 = (L^{-1})^2 (K_*)^2$

```python
# points we're going to make predictions at.
Xtest = np.linspace(-5, 5, n).reshape(-1,1)

# compute the mean at our test points.
Lk = np.linalg.solve(L, se_kernel(X, Xtest))
mu = np.dot(Lk.T, np.linalg.solve(L, y))

# compute the variance at our test points
# Lk is size (10, 50), sum returns (50,)
K_ = se_kernel(Xtest, Xtest)
s2 = np.diag(K_) - np.sum(Lk**2, axis=0) #math trick
s = np.sqrt(s2)

# Compute full covariance matrix, then take diagonal to get variance
# Produces same results as calculation above
# cov(all points) = K** - (K* K^{-1} K*)
# K is covariance of GP, defined by kernel
# Note that this is less inefficient and inaccurate; doesn't use Cholesky
# See below for correct way to compute this
Kinv = np.linalg.inv(K)
K2 = K_ - (np.matmul(se_kernel(Xtest, X), np.matmul(Kinv, se_kernel(X, Xtest))))

# Note that doing it the first way results in some numerical errors
K2 = K_ - np.dot(Lk.T, Lk) + 1e-6*np.eye(n)

s2_2 = np.diag(K2)
s_2 = np.sqrt(s2_2)
```

```python
# Trying to make plot that is similar to plots in papers for DMPs and POMDPs, by summing up Gaussian basis functions
from scipy.stats import norm

"utility function for plotting gaussians around each data point"
def plot_gaussian(coord, mean, std, zero_height=0):
    x = np.linspace(norm.ppf(0.0001), norm.ppf(0.9999), 100)
    plt.plot(x+coord, mean*norm.pdf(x, 0, std)+zero_height, 'g-', lw=0.5)
```

```python
# PLOTS:
plt.plot(X, y, 'r.', ms=20)
plt.plot(Xtest, f(Xtest), 'r-')
plt.plot(Xtest, mu, 'b--', lw=2)
plt.gca().fill_between(Xtest.flat, mu-3*s, mu+3*s, color="#dddddd") #variances from shortcut trick
# plt.gca().fill_between(Xtest.flat, mu-3*s_2, mu+3*s_2, color="#ddddd0") #variances from full covariance matrix
plt.title('Mean predictions plus 3 st.deviations')
plt.axis(plt_range)

# Plot Gaussians at each test point, to show interpretation that GP mean is
# weighted sum of Gaussians from each data point
xg_total = 0
yg_total = 0
for i,x in enumerate(X):
    coord = x[0]
    mean = y[i]
    var = K2[i, i]
    plot_gaussian(coord, mean, np.sqrt(var), 0)

plt.legend(['Training points', 'Ground truth function', 'Predicted mean'])
```

Note that each test point (red plus signs) becomes a Gaussian "basis function" for the mean of the GP (green bell curves). The mean is the sum of all of these Gaussians (this is how DMPs, ProMPs work).

---

##### Aside: sampling from distribution specified by covariance

The covariance matrix is $\sigma^2$, and the Cholesky decomposition gives you the "square root" of the (positive definite) original matrix.

$$A = LL^T$$
Where $L$ is lower triangular with positive diagonal elements.

Recall that the definition of the covariance matrix is the outer product $\Sigma = (x-\mu)(x-\mu)^T$

From wikipedia: "The Cholesky decomposition is commonly used in the Monte Carlo method for simulating systems with multiple correlated variables. The correlation matrix is decomposed, to give the lower-triangular L. Applying this to a vector of uncorrelated samples u produces a sample vector Lu with the covariance properties of the system being modeled."

You can sample random points from a Gaussian distribution given $(\mu, \Sigma)$, using the Cholesky decomposition. The random points can be sampled with:
$$ p = \mu + LX$$
Where $X$ is a matrix of size $n \times k$, $k$ being the number of samples. Each column of $X$ is drawn from a normal distribution $N(0,1)$

The covariance of the points generated by this method are calculated as:
$$\Sigma_\text{sample} = (LX)(LX)^T = LXX^TL^T = L(XX^T)L^T$$
$XX^T \approx I$ because the vectors are taken from a normal distribution; therefore, $\Sigma_\text{sample} = \Sigma$

```python
# Testing random samples from distribution with Cholesky
np.set_printoptions(suppress=True, precision=4)

Sigma = np.array([[2, 0.5], [0.5, 1]])
print('Original covariance:\n', Sigma)
L_rand = np.linalg.cholesky(Sigma)

np.random.seed()
num_samples = 100000
Rand = np.random.normal(size=(2, num_samples))
p = np.matmul(L_rand, Rand) # assume zero mean

print('Sample covariance:\n', np.cov(p))
print('Random seed covariance (expect identity): \n', np.cov(Rand))
```

    Original covariance:
     [[2.  0.5]
     [0.5 1. ]]
    Sample covariance:
     [[1.998  0.5019]
     [0.5019 1.0052]]
    Random seed covariance (expect identity):
     [[0.999  0.0018]
     [0.0018 1.0047]]

---

##### Multivariate Gaussian Theorem

See Rasmussen chapter 2

$$\mu=\begin{bmatrix}\mu_1 \\ \mu_2 \end{bmatrix}, \Sigma = \begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22}\end{bmatrix}$$

Marginals:
$$p(x_1|x_2)=N(x_1 \vert \mu_{1|2}, \Sigma_{1|2})$$
$$\mu_{1|2} = \mu_1 + \Sigma_{12}\Sigma_{22}^{-1}(x_2-\mu_2)$$
$$\Sigma_{1|2} = \Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}$$

This section from: http://bridg.land/posts/gaussian-processes-1

Now recalling the definition of a GP, we will model the joint distribution of all of $\mathbf{y}$ and $\mathbf{y}_*$ as:

$$ \begin{pmatrix}\mathbf{y} \\ \mathbf{y}_*\end{pmatrix} \sim\mathcal{N}\left(\begin{pmatrix}m(\mathbf{x})\\ m(\mathbf{x}_*)\end{pmatrix}, \begin{pmatrix}\mathbf{K} & \mathbf{K}_* \\ \mathbf{K}_*^T & \mathbf{K}_{**}\end{pmatrix}\right) $$

where $\mathbf{K}=\kappa(\mathbf{x}, \mathbf{x})$, $\mathbf{K}_* = \kappa(\mathbf{x}, \mathbf{x}_*)$ and $\mathbf{K}_{**}=\kappa(\mathbf{x}_*, \mathbf{x}_*)$. As before we are going to stick with a zero mean.

However this is modeling $p(\mathbf{y}, \mathbf{y}_*\vert \mathbf{x}, \mathbf{x}_*)$ and we only want a distribution over $\mathbf{y}_*$!

We want: $p(\mathbf{y}_* \ \vert \  m(\mathbf{x}_*), \mathbf{K}_*)$

Unknowns: $m(\mathbf{x}_*), \mathbf{K}_*$

Knowns: $\mathbf{x}, \mathbf{x}_*, m(\mathbf{x}), \mathbf{y}, \mathbf{K}, \mathbf{K}_*, \mathbf{K}_{**}$.

TODO try to color equation above to distinguish knowns/unknowns

### Conditioning Multivariate Gaussians

Rather than deriving it from scratch we can just make use of [this standard result](https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Conditional_distributions). If we have a joint distribution over $\mathbf{y}$ and $\mathbf{y}_*$ as above, and we want to condition on the data we have for $\mathbf{y}$ then we have the following:

$$\begin{align}
p(\mathbf{y}_*\vert \mathbf{x}_*, \mathbf{x}, \mathbf{y})&=\mathcal{N}(\mathbf{y}_*\vert \mu_*, \Sigma_*)\\
\mu_*&=m(\mathbf{x}_*)+\mathbf{K}_*^T\mathbf{K}^{-1}(\mathbf{y}-m(\mathbf{x}))\\
\Sigma_*&=\mathbf{K}_{**}-\mathbf{K}_*^T\mathbf{K}^{-1}\mathbf{K}_*
\end{align}$$

Now we have a posterior distribution over $\mathbf{y}_*$ using a prior distribution and some observations!

```python
# draw samples from the prior at our test points.
L = np.linalg.cholesky(K_ + 1e-6*np.eye(n))
f_prior = np.dot(L, np.random.normal(size=(n,10)))

plt.plot(Xtest, f_prior)
plt.title('Ten samples from the GP prior')
plt.axis(plt_range)
```

```python
# Draw samples from the posterior at our test points.

# Calculate posterior covariance matrix
# K_posterior = K(Xtest,Xtest) - K(Xtest,X)K^{-1}K(X,Xtest)
# Add small value on diagonal to make sure it's positive definite
L = np.linalg.cholesky(K_ - np.dot(Lk.T, Lk) + 1e-6*np.eye(n))
f_post = mu.reshape(-1,1) + np.dot(L, np.random.normal(size=(n,10)))

plt.plot(Xtest, f(Xtest), 'k-', lw=4) # Real function
plt.plot(Xtest, f_post)
plt.title('Ten samples from the GP posterior')
plt.axis(plt_range)
plt.legend(['Ground truth'])
```

### Moving forward: stuff to consider

* How to pick kernel? This always needs to be hand-picked
* How to pick kernel parameters? These reflect our prior belief on the data. There are methods to do this
* We also have an assumed noise variance in the GP formulation

### Modeling the Rastrigin function

The Rastrigin function is commonly used as an example for optimization problems where gradient descent fails to find a good minimum (global). Here, we'll just use it as a "wavy shape" to see if a Gaussian process can fit it.

In n-dimensional domain, it is defined by:
$$f(x) = An + \sum^n_{i=1}[x_i^2-Acos(2\pi x_i)]$$

Global minimum is $x=0$, where $f(x)=0$.

Set $A=10$, $n=2$, $x\in[-5.12, 5.12]$ to copy wikipedia's example

```python
from matplotlib import cm #colormap
import math #pi
from mpl_toolkits.mplot3d import Axes3D

%matplotlib inline
# so 3d plots are interactive
# %matplotlib notebook
# from ipykernel.pylab.backend_inline import flush_figures
```

```python
def rastrigin(X, A=3):
    return A + sum([(x**2 - A * np.cos(2 * math.pi * x)) for x in X])
```

```python
# TODO find a better way to do this. meshgrid?
def xy_combinations(X, Y):
    # Make a matrix of (x,y) coordinates for scatter3d for two vectors X,Y
    XY_plot = np.zeros((X.size*Y.size, 2)) # each row will be one point, so Nx2
    i = 0
    for x in X:
        for y in Y:
            XY_plot[i, :] = np.array([x, y])
            i += 1
    return XY_plot
```

```python
# Plot reference function
X_plot = np.linspace(-3, 3, 100)
Y_plot = np.linspace(-3, 3, 100)
X2, Y2 = np.meshgrid(X_plot, Y_plot)
Z = rastrigin(X2, Y2)

fig = plt.figure()
ax = fig.gca(projection='3d')
# ax.plot_surface(X2, Y2, Z, rstride=1, cstride=1, linewidth=0, antialiased=False)#, cmap=cm.plasma)

XY_plot = xy_combinations(X_plot, Y_plot)
Z = rastrigin(XY_plot.T)
ax.scatter3D(XY_plot[:,0], XY_plot[:,1], Z, c=Z, cmap='viridis')
```

```python
# Alternative plot: trisurf
# Creates surface by finding a set of triangles formed between adjacent points
# (x,y,z) can be 1D arrays instead of rectilinear. samples can be random
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_trisurf(XY_plot[:,0], XY_plot[:,1], Z, cmap='viridis')
```

```python
# grid training points
# N = 5
# X = np.linspace(-3, 3, N)
# Y = np.linspace(-3, 3, N)

# random training points
n_train = 30
X = np.random.uniform(-3, 3, n_train)
Y = np.random.uniform(-3, 3, n_train)

XY_train = xy_combinations(X, Y)
print(XY_train.shape) # number of points will be n_train^2
```

    (900, 2)

```python
# "Fit" GP
s = 5e-5    # scale parameter / noise variance prior
K = se_kernel(XY_train, XY_train)
L = np.linalg.cholesky(K + s*np.eye(K.shape[0]))

y = rastrigin(XY_train.T)

print(K.shape, L.shape, y.shape)

# Note that we only need to compute covariance because we're assuming zero mean
# Is this common practice?
# TODO check if Limbo does the same thing
```

    (900, 900) (900, 900) (900,)

Note that although the training data (XY_train) is multi-dimensional: $(900,2)$, the covariance is still $(900,900)$. The kernel function takes each pair of 2D vectors and outputs one value, which in this case is squared exponential distance.

```python
# Generate test points on grid (for querying GP predictions)
n_test = 70
X_test = np.linspace(-3, 3, n_test)
Y_test = np.linspace(-3, 3, n_test)
XY_test = xy_combinations(X_test, Y_test)
XY_test.shape # number of points will be 70^2
```

    (4900, 2)

```python
# Find function mean at test points
Lk = np.linalg.solve(L, se_kernel(XY_train, XY_test))
mu = np.dot(Lk.T, np.linalg.solve(L, y))
```

```python
# Plot
X3 = XY_test[:,0].flatten()
Y3 = XY_test[:,1].flatten()

fig = plt.figure()
ax = fig.gca(projection='3d')
# ax.plot_surface(X3, Y3, mu, rstride=1, cstride=1, linewidth=0, antialiased=False, color='g') #cmap=cm.plasma
ax.scatter3D(X3, Y3, mu, color='g', s=1)
ax.scatter3D(XY_plot[:,0].flatten(), XY_plot[:,1].flatten(), Z, color='r', s=0.2) # c=Z, cmap='viridis'
plt.legend(['prediction', 'ground truth'])
```

With $n_\text{samples}=400$, it finds a decent fit. With $n_\text{samples}=900$, the GP seems to be able to fit the function well in the region that training points were provided!

---

```python
# Testing jupter notebook built-in commands
# http://ipython.readthedocs.io/en/stable/interactive/magics.html
%mkdir test
%cd test
```

    mkdir: test: File exists
    /Users/Kazu/Code/notebooks/test

```python
%%time
import time
for _ in range(200):
    time.sleep(0.01) #[s]
```

    CPU times: user 5.05 ms, sys: 4.93 ms, total: 9.99 ms
    Wall time: 2.29 s

```python
#TODO: plot some error metrics across n_samples
#TODO: take cross sections on a slice of x or y dimension to check fit closely
```

---

### Function approximation with neural networks

Related? https://ac.els-cdn.com/S1877050913002093/1-s2.0-S1877050913002093-main.pdf?_tid=03a1388a-089c-11e8-8dc0-00000aab0f26&acdnat=1517632860_5ae101a5d49e99d554ad7a1035da86a9

```python
import tensorflow as tf
from tqdm import tqdm #progress bar
```

```python
xy_min, xy_max = -5, 5
n_samples = 1000
```

```python
def create_data():
    X = np.random.uniform(xy_min, xy_max, n_samples).reshape((n_samples,1))
    Y = np.random.uniform(xy_min, xy_max, n_samples).reshape((n_samples,1))

    XY = np.concatenate((X,Y), axis=1)
    Z = rastrigin(XY.T)

    return XY.astype(np.float32), Z.astype(np.float32)
```

```python
def create_model():
    xy_coords = tf.placeholder(tf.float32, shape=[None, 2], name="xy-input")

    with tf.variable_scope("layer1"):
        hidden = tf.layers.dense(xy_coords, 512, activation=tf.nn.relu)
    with tf.variable_scope("layer2"):
        hidden = tf.layers.dense(hidden, 512, activation=tf.nn.relu)
    with tf.variable_scope("layer3"):
        hidden = tf.layers.dense(hidden, 512, activation=tf.nn.relu)
    with tf.variable_scope("output"):
        z_val_pred = tf.layers.dense(hidden, 1)
        z_val_pred = tf.reshape(z_val_pred, [-1])

    return xy_coords, z_val_pred
```

```python
def create_training(z_val_pred):
    label_ph = tf.placeholder(tf.float32, shape=[None], name="z-labels")

    with tf.variable_scope("loss"):
        loss = tf.losses.mean_squared_error(labels=label_ph, predictions=z_val_pred)
        tf.summary.scalar('loss', loss)

    learning_rate = tf.placeholder(tf.float32, shape=[], name="learning_rate")
    with tf.variable_scope("training"):
        optimizer = tf.train.AdamOptimizer(learning_rate=learning_rate)
        train_op = optimizer.minimize(loss=loss)
        tf.summary.scalar('learning_rate', learning_rate)

    return train_op, loss, label_ph, learning_rate
```

```python
tf.reset_default_graph()
xy_coords, z_val_pred = create_model()
train, loss, labels, learning_rate = create_training(z_val_pred)
summary_op = tf.summary.merge_all()
```

```python
sess = tf.Session()
sess.run(tf.global_variables_initializer())

lr = 1e-3

n_samples_total = 0
for i in tqdm(range(600)):
    xy_train, z_train = create_data()
    n_samples_total += n_samples

    _, cur_loss, cur_summaries = sess.run([train, loss, summary_op],
        feed_dict={ xy_coords: xy_train, labels: z_train, learning_rate: lr})

    if i%50==0:
        print(cur_loss)

    if(cur_loss < 1):
        break
print(cur_loss)
```

```python
query = np.array([1, 1], dtype=np.float32)
sess.run(z_val_pred, feed_dict={xy_coords: [query.flatten()]})
```

```python
# Plot learned function
N = 50
X = np.linspace(xy_min, xy_max, N)
XY_test = xy_combinations(X, X).astype(np.float32)

Z_test = sess.run(z_val_pred, feed_dict={xy_coords: XY_test})
Z_tgt = rastrigin(XY_test.T)
```

```python
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.scatter3D(XY_test[:,0], XY_test[:,1], Z_test, color='g', s=1)
# ax.plot_trisurf(XY_test[:,0], XY_test[:,1], Z_test, color='g')

# ax.scatter3D(XY_test[:,0], XY_test[:,1], Z_tgt, color='r', s=1)
# ax.plot_trisurf(XY_test[:,0], XY_test[:,1], Z_tgt, color='r')
```

Neural network can't seem to learn the wavy-ness of a Rastrigin function. This is partly because of the mse loss.. Is this even a fair comparison?

### Bayesian optimization

Bayesian optimization is an optimization strategy that treats the objective function as a random function (here, a GP) and places a prior over it. It iteratively picks a sample point (based on the GP) and updates its model, in an effort to find the minimum point. It can be used for black-box functions (don't need to know structure, just have to be able to query at specific points), and doesn't require derivatives. Generally good for functions that are expensive to evaluate, since the uncertainty encoded in a GP model allows for balancing explorating-exploitation behavior.

```
initialize D as empty prior
for t = 1, 2, ... do
    Find x_t by combining attributes of posterior distribution
    in a utility function u and maximizing:
    x_t = argmax_x u(x, D)

    Sample the objective function:
    y_t = f(x_t) + z_t

    Update the GP with new data (x_t, y_t_)
```

Utility function `u(x, D)` is also called acquisition function. We want to choose the next point $x$ where the mean is high (exploitation) and the variance is high (exploration, and assuming we're maximizing). Balance this tradeoff in a function:
$$\mu(x) + \kappa \sigma(x)$$

Other acquisition functions: Bayesian expected losses, upper confidence bounds (UCB), Thompson sampling

For many applications, the acquisition function is much cheaper to evaluate than the real function. By using the acquisition function as a cheap proxy, and optimizing that, we can try to optimize the true objective function.

Try Bayesian optimization for:

* function minimization
* multi-arm bandit?
* Also used in hyperparameter optimization for ML algorithms

References:

* http://www.cs.ubc.ca/~nando/540-2013/lectures/l7.pdf
* https://en.wikipedia.org/wiki/Bayesian_optimization
* https://github.com/fmfn/BayesianOptimization
* https://www.iro.umontreal.ca/~bengioy/cifar/NCAP2014-summerschool/slides/Ryan_adams_140814_bayesopt_ncap.pdf

Example of Bayesian Optimization on Rastrigin function is shown below.

Two acquisition functions. Note that we subtract term from $\mu$, because we're trying to minimize

* Gaussian Lower Confidence Bound (GLCB): $\mu(x) - \kappa \sigma(x)$
* UCB: $\mu(x) - \sqrt{\alpha \log(t+1) / 2T(i)}$, where $t$ is timestep and $T(i)$ is number of visits for a state $i$
    * Typically used in discrete bandit case; trying it here
    * See [here](http://banditalgs.com/2016/09/18/the-upper-confidence-bound-algorithm/) for theory behind this

Other options for acquisition functions: Gaussian Probability of Improvement, Gaussian Expected Improvement. See [here](http://efavdb.com/gaussian-processes/)

BO finds global minimum within a few queries. TODO compare to other people's implementations

Note that in this example implementation, I just evaluate the acquisition function over a grid and pick the best one. In practice, people usually run a global optimization algorithm (e.g. CMA-ES) on the acquisition function (i.e. an inner loop optimization): https://stats.stackexchange.com/questions/194166/practical-implementation-detail-of-bayesian-optimization

```python
import scipy

#Test points on grid (for querying acquisition function)
n_probe = 51 # make sure this is odd [0, 0] is included
X_probe = np.linspace(-3, 3, n_probe)
Y_probe = np.linspace(-3, 3, n_probe)
XY_probe = xy_combinations(X_probe, Y_probe)
XY_probe.shape # number of points will be 70^2

s = 5e-5    # scale parameter / noise variance prior

# For UCB. Initialize with small number, to avoid divide by zero
# This still encourages UCB to query unvisited states first, but it's a hacky way to do this
# TODO how is it supposed to be done?
visit_counts = np.ones(XY_probe.shape[0]) * 1e-3

for t in range(300):
    if t==0:
        query_point = np.random.uniform(-3, 3, 2).reshape((1,2)) # initial guess
        XY_train = query_point
        y_train = rastrigin(query_point.T)
    else:
        y_query = rastrigin(query_point.T) # "Expensive" evaluation
        XY_train = np.concatenate((XY_train, query_point), 0)
        y_train = np.concatenate((y_train, y_query), 0)

    ### Fit a GP to the data we have so far. # TODO is there a way to do this incrementally?
    K = se_kernel(XY_train, XY_train)
    K_ = se_kernel(XY_probe, XY_probe)
    L = np.linalg.cholesky(K + s*np.eye(XY_train.shape[0]))
    Lk = np.linalg.solve(L, se_kernel(XY_train, XY_probe))
    s2 = np.diag(K_) - np.sum(Lk**2, axis=0)

    ### "Cheap" GP evaluation
    mu = np.dot(Lk.T, np.linalg.solve(L, y_train))
    sigma = np.sqrt(s2)

    ### Find acquisition function values at test points

    # Function 1: GLCB
    kappa = 0.5
    acq = mu - kappa*sigma

    # Function 2: UCB
#     alpha = 1.0
#     acq = mu - np.sqrt((alpha*np.log(t+1)) / (2*visit_counts))

    # Function 3: Gaussian Probability of Improvement
#     acq = -scipy.stats.norm.cdf(np.min(y_train), mu, sigma) # TODO make sure this works

    # Function 4: Expected improvement
    # TODO implement this "Practical Bayesian Optimization of Machine Learning Algorithms" claims this is empirically better behaved

    ### Pick next point to evaluate

    # Method 1: greedily pick min. prone to get stuck in local minima?
#     sample_idx = np.argmin(acq)

    # Method 2: sample over softmax (Boltzmann?) distribution
    # -acq because the lower the acquisition function (some measure of expected value) is, we want higher prob
    sample_probs = np.exp(-acq - scipy.special.logsumexp(-acq))
    sample_idx = np.random.multinomial(1, sample_probs).argmax()

    visit_counts[sample_idx] += 1

    query_point = XY_probe[sample_idx, :].reshape(1, 2)
    if t%20==0:
        print("Iter {}, best so far: {:.2f} at {}".format(t, y_train.min(), XY_train[np.argmin(y_train), :]))
#         print("Acquistion function at [0, 0]: {:.3f}. min acq: {:.3f}".format(acq[XY_probe.shape[0]//2], np.min(acq)))

print("Best: {:.2f} at {}".format(y_train.min(), XY_train[np.argmin(y_train), :]))
print("Target: {:.2f} at [0, 0]".format(rastrigin(np.array([0, 0]))))
```

    Iter 0, best so far: 3.69 at [0.7992 0.727 ]
    Iter 20, best so far: 3.16 at [-0.84  1.92]
    Iter 40, best so far: 0.64 at [1.08 0.84]
    Iter 60, best so far: 0.11 at [0.96 0.84]
    Iter 80, best so far: -1.46 at [ 0.   -1.08]
    Iter 100, best so far: -1.98 at [ 0.   -0.96]
    Iter 120, best so far: -1.98 at [ 0.   -0.96]
    Iter 140, best so far: -1.98 at [ 0.   -0.96]
    Iter 160, best so far: -1.98 at [ 0.   -0.96]
    Iter 180, best so far: -1.98 at [ 0.   -0.96]
    Iter 200, best so far: -1.98 at [ 0.   -0.96]
    Iter 220, best so far: -1.98 at [ 0.   -0.96]
    Iter 240, best so far: -1.98 at [ 0.   -0.96]
    Iter 260, best so far: -1.98 at [ 0.   -0.96]
    Iter 280, best so far: -1.98 at [ 0.   -0.96]
    Best: -1.98 at [ 0.   -0.96]
    Target: -3.00 at [0, 0]

```python
# Visualize what our model has learned so far about the underlying function
fig = plt.figure()
ax = fig.gca(projection='3d')

Z_probe = mu
# ax.scatter3D(XY_test[:,0], XY_test[:,1], Z_test, color='g', s=1)
ax.plot_trisurf(XY_probe[:,0], XY_probe[:,1], Z_probe, cmap='viridis')
```

On why BOs don't do well in high-dimensional spaces, from "High Dimensional Bayesian Optimization with Elastic Gaussian Process":

> However, it has a serious limitation when the parameter space is high-dimensional as Bayesian optimization crucially depends on solving a global optimization of a surrogate utility function in the same sized dimensions. The surrogate utility function, known commonly as acquisition function is a continuous function but can be extremely sharp at high dimension - having only a few peaks marooned in a large terrain of almost flat surface.

Wang, Zoghi, Hutter, Matheson & de Freitas, IJCAI 2013 suggest a way to deal with many dimensions. Approach: choose a random embedding of the input dimensions and then optimize in the subspace

**Would this work for high-dim robot control? Or should we make some hand-made features to embed from high-dim to low-dim, like Akshara Rai paper? **

```python

```


text:
```python
import numpy as np
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import cv2
import math
```

```python
%%HTML
<style> code {background-color : orange !important;} </style>
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
```

<style> code {background-color : orange !important;} </style>
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))

```python
# Load our image
binary_warped = mpimg.imread('test_images/warped_example.jpg')

def find_lane_pixels(binary_warped):
    # Take a histogram of the bottom half of the image
    histogram = np.sum(binary_warped[binary_warped.shape[0]//2:,:], axis=0)

    # Create an output image to draw on and visualize the result
    out_img = np.dstack((binary_warped, binary_warped, binary_warped))

    # Find the peak of the left and right halves of the histogram
    # These will be the starting point for the left and right lines
    midpoint = np.int(histogram.shape[0]//2)
    leftx_base = np.argmax(histogram[:midpoint])
    rightx_base = np.argmax(histogram[midpoint:]) + midpoint

    # HYPERPARAMETERS
    # Choose the number of sliding windows
    nwindows = 9
    # Set the width of the windows +/- margin
    margin = 100
    # Set minimum number of pixels found to recenter window
    minpix = 50

    # Set height of windows - based on nwindows above and image shape
    window_height = np.int(binary_warped.shape[0]//nwindows)

    # Identify the x and y positions of all nonzero pixels in the image
    nonzero = binary_warped.nonzero()
    nonzeroy = np.array(nonzero[0])
    nonzerox = np.array(nonzero[1])

    #print("nonzeroxy", nonzero)
    #print("nonzerox", nonzerox)
    #print("nonzeroy", nonzeroy)
    # Current positions to be updated later for each window in nwindows
    leftx_current = leftx_base
    rightx_current = rightx_base

    # Create empty lists to receive left and right lane pixel indices
    left_lane_inds = []
    right_lane_inds = []

    # Step through the windows one by one
    for window in range(nwindows):
        #print("======================= processing window: ", window)
        # Identify window boundaries in x and y (and right and left)
        win_y_low = binary_warped.shape[0] - (window+1)*window_height
        win_y_high = binary_warped.shape[0] - window*window_height

        ### TO-DO: Find the four below boundaries of the window ###
        win_xleft_low = leftx_current -  margin# Update this
        win_xleft_high = leftx_current +  margin # Update this
        win_xright_low = rightx_current -  margin  # Update this
        win_xright_high = rightx_current +  margin  # Update this

        # Draw the windows on the visualization image
        lineWidth = 3
        color = (0,255,0)
        cv2.rectangle(out_img, (win_xleft_low, win_y_low), (win_xleft_high, win_y_high), color, lineWidth)
        cv2.rectangle(out_img, (win_xright_low, win_y_low), (win_xright_high, win_y_high), color, lineWidth)

        leftWindow = binary_warped[win_y_low : win_y_high, win_xleft_low : win_xleft_high]
        rightWindow = binary_warped[win_y_low : win_y_high, win_xright_low : win_xright_high]

        ### TO-DO: Identify the nonzero pixels in x and y within the window ###
        good_left_inds = np.transpose(np.nonzero(leftWindow))
        good_right_inds = np.transpose(np.nonzero(rightWindow))
        #good_left_inds = np.nonzero(leftWindow)
        #good_right_inds = np.nonzero(rightWindow)

        # Identify the nonzero pixels in x and y within the window #
        good_left_inds = ((nonzeroy >= win_y_low) & (nonzeroy < win_y_high) & (nonzerox >= win_xleft_low) &  (nonzerox < win_xleft_high)).nonzero()[0]
        good_right_inds = ((nonzeroy >= win_y_low) & (nonzeroy < win_y_high) & (nonzerox >= win_xright_low) &  (nonzerox < win_xright_high)).nonzero()[0]

        # Append these indices to the lists
        left_lane_inds.append(good_left_inds)
        right_lane_inds.append(good_right_inds)

        nonzeroLeft = leftWindow.nonzero()
        nonzeroxLeft = nonzeroLeft[1]

        nonzeroRight = rightWindow.nonzero()
        nonzeroxRight = nonzeroRight[1]

        if (len(nonzeroxLeft) > minpix):
            leftGoodPixelAvgX = int(np.mean(nonzeroxLeft))
        else:
            leftGoodPixelAvgX = margin

        if (len(nonzeroxRight) > minpix):
            rightGoodPixelAvgX = int(np.mean(nonzeroxRight))
        else:
            rightGoodPixelAvgX = margin

        #print("leftGoodPixelAvgX= ", leftGoodPixelAvgX)
        #print("rightGoodPixelAvgX= ", rightGoodPixelAvgX)

        left_peak = leftGoodPixelAvgX + win_xleft_low
        right_peak = rightGoodPixelAvgX + win_xright_low
        #print("peaksLR= ({}, {})".format(left_peak, right_peak))

        colorLine = (255, 130, 0)
        Pt0 = (int(left_peak), win_y_low)
        Pt1 = (int(left_peak), win_y_high)
        cv2.line(out_img, Pt0, Pt1, colorLine, lineWidth)

        Pt0 = (int(right_peak), win_y_low)
        Pt1 = (int(right_peak), win_y_high)
        cv2.line(out_img, Pt0, Pt1, colorLine, lineWidth)

        leftx_current = left_peak
        rightx_current = right_peak

    # Concatenate the arrays of indices (previously was a list of lists of pixels)
    try:
        left_lane_inds = np.concatenate(left_lane_inds)
        right_lane_inds = np.concatenate(right_lane_inds)
    except ValueError:
        # Avoids an error if the above is not implemented fully
        pass

    # Extract left and right line pixel positions
    leftx = nonzerox[left_lane_inds]
    lefty = nonzeroy[left_lane_inds]
    rightx = nonzerox[right_lane_inds]
    righty = nonzeroy[right_lane_inds]

    return leftx, lefty, rightx, righty, out_img

#%matplotlib qt4
#%matplotlib inline
#%matplotlib notebook
#leftx, lefty, rightx, righty, out_img = find_lane_pixels(binary_warped)
#plt.figure(figsize=(12,8))
#plt.imshow(out_img)
#plt.tight_layout()

```

```python
#%matplotlib qt4
#%matplotlib notebook
#%matplotlib inline

def fit_polynomial(binary_warped):
    # Find our lane pixels first
    leftx, lefty, rightx, righty, out_img = find_lane_pixels(binary_warped)

    ### TO-DO: Fit a second order polynomial to each using `np.polyfit` ###
    left_fit = np.polyfit(lefty, leftx, 2)
    right_fit = np.polyfit(righty, rightx, 2)

    # Generate x and y values for plotting
    ploty = np.linspace(0, binary_warped.shape[0]-1, binary_warped.shape[0] )
    try:
        left_fitx = left_fit[0]*ploty**2 + left_fit[1]*ploty + left_fit[2]
        right_fitx = right_fit[0]*ploty**2 + right_fit[1]*ploty + right_fit[2]
    except TypeError:
        # Avoids an error if `left` and `right_fit` are still none or incorrect
        print('The function failed to fit a line!')
        left_fitx = 1*ploty**2 + 1*ploty
        right_fitx = 1*ploty**2 + 1*ploty

    ## Visualization ##
    # Colors in the left and right lane regions
    out_img[lefty, leftx] = [255, 0, 0]
    out_img[righty, rightx] = [0, 0, 255]

    plt.imshow(out_img)
    # Plots the left and right polynomials on the lane lines
    plt.plot(left_fitx, ploty, color='yellow', linewidth=3)
    plt.plot(right_fitx, ploty, color='yellow', linewidth=3)
    plt.show()
    return out_img

out_img = fit_polynomial(binary_warped)

```

```python
# Polynomial fit values from the previous frame
# Make sure to grab the actual values from the previous step in your project!
left_fit = np.array([ 2.13935315e-04, -3.77507980e-01,  4.76902175e+02])
right_fit = np.array([4.17622148e-04, -4.93848953e-01,  1.11806170e+03])

def fit_poly(img_shape, leftx, lefty, rightx, righty):
     ### TO-DO: Fit a second order polynomial to each with np.polyfit() ###
    left_fit = np.polyfit(lefty, leftx, 2)
    right_fit = np.polyfit(righty, rightx, 2)

    # Generate x and y values for plotting
    ploty = np.linspace(0, img_shape[0]-1, img_shape[0])

    ### TO-DO: Calc both polynomials using ploty, left_fit and right_fit ###
    left_fitx = left_fit[0]*ploty**2 + left_fit[1]*ploty + left_fit[2]
    right_fitx = right_fit[0]*ploty**2 + right_fit[1]*ploty + right_fit[2]

    return left_fitx, right_fitx, ploty
```

```python
%matplotlib inline
def search_around_poly(binary_warped):
    # HYPERPARAMETER
    # Choose the width of the margin around the previous polynomial to search
    # The quiz grader expects 100 here, but feel free to tune on your own!
    margin = 100

    # Grab activated pixels
    nonzero = binary_warped.nonzero()
    nonzeroy = np.array(nonzero[0])
    nonzerox = np.array(nonzero[1])

    # Set the area of search based on activated x-values ###
    left_lane_inds = ((nonzerox > (left_fit[0]*(nonzeroy**2) + left_fit[1]*nonzeroy + left_fit[2] - margin))
                    & (nonzerox < (left_fit[0]*(nonzeroy**2) + left_fit[1]*nonzeroy + left_fit[2] + margin)))

    right_lane_inds = ((nonzerox > (right_fit[0]*(nonzeroy**2) + right_fit[1]*nonzeroy +right_fit[2] - margin))
                    &  (nonzerox < (right_fit[0]*(nonzeroy**2) +  right_fit[1]*nonzeroy + right_fit[2] + margin)))

    # Again, extract left and right line pixel positions
    leftx = nonzerox[left_lane_inds]
    lefty = nonzeroy[left_lane_inds]
    rightx = nonzerox[right_lane_inds]
    righty = nonzeroy[right_lane_inds]

    # Fit new polynomials
    left_fitx, right_fitx, ploty = fit_poly(binary_warped.shape, leftx, lefty, rightx, righty)

    ## Visualization ##
    # Create an image to draw on and an image to show the selection window
    out_img = np.dstack((binary_warped, binary_warped, binary_warped))*255
    window_img = np.zeros_like(out_img)

    # Color in left and right line pixels
    out_img[nonzeroy[left_lane_inds], nonzerox[left_lane_inds]] = [255, 0, 0]
    out_img[nonzeroy[right_lane_inds], nonzerox[right_lane_inds]] = [0, 0, 255]

    # Generate a polygon to illustrate the search window area
    # And recast the x and y points into usable format for cv2.fillPoly()
    left_line_window1 = np.array([np.transpose(np.vstack([left_fitx-margin, ploty]))])
    left_line_window2 = np.array([np.flipud(np.transpose(np.vstack([left_fitx+margin, ploty])))])
    left_line_pts = np.hstack((left_line_window1, left_line_window2))

    right_line_window1 = np.array([np.transpose(np.vstack([right_fitx-margin, ploty]))])
    right_line_window2 = np.array([np.flipud(np.transpose(np.vstack([right_fitx+margin, ploty])))])
    right_line_pts = np.hstack((right_line_window1, right_line_window2))

    # Draw the lane onto the warped blank image
    cv2.fillPoly(window_img, np.int_([left_line_pts]), (0,255, 0))
    cv2.fillPoly(window_img, np.int_([right_line_pts]), (0,255, 0))
    result = cv2.addWeighted(out_img, 1, window_img, 0.3, 0)

    # Plot the polynomial lines onto the image
    plt.plot(left_fitx, ploty, color='yellow')
    plt.plot(right_fitx, ploty, color='yellow')
    ## End visualization steps ##

    return result

# Run image through the pipeline
# Note that in your project, you'll also want to feed in the previous fits
result = search_around_poly(binary_warped)

# View your output
#plt.imshow(result)
```

```python

# Read in a thresholded image
warped = binary_warped
# window settings
window_width = 50
window_height = 80 # Break image into 9 vertical layers since image height is 720
margin = 100 # How much to slide left and right for searching

def window_mask(width, height, img_ref, center,level):
    output = np.zeros_like(img_ref)
    output[int(img_ref.shape[0]-(level+1)*height):int(img_ref.shape[0]-level*height),max(0,int(center-width/2)):min(int(center+width/2),img_ref.shape[1])] = 1
    return output

def find_window_centroids(image, window_width, window_height, margin):

    window_centroids = [] # Store the (left,right) window centroid positions per level
    window = np.ones(window_width) # Create our window template that we will use for convolutions

    # First find the two starting positions for the left and right lane by using np.sum to get the vertical image slice
    # and then np.convolve the vertical image slice with the window template

    # Sum quarter bottom of image to get slice, could use a different ratio
    l_sum = np.sum(image[int(3*image.shape[0]/4):,:int(image.shape[1]/2)], axis=0)
    l_center = np.argmax(np.convolve(window,l_sum))-window_width/2
    r_sum = np.sum(image[int(3*image.shape[0]/4):,int(image.shape[1]/2):], axis=0)
    r_center = np.argmax(np.convolve(window,r_sum))-window_width/2+int(image.shape[1]/2)

    # Add what we found for the first layer
    window_centroids.append((l_center,r_center))

    # Go through each layer looking for max pixel locations
    for level in range(1,(int)(image.shape[0]/window_height)):
	    # convolve the window into the vertical slice of the image
	    image_layer = np.sum(image[int(image.shape[0]-(level+1)*window_height):int(image.shape[0]-level*window_height),:], axis=0)
	    conv_signal = np.convolve(window, image_layer)
	    # Find the best left centroid by using past left center as a reference
	    # Use window_width/2 as offset because convolution signal reference is at right side of window, not center of window
	    offset = window_width/2
	    l_min_index = int(max(l_center+offset-margin,0))
	    l_max_index = int(min(l_center+offset+margin,image.shape[1]))
	    l_center = np.argmax(conv_signal[l_min_index:l_max_index])+l_min_index-offset
	    # Find the best right centroid by using past right center as a reference
	    r_min_index = int(max(r_center+offset-margin,0))
	    r_max_index = int(min(r_center+offset+margin,image.shape[1]))
	    r_center = np.argmax(conv_signal[r_min_index:r_max_index])+r_min_index-offset
	    # Add what we found for that layer
	    window_centroids.append((l_center,r_center))

    return window_centroids

def DrawWindowCentroids(window_centroids):
    # If we found any window centers
    if len(window_centroids) > 0:

        # Points used to draw all the left and right windows
        l_points = np.zeros_like(warped)
        r_points = np.zeros_like(warped)

        # Go through each level and draw the windows
        for level in range(0,len(window_centroids)):
            # Window_mask is a function to draw window areas
            l_mask = window_mask(window_width,window_height,warped,window_centroids[level][0],level)
            r_mask = window_mask(window_width,window_height,warped,window_centroids[level][1],level)
            # Add graphic points from window mask here to total pixels found
            l_points[(l_points == 255) | ((l_mask == 1) ) ] = 255
            r_points[(r_points == 255) | ((r_mask == 1) ) ] = 255

        # Draw the results
        template = np.array(r_points+l_points,np.uint8) # add both left and right window pixels together
        zero_channel = np.zeros_like(template) # create a zero color channel
        template = np.array(cv2.merge((zero_channel,template,zero_channel)),np.uint8) # make window pixels green
        warpage= np.dstack((warped, warped, warped))*255 # making the original road pixels 3 color channels
        output = cv2.addWeighted(warpage, 1, template, 0.5, 0.0) # overlay the orignal road image with window results

    # If no window centers found, just display orginal road image
    else:
        output = np.array(cv2.merge((warped,warped,warped)),np.uint8)

window_centroids = find_window_centroids(warped, window_width, window_height, margin)
DrawWindowCentroids(window_centroids)

# Display the final results
plt.imshow(output)
plt.title('window fitting results')
plt.show()

```

```python
def GererateFakeLanePixels():
    # Generate some fake data to represent lane-line pixels
    ploty = np.linspace(0, 719, num=720)# to cover same y-range as image
    quadratic_coeff = 3e-4 # arbitrary quadratic coefficient
    # For each y position generate random x position within +/-50 pix
    # of the line base position in each case (x=200 for left, and x=900 for right)
    leftx = np.array([200 + (y**2)*quadratic_coeff + np.random.randint(-50, high=51)
                                  for y in ploty])
    rightx = np.array([900 + (y**2)*quadratic_coeff + np.random.randint(-50, high=51)
                                    for y in ploty])

    leftx = leftx[::-1]  # Reverse to match top-to-bottom in y
    rightx = rightx[::-1]  # Reverse to match top-to-bottom in y

    # Fit a second order polynomial to pixel positions in each fake lane line
    left_fit = np.polyfit(ploty, leftx, 2)
    left_fitx = left_fit[0]*ploty**2 + left_fit[1]*ploty + left_fit[2]
    right_fit = np.polyfit(ploty, rightx, 2)
    right_fitx = right_fit[0]*ploty**2 + right_fit[1]*ploty + right_fit[2]

    # Plot up the fake data
    mark_size = 3
    plt.plot(leftx, ploty, 'o', color='red', markersize=mark_size)
    plt.plot(rightx, ploty, 'o', color='blue', markersize=mark_size)
    plt.xlim(0, 1280)
    plt.ylim(0, 720)
    plt.plot(left_fitx, ploty, color='green', linewidth=3)
    plt.plot(right_fitx, ploty, color='green', linewidth=3)
    plt.gca().invert_yaxis() # to visualize as we do the images

GererateFakeLanePixels()
```

### Formula for Radius for curvature
See (http://www.intmath.com/applications-differentiation/8-radius-curvature.php)

$$
\Large
\begin{align}
R_{curve} &=  \frac{(1 + (2Ay+B)^2 )^{3/2}}{|2A|}
\end{align}
$$

```python
# See http://www.intmath.com/applications-differentiation/8-radius-curvature.php
def CalcRadiusOfCurvatureParabola(funcCoeffients, pixY, metersPerPixX=1, metersPerPixY=1):
    # Convert parabola from pixels to meters
    coA = funcCoeffients[0] * metersPerPixX / (metersPerPixY**2)
    coB = funcCoeffients[1] * metersPerPixY
    #coC = funcCoeffients[2] * mPerPixY # Not used

    mY = pixY * metersPerPixY

    radiusNumerator = pow((1 + (2 * coA * mY + coB)**2), 3/2)
    radiusDenominator = abs(2*coA)
    radius = radiusNumerator/radiusDenominator
    return(radius)

```

```python

def generate_data():
    '''
    Generates fake data to use for calculating lane curvature.
    In your own project, you'll ignore this function and instead
    feed in the output of your lane detection algorithm to
    the lane curvature calculation.
    '''
    # Set random seed number so results are consistent for grader
    # Comment this out if you'd like to see results on different random data!
    np.random.seed(0)

    # Generate some fake data to represent lane-line pixels
    ploty = np.linspace(0, 719, num=720)# to cover same y-range as image
    quadratic_coeff = 3e-4 # arbitrary quadratic coefficient

    # For each y position generate random x position within +/-50 pix
    # of the line base position in each case (x=200 for left, and x=900 for right)
    leftx = np.array([200 + (y**2)*quadratic_coeff + np.random.randint(-50, high=51) for y in ploty])
    rightx = np.array([900 + (y**2)*quadratic_coeff + np.random.randint(-50, high=51)  for y in ploty])

    leftx = leftx[::-1]  # Reverse to match top-to-bottom in y
    rightx = rightx[::-1]  # Reverse to match top-to-bottom in y

    # Fit a second order polynomial to pixel positions in each fake lane line
    left_fit = np.polyfit(ploty, leftx, 2)
    right_fit = np.polyfit(ploty, rightx, 2)

    return ploty, left_fit, right_fit

def measure_curvature_pixels():
    '''
    Calculates the curvature of polynomial functions in pixels.
    '''
    # Start by generating our fake example data
    # Make sure to feed in your real data instead in your project!
    ploty, left_fit, right_fit = generate_data()

    # Define y-value where we want radius of curvature
    # We'll choose the maximum y-value, corresponding to the bottom of the image
    y_eval = np.max(ploty) # 719.0

    metersPerPixX = 3.7/700 # 0.005285714285714286 meters per pixel in x dimension 189.189 p/m
    metersPerPixY = 30/720 # 0.041666666666666664 meters per pixel in y dimension 24 p/m

    left_curverad = CalcRadiusOfCurvatureParabola(left_fit, y_eval, metersPerPixX, metersPerPixY)
    right_curverad = CalcRadiusOfCurvatureParabola(right_fit, y_eval, metersPerPixX, metersPerPixY)

    print("left_fit", left_fit)
    print("right_fit", right_fit)
    return left_curverad, right_curverad

# Calculate the radius of curvature in pixels for both lane lines
#left_curverad, right_curverad = measure_curvature_pixels()
left_curverad, right_curverad = measure_curvature_pixels()

# If using the default `generate_data` function with given seed number
# ExpPixels= (1625.06 and 1976.30
# ExpMeters= (533.75 and 648.16)
print(left_curverad, right_curverad)

```

    left_fit [  3.07683280e-04  -4.44713275e-01   3.59067572e+02]
    right_fit [  2.53380891e-04  -3.96103079e-01   1.04908806e+03]
    534.884701658 649.000506718

```python

```

max_stars_repo_name: cielsys/CarND-Advanced-Lane-Lines

max_issues_repo_name: cielsys/CarND-Advanced-Lane-Lines

max_forks_repo_name: cielsys/CarND-Advanced-Lane-Lines


text:
```python
import theme
theme.load_style()
```

# Lesson 12: Numerical Integration

This lecture by Tim Fuller is licensed under the
Creative Commons Attribution 4.0 International License.  All code examples are also licensed under the [MIT license](http://opensource.org/licenses/MIT).

# Topics

- [Numerical Integration](#num_int)
- [Gauss Legendre Quadrature](#gl_quad)
- [Natural Coordinates](#nat_coords)
  - [Approximation of Geometry](#appx_geom)
  - [\*Parametric Forms](#parametric_forms)
- [Examples](#examples)
  - [Example 1: Element Stiffness Integration](#ex-1)
  - [Example 2: Reduced Integration](#ex-2)

#   Numerical Integration[](#top)

Consider the element stiffness

$$
K_{ij}^{(e)} = \int_{x_1^{(e)}}^{x_2^{(e)}} a^{(e)} \phi'_i \phi'_j dx
$$

In the event that $a_e$ is a complicated function of $x$, the exact evaluation of the stiffness may not be possible.  In such cases, numerical integration is necessary.  Likewise, with other integrals encountered in the FEM.

Numerical evaluation of integrals is called numerical quadrature and typically involves approximating the integrand by a polynomial of sufficient degree.

For example, consider the integral

$$
J[f] = \int_{x_A}^{x_B} f(x) dx
$$

We approximate $f(x)$ by

$$
f(x) \approx \sum_i^n f(x_i) \phi_i(x)
$$

where the $x_i$ are the $i^{\rm th}$ point of the interval $[x_A, x_B]$ and $\phi_i(x_i)$ are polynomials of degree $n-1$.  Substitution of the approximate solution for $f$ in to the integral equation gives

$$
\begin{align}
J[f] &= \int_{x_A}^{x_B} f(x) dx \\
 &\approx \sum_{i=1}^n f(x_i) \int_{x_A}^{x_B} \phi_i(x) dx \\
& = f(x_1) \int_{x_A}^{x_B} \phi_1(x) dx + \ldots + f(x_N) \int_{x_A}^{x_B} \phi_n(x) dx
\end{align}
$$

Suppose we choose Lagrange linear interpolation polynomials for the approximation of $f$ with $n=2$.  Then

$$
J[f] \approx f(x_A) \int_{x_A}^{x_B} \frac{x_B - x}{h} dx + f(x_B) \int_{x_A}^{x_B} \frac{x - x_A}{h}  dx = \frac{1}{2}h\left(f(x_A) + f(x_B)\right)
$$

Graphically, this represents the area of a trapezoid used to approximate the area under the function $f(x)$

This is known as the trapzoidal rule.  If we use quadratic Lagrange interpolation of $f(x)$ we obtain

$$
J[f] \approx \frac{1}{6}h \left(f(x_A) + 4 f(x_A + \frac{1}{2}h) + f(x_B)\right)
$$

The previous examples demonstrate the basic idea behind numerical quadrature.  In the trapezoidal type rules, we choose the $x_i$ and the weights shake out naturally.  It can be shown that this form of numerical integration is exact for polynomials of up to degree $n-1$.  The trapezoidal rules also demonstrate a more general form of numerical quadrature:

$$
\begin{align}
J[f] &= \int_{x_A}^{x_B} f(x) w(x) dx \\
J_n[f] &= \int_{x_A}^{x_B}w(x) \sum_{i=1}^n f(x_i) L_i(x)
= \sum_{i=1}^n w_i f(x_i)
\end{align}
$$

with

$$
w_i = \int_{x_A}^{x_B}L_i(x)w(x) dx, \quad
L_i^n(x)=\prod_{\stackrel{j=1}{j \neq i}}^n \frac{x-x_j}{x_j-x_i}
$$

where the $x_i$ are called quadrature points, $w_i$ are quadrature weights, and $L_i(x)$ are $n^{\text{th}}$ order Lagrange polynomials.  The goal of numerical quadrature is to choose the $x_i$ and the $w_i$ such that $J_n$ is exact for polynomials up to degree $2n-1$.

#  Gauss-Legendre Quadrature[](#top)

So, how do we choose the $x_i$ and $w_i$?  That question is the subject of numerical analysis, and the exact details are beyond the scope of this class.  Here, we seek to determine the quadrature points and weights for small $n$ and appeal to reason and mathematical proofs to furnish the results for all $n$.  To that end, we consider only the special case that $w(x)=1$, $a=-1$, $b=1$, and $f(x) \in P_n$, where $P_n$ is the space of polynomials of degree $n$.  This quadrature rule is known as Gauss-Legendre quadrature.  Other quadrature rules are distinguished by their choice of $w$, $a$, and $b$.

So, how do we go about finding the $x_i$ and $w_i$?  Suppose we want to determine the $w_i$ and $x_i$ for $n=2$:

$$
\begin{align}
J_2[f] &= \sum_{i=1}^2 f(x_i) \int_{-1}^{1} L_i(x) = w_1 f(x_1) + w_2 f(x_2)
\end{align}
$$

Evidently, to solve for the unkown $w_i$ and $x_i$, we need 4 equations.  Let us form them by considering the first four polynomials $f=1, x, x^2, x^3$:

$$
\begin{align}
w_1 \cdot 1 + w_2 \cdot 1 = \int_{-1}^1 dx = 2 &\qquad
w_1 \cdot x_1 + w_2 \cdot x_2 = \int_{-1}^1 x dx = 0 \\
w_1 \cdot x_1^2 + w_2 \cdot x_2^2 = \int_{-1}^1 x^2 dx = \frac{2}{3} &\qquad
w_1 \cdot x_1^3 + w_2 \cdot x_2^3 = \int_{-1}^1 x^3 dx = 0
\end{align}
$$

After some algebra and using the methods of undetermined coefficients to solve the nonlinear system of equations, we find that

$$
w_1 = w_2 = 1, \qquad x_1 = -x_2 = -\frac{\sqrt{3}}{3}
$$

and

$$
J_2[f] = f\left(-\frac{\sqrt{3}}{3}\right) + f\left(\frac{\sqrt{3}}{3}\right)
$$

This is called a 2 point Gauss-Legendre quadrature rule.  **The Gauss-Legendre quadrature rule gives exact results if $f(x)$ is a polynomial of degree $2n-1$.**

How can we find the $w_i$ and $x_i$ for higher order Gauss-Legendre quadrature?  We can go through the same process as with $n=2$, but it is slow and tedious.  Instead, using properties of orthogonal polynomials, the $x_i$ and $w_i$ can be derived for many $w_x$ and domains of integration $[a, b]$.  For Gauss-Legendre quadrature, it turns out that the $x_i$ are the roots of the $n^{\text{th}}$ order Legendre polynomial.  The $w_i$ can then be found from knowledge of the $x_i$ and properties of the Legendre polynomials.  The details are beyond the scope of this class, but tabulated values of the $x_i$ and $w_i$ can be found in Table 4.1 of the text.

#  Natural Coordinates[](#top)

To make use of the results from the previous section, the domain of integration for an integral of interest must be $[-1,1]$.  For integration of element functions, this implies that we must perform a change of coordinates from the element domain $[x_1^{(e)},x_2^{(e)}]$ to $[-1,1]$.  The transformed domain $[-1, 1]$ is referred to as the **natural coordinates**.  The actual element domain will be referred to as the **physical coordinates**.  To distinguish natural coordinates from the physical coordinates, the natural coordinates will be labeled $\xi$.

Mathematically, we seek a transformation such that

$$
\xi = \begin{cases}
-1 & \text{for } x = x_1^{(e)} \\
1 & \text{for } x = x_2^{(e)} \\
\end{cases}
$$

Then the transformation between $x$ and $\xi$ can be represented by

$$
x = a + b\xi
$$

which is known as a "stretch" transformation.  We determine the constants $a$ and $b$ be requiring that

$$
\begin{align}
x_1^{(e)} &= a + b(-1) \\
x_2^{(e)} &= a + b(1)
\end{align}
$$

Solving for $a$ and $b$ we get

$$
a = b = \frac{1}{2} h^{(e)}
$$

Then, the transformation becomes

$$
x = x_1^{(e)} + \frac{1}{2}h^{(e)} \left( 1 + \xi \right)
$$

where $x_1^{(e)}$ denotes the left end node of the element and $h^{(e)}$ is the element length.

The value of the natural coordinate always lies in $[-1,1]$ with origin at the element center, as shown:

##  Approximation of Geometry

Recall that the element stiffness and force arrays (for second order boundary value problems) are given by

$$
k_{ij} = \int_{x_A}^{x_B} a\phi'_i\phi'_j + c\phi_i\phi_j dx
$$
$$
f_i = \int_{x_A}^{x_B} q\phi_i + P\delta(x - x_i) dx
$$

These element matrices involve derivatives of the interpolation functions with respect to the physical coordinates $x$.  Thus, we need to find a transformation of the form

$$
x = \chi(\xi)
$$

so that the integrals can be cast in terms of $\xi$ on $[-1, 1]$.  The transformation $\chi$ is assumed to be 1 to 1.  An example of one such transformation is

$$
\chi(\xi) = x_1^{(e)} + \frac{1}{2}h^{(e)}\left(1+\xi\right)
$$

which is the same as our original transformation.  In this transformation $\chi(\xi)$ is a linear function of $\xi$.  It will map a straight line to a straight line.  If $\chi$ were non-linear, it would map a straight line to a curved line.

It is natural to think of approximating the geometry in the same way as we approximated a dependent variable.  In other words, the transformation $x=\chi(\xi)$ can be written as

$$
x = \chi(\xi) = \sum_{i=1}^m x_i^{(e)} \phi_i^{(e)}(\xi)
$$

where $x_i^{(e)}$ is the physical coordinate of the $i^{\rm th}$ node of element $(e)$.  $\phi_i^{(e)}(\xi)$ are interpolating polynomials of degree $m-1$. This equation represents the shape of the element and the $\phi_i^{(e)}(\xi)$ are called the **shape functions**, aften denoted as $\{N(\xi)\}$.

Physically, this mapping maps a geometric shape from $\xi$ to $x$.  In other words, for any given value of $\xi$ we can determine the corresponding value of $x$.  When the element is a straight line, the preceding transformation is the same as that for the linear case.

The transformation is not used to change the actual geometry of the element - the transformation is used to write

$$
\int_{x_1^{(e)}}^{x_2^{(e)}}f(x)dx = \int_{-1}^{1}f(\xi)d\xi
$$

so that Gauss quadrature can be used to approximate the integral.

To complete the transformation, we need $dx$ in terms of $d\xi$.  $dx$ in the global coordinates is related to $d\xi$ by

$$
dx = \frac{\partial x}{\partial \xi} d\xi = J^{(e)} d\xi
$$

where $J^{(e)}$ is called the "Jacobian" of the transformation.  Using our transformation, we have

$$
J^{(e)} = \frac{\partial x}{\partial \xi} =
\frac{\partial }{\partial \xi} \sum_{i=1}^{m} x_i^{(e)}\frac{\partial \phi_i^{(e)}(\xi)}{\partial \xi}
$$

##  \*Parametric Formulations

Recall that our approximate solution is of the form

$$
u_N(x) = \sum_{j=1}^{n} u_j \phi_j(x)
$$

In general, the degree of the transformation $m$ does not have to equal the degree of the approximation $n$.  In other words, $\phi_i(x) \neq \phi_i(\xi)$.

Two elements can be used in the finite element analysis: one for the approximation of the geometry and the other for the interpolation of $u$.  Finite element formulations are classified as

- Subparametric: m < n
- Isoparametric: m = n
- Superparametric: m > n

Examples:

- Subparametric: In the Euler-Bernouli beam element the transverse deflection is approximated by the cubic Hermite polynomials and straight beams are approximated by linear polynomials.
- Isoparametric: The same element is used to approximate the geometry and dependent variable.
- Superparametric: The geometry is approximated by higher order polynomial than the dependent variable.  Rarely used in practice.

#  Examples[](#top)

##  Example 1: Element Stiffness

We now piece together what we have done so far.  Consider the FE stiffness matrix in global coordinates

$$
k_{ij}^{(e)} = \int_{x_A}^{x_B} a_e \frac{d\phi_i^{(e)}}{dx} \frac{d\phi_j^{(e)}}{dx} dx
$$

To transform the the integral from $x$ to $\xi$, we note that

$$
a(x) = a(\chi(\xi)), \quad dx = Jd\xi
$$

$$
\frac{d\phi}{dx} = \frac{d\phi}{d\xi}\frac{d\xi}{dx} = \frac{1}{J}\frac{d\phi}{d\xi}
$$

Therefore, the stiffness can be written as

$$\begin{align}
k_{ij}^{(e)} &= \int_{-1}^{1} a^{(e)}(\chi(\xi))
                \frac{1}{J^{(e)}}\frac{d\phi_i^{(e)}}{d\xi}\frac{1}{J^{(e)}}\frac{d\phi_j^{(e)}}{d\xi}
                J d\xi \\
             &= \int_{-1}^{1} \frac{a^{(e)}(\chi(\xi))}{J^{(e)}} \frac{d\phi_i^{(e)}}{d\xi}\frac{d\phi_j^{(e)}}{d\xi}d\xi \\
             &= \sum_{I=1}^n Q_{ij}^{(e)}\left(\xi_I\right)w_I
\end{align}
$$

where

$$
Q_{ij}^{(e)} = \frac{a^{(e)}}{J^{(e)}} \frac{d\phi_i^{(e)}}{d\xi}\frac{d\phi_j^{(e)}}{d\xi}
$$

For linear interpolation elements and an isoparametric mapping, the shape functions are

$$
\{\phi(\xi)\}^{(e)}
= \frac{1}{2}
\begin{Bmatrix}
1 - \xi \\
1 + \xi
\end{Bmatrix}
$$

and

$$
J = \sum_{i=1}^m x_i^{(e)} \frac{d\phi_i^{(e)}}{d\xi} = \frac{1}{2}\left(x_1(-1) + x_2(1)\right) = \frac{1}{2}h^{(e)}
$$

For elementwise constant values of $a$ and $h$

$$\begin{align}
k_{ij}^{(e)} &= Q_{ij}^{(e)}(\xi_1) w_1 + Q_{ij}^{(e)}(\xi_2) w_2 \\
       &= \frac{2a^{(e)}}{h^{(e)}}\left(2\frac{d\phi_i^{(e)}}{d\xi}\frac{d\phi_j^{(e)}}{d\xi}\right)
\end{align}
$$

from which

$$ k_{11}^{(e)} = \frac{2a^{(e)}}{h^{(e)}} \left(2(-.5)(-.5)\right) = \frac{a}{h}^{(e)} $$
$$ k_{12}^{(e)} = \frac{2a^{(e)}}{h^{(e)}} \left(2(-.5)(.5)\right) = -\frac{a}{h}^{(e)} $$
$$ k_{21}^{(e)} = \frac{2a^{(e)}}{h^{(e)}} \left(2(.5)(-.5)\right) = -\frac{a}{h}^{(e)} $$
$$ k_{22}^{(e)} = \frac{2a^{(e)}}{h^{(e)}} \left(2(.5)(.5)\right) = \frac{a}{h}^{(e)} $$

and

$$
[k]^{(e)} = \frac{a^{(e)}}{h^{(e)}} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}
$$

As noted before, the Jacobian matrix will be the same - $J^{(e)} = h^{(e)} / 2$ - when the element is straight, even when the coordinate transformation is not linear.  However, when the element is curved, the Jacobian is a function of $\xi$ for transformations other than linear.

##  Example 2: Reduced Integration

It is often beneficial to reduce the number of integration points that are evaluated.  In the case of the element stiffness with constant properties, single point Gauss-Legendre quadrature can be used to find $k_{ij}^{(e)}$ exactly:

$$
\begin{align}
k_{ij} &= Q_{ij}^{(e)}(\xi_1)w_1 \\
&= \frac{2a^{(e)}}{h^{(e)}}\left(\frac{d\phi_i^{(e)}}{d\xi}\frac{d\phi_j^{(e)}}{d\xi}\right) 2 \\
& = \frac{a^{(e)}}{h^{(e)}} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}
\end{align}
$$

where we used the fact that $w_1=2$ for single point integration.

It turns out that in higher dimensions, using reduced integration can have a dramatic influence on the finite element solution - both good and bad!

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: pinkieli/fem-with-python

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: pinkieli/fem-with-python

max_forks_repo_licenses: MIT Unlicense