text:
 **Chapter 2: [Diffraction](CH2_00-Diffraction.ipynb)**

# Bethe Theory

[Download](https://raw.githubusercontent.com/gduscher/MSE672-Introduction-to-TEM//main/Diffraction/CCH2_D01-Bloch.ipynb)

[](
    https://colab.research.google.com/github/gduscher/MSE672-Introduction-to-TEM/blob/main/Diffraction/CH2_D01-Bloch.ipynb)

part of

 **[MSE672:  Introduction to Transmission Electron Microscopy](../_MSE672_Intro_TEM.ipynb)**

by Gerd Duscher, Spring 2021

Microscopy Facilities

Joint Institute of Advanced Materials

Materials Science & Engineering

The University of Tennessee, Knoxville

Background and methods to analysis and quantification of data acquired with transmission electron microscopes

## Import numerical and plotting python packages
Import the python packages that we will use:

Beside the basic numerical (numpy) and plotting (pylab of matplotlib) libraries,

```python
import sys
if 'google.colab' in sys.modules:
    %pylab --no-import-all inline
else:
    %pylab --no-import-all notebook
```

    Populating the interactive namespace from numpy and matplotlib

## Notation

In the following I will introduce this Bethe diffraction theory in Dirac's  bra--ket notation.
\begin{eqnarray*}
|\psi\rangle &=& \mbox{ket} = \psi_1 |1\rangle + \psi_2 |2\rangle\\
&=& \left( \begin{array}{c} \psi_1 \\\psi_2  \end{array} \right)= \mbox{column vector}\\
\langle \phi| &=& \mbox{bra} = \phi^*_1\langle 1|  + \phi^*_2 \langle 2|\\
&=& \left( \phi^*_1 ,\phi^*_2 \right)= \mbox{row vector}
\end{eqnarray*}

Consequently the Dirac bracket is:
\begin{equation*}
\langle \phi|\psi\rangle = \left( \phi^*_1 ,\phi^*_2 \right) \left( \begin{array}{c} \psi_1 \\\psi_2  \end{array} \right) = \phi^*_1 \psi_1 +\phi^*_2 \psi_2
\end{equation*}

Similarly

\begin{equation*}
|\psi\rangle \langle \phi| = \mbox {a matrix}
\end{equation*}
and this will make it easier to follow the equations of many beams in diffraction, which otherwise can become very messy.

Throughout this chapter we will use the notation for the wave vector common in solid state physics and not the one found in many crystallography books ($k=1/\lambda$):
\begin{equation}\label{vector}
k = \frac{2\pi}{\lambda} = \frac{\sqrt{2meE}}{\hbar}
\end{equation}
with:
- $\lambda$: wavelength
- $m=\gamma_0$: relativistically corrected mass of the electron
- $\gamma = 1+\frac{e^2E^2}{2m_0c^2}$: relativistic correction necessary for $E$ above 2kV
- $e$: electron charge
- $E$: acceleration voltage
- $c$: speed of light in vacuum

This is the short version of quantum mechanics. The Axioms of quantum mechanics are:

1. The state of a system is described by its state vector $|\psi\rangle$.

2. An observables are expressed by hermitic operators $A$
3. The mean value of an observable is given by $\langle A \rangle$
4. The time dependence is given by the time-dependent Schrödinger
    equation: $ H|\psi,t\rangle = i\hbar \frac{\partial}{\partial t} |\psi,t\rangle $
5. If you measure $A$ the system changes to $|n\rangle$ if $a_n$ was measured

Axioms 2. and 3. give:

If a system is in the state $|\psi\rangle = \sum_n c_n |n\rangle$

with $c_n = \langle n|\psi\rangle$, where $|n\rangle$ are the Eigen states of $A$, meaning $A|n\rangle = a_n |n\rangle$. The probability to find the value $a_n$ when measuring $A$ is given by $|c_n|^2$.

## Introduction to Bethe Theory

The dynamic theory calculates the probability of a transition from an initial state $|i\rangle$ to a final state $|f\rangle$. In our case, the initial state $|i\rangle$ is the incoming beam which scatters to a final state the Bragg reflection $|f\rangle$.
What we want to know is the transition probability from the initial to the final state:
\begin{equation}
\omega_{i\rightarrow f} =  \langle i |T|f\rangle
\end{equation}
the initial state is the beam with wave vector $\vec{k}_0$ and the end state is a wave-vector of diffracted beam according to Bragg's law $\vec{q} = \vec{g}$.
The above equation now looks like this:
\begin{equation}
\omega_{0\rightarrow g} = \langle 0 |T|g\rangle
\end{equation}

For the stationary problems the states $|\psi,t\rangle$ can be expressed as:
\begin{equation}
|\psi_n,t\rangle = \exp(-E_nt/\hbar) |\psi_n\rangle
\end{equation}
and we get the time-independent Schröinger equation for $|\psi_n\rangle$:
\begin{equation}
H|\psi_n\rangle = E_n |\psi_n \rangle
\end{equation}

For the incoming wave, we get:

\begin{equation}
H_0 = { -\frac{h^2}{8\pi^2} \nabla^2 }
\end{equation}
which expresses the kinetic energy.

Within a crystal the Hamiltonian will change to:
\begin{equation}
H = H_0 + V
\end{equation}
We have  the Schrödinger equation for the incoming wave:
\begin{equation}  \label{IncomingWave}
(H_0 ) |\vec{k}\rangle = E_k |\vec{k}\rangle
\end{equation}

and we want to solve:
\begin{equation}
(H_0 +V) |\psi\rangle = E_k |\psi\rangle
\end{equation}

which we transform with equation \ref{IncomingWave}:

> Effectively, we changed in the integral equation.

\begin{equation*}
|\psi\rangle = |\vec{k}\rangle + \frac{1}{E_k -H_0}V|\psi\rangle
\end{equation*}

### Hamiltonian in Bethe Theory
For our diffraction experiment is often better to use a Hamiltonian that contains the wave vector:
\begin{equation}
H_0 = { -\Delta_\rho -\xi^2 }
\end{equation}

OOps, where is the wavevector $\vec{k}$?

I replaced it by $\xi$ which is the effective deviation from wavevector $\vec{k}$

We also take into account that our electrons are very fast and distort the space, reducing the problem to two dimensions ($\rho = (x,y)$.
The crystal potential V has then to be changed as well to:
\begin{equation}
V = V(\rho,z)
\end{equation}
and the full Hamiltonian is changed to:
\begin{equation}
H = \frac{1}{k_z} (H_0 + V(\rho,z)
\end{equation}

The time- independent wave equation is then:
\begin{equation}
\nabla^2 \psi  + k^2 \psi = 0
\end{equation}
with the plane wave solution:
\begin{equation}
\psi  = \exp(\pm\vec{k}\bullet\vec{r}) = 0
\end{equation}

## Schrödinger  Equation of Bethe Theory

The Bethe theory is based on the (time independent, non relativistic) Schrödinger  equation:

\begin{equation}  \label{Schrodinger}
\Big[ \underbrace{ -\frac{h^2}{8\pi^2} \nabla^2 }_{\mbox{ kinetic energy}} +\underbrace{ \mathcal{V}(\vec{r})}_{\mbox{ pot. energy}} \Big]\, |\psi(\vec{r}) \rangle = \underbrace{\mathcal{E}}_{\mbox{ total energy}} \underbrace{|\psi(\vec{r})\rangle }_{\mbox{ wave function}}
\end{equation}

What does that mean for the TEM?
>
> We have a acceleration voltage (electric field potential) $E$ of 100kV
>
> We have a charge of the electron $q$ with the value $e$.
>
> We have a total Energy $-\mathcal{E} = E\cdot q$ which is just $E$ in the units of $[eV]$.
>
> We have a crystal with the potential $V(\vec{r})$, which we declare positive inside the crystal and zero outside.
>
> We have a potential Energy $\mathcal{V}(\vec{r}) = q\cdot V(\vec{r})$.

Now, that we declared all our variables we can transform the Schrödinger equation \ref{Schrodinger} to:

$$
\nabla^2  |\psi(\vec{r})\rangle = -\frac{8\pi m e}{h^2}\,  [E+V]\, |\psi(\vec{r})\rangle
$$

The left hand part of this equation is the impulse of the electron and the right hand part consists of a total energy part, which is boring and a part which originates from the crystal (interesting!).

## Bloch Waves in Bethe Theory
Well, if the potential is periodic, then the solution (wave function) must be periodic, too.

First we make a substitution in case our wave function is complicated: we define it as a linear combination of other waves. That is a useful trick, which makes the mathematics easier as we'll see in a bit.
\begin{eqnarray} \label{Bloch}
|\psi(\vec{r})\rangle &=& \sum_j |b_j\rangle \\
|b_j\rangle  &=&  = |b(\vec{k}^{(j)}, \vec{r})\rangle
\end{eqnarray}
The $ |b(\vec{k}^{(j)}, \vec{r})\rangle $ are called Bloch waves and are only defined for specific $\vec{k}$--vectors,
because the  $|b(\vec{k}^{(j)}, \vec{r})\rangle$ are plane waves, each traveling in $k^{(j)}$ direction. For the $|\psi(\vec{r})\rangle$ we did not and could not have made any assumption like that.\\

Now, we express the fact that these Bloch waves are indeed plane waves mathematically:

\begin{equation} \label{FourierExpand}
|b^{(j)}(\vec{r})\rangle  = b(\vec{k}^{(j)}, \vec{r})=\mu(\vec{k}^{(j)}, \vec{r})\cdot e^{2\pi i \vec{k}^{(j)} \vec{r}} = \underbrace{\mu^{(j)}(\vec{r})}_{\mbox{ Bloch function}}   e^{2\pi i \vec{k}^{(j)} \vec{r}}
\end{equation}

by dividing it in a plane wave part (the exponential function) and a amplitude part (the Bloch function).
Because of the periodicity which we assume for the solution, we expand the Bloch waves in a into a Fourier series, again (the same as in equation \ref{FourierExpand} for the potential).

\begin{equation}
b^{(j)}(\vec{r}) = \sum_g C_g^{(j)} e^{2\pi i (\vec{k}^{(j)} + \vec{g}) \vec{r}}
\end{equation}
The sum in this equation goes over all excited (aha!) points in the reciprocal lattice, including the incident direction $g_1 = 0$. (  \vec{g} are defined through the Milller indices as (h/a, k/b, l/c), where the (a,b,c) are the real space lattice vectors).

Theoretically, thre are an infinite number of $\vec{g}$ vectors, but only few are allowed and only a few have a small excitation error.

So in practice there are only a few $\vec{g}$ vectors to consider.

## Crystal Potential in Bethe Theory

The crystal potential is periodic and so we also make a Fourier expansion of that potential
$$
V = V(x,y,z) = V(\vec{r}) = \sum_g V_g \exp(2\pi i\vec{g}\cdot\vec{r})
$$

The Fourier component of the crystal potential (in Volts) consists of several atoms $j$ over which we sum:
\begin{eqnarray}
V_g &=& \frac{h^2}{2 \pi m_0 e} \frac{1}{\Omega} \sum_j  f_{e_j}(\vec{g}) \exp(-2\pi i\vec{g}\cdot\vec{r_j})\\
&=& \frac{2 \pi  e a_0}{\Omega} \sum_j  f_{e_j}(\vec{g}) \exp(-2\pi i\vec{g}\cdot\vec{r_j})\\
\end{eqnarray}

where:

- $f_{e_j}$: atomic form (atomic scattering) factor of the $j$th atom
- $e$: charge of electron
- $m_0$: rest mass of electron
- $a_0$: Bohr radius
- $\Omega$: Volume of the unit cell

## Solution of Bethe Theory

Now, so far we haven't done anything, but substitute and expand. Let's put all this into the Schrödinger equation \ref{Schrodinger}
\begin{equation}
4\pi \left[ K^2 - (k_0^{(j)} + g)^2 +\sum_{h \neq 0} U_h e^{2\pi i \vec{h}\vec{r} }  \right] \cdot C_g^{(j)} e^{2\pi i (\vec{k}_o^{(j)} +\vec{g})\cdot \vec{r} }= 0
\end{equation}
This can only be zero, if all coefficients with same exponential function simultaneous become zero; this results in a set of equations, after collecting up terms containing the factor
$e^{2\pi i (\vec{k}_o^{(j)} +\vec{g})\cdot \vec{r} }$:
\begin{equation} \label{setEquat}
\left[ K^2 -(\vec{k}_o^{(j)} +\vec{g})^2 \right] C_g^{(j)} + \sum_{h\neq 0} U_h C_{g-h}^{(j)}=0; \qquad \vec{g}=\vec{g}_1, \vec{g}_2, ..., \vec{g}_n
\end{equation}

I made use of an abbreviation:
\begin{equation}
K=\frac{1}{h}\left[ 2m_0 E (1+\frac{E}{2E_0}) +2m_0 e U_0(1+\frac{E}{E_0}) \right]^{\frac{1}{2}}
\end{equation}
for the wave vector inside the the crystal which are not identical to the magnitude of the wave vectors of the Bloch waves $\vec{k}_g^{(j)}=\vec{k}_o^{(j)} +\vec{g}$.

Please note, that I introduced relativistic corrections (the terms in the round brackets in the equation above), too. It is enough to add this corrections for the energy at this point; it is not necessary to solve the Dirac equation (relativistic  Schrödinger equation).

The set of equations defined in \ref{setEquat} are essential for the  understanding of dynamic diffraction. Let's look at it a little more closely.

We get for each $j$ one equation; and this means we get for each Bloch wave one equation.

The second term in equation \ref{setEquat} (the term with the sum) mixes the Bloch waves ($C_{g-h}^{(j)}$). Effectively, we state that the inner potential $U_h$ mixes the Bloch waves; this is called dynamical coupling.

> **In summation:**
>
> We separated the problem!

## Two Beam Case
We rewrite the matrix expression for the boundary condition in the two beam case:
\begin{equation}
\left(
\begin{array}{*{2}{c}}
C_{0}^{(1)} & C_{0}^{(2)} \\
C_{g}^{(1)} & C_{g}^{(2)} \\
\end{array}
\right) \cdot
\left(
\begin{array}{*{1}{c}}
\epsilon^{(1)} \\
\epsilon^{(2)} \\
\end{array}
\right) =
\left(
\begin{array}{*{1}{c}}
\phi_0^{(0)} \\
\phi_g^{(0)} \\
\end{array}
\right) =
\left(
\begin{array}{*{1}{c}}
1 \\
0\\
\end{array}
\right)
\end{equation}

In the kinematic case, the centers M of the various Ewald spheres (for the various incident directions) lay on a sphere of radius $k=\frac{1}{\lambda}$ around the origin of the reciprocal lattice. At some point the intensity in the diffracted beam will be more intense than in the incident beam, and therefore, we have to treat this scattered beam now in the same way we have treated the incident beam before: we need an Ewald sphere of radius $k$ for this direction/ for this reciprocal lattice point $\vec{g}$. Now we have two Ewald spheres, one around the origin $0$ and one around the reciprocal lattice point $\vec{g}$. The two spheres are not allowed to intersect each other, but a smooth kind of complicated surface has to be constructed.

The fundamental equations of the dynamic theory for the two beam case are:

\begin{eqnarray}
-\gamma^{(j)} C_0^{j} +& \frac{U_g}{2K} C_g^{(j)} &= 0\\
\frac{U_g}{2K} C_0^{j} +& (-\gamma^{(j)}+s) C_g^{(j)} &= 0
\end{eqnarray}

Such a homogeneous linear equation system for the $C_g^{(j)}$ has a non-zero solution of and only if the determinant of the coefficients is zero:

\begin{equation}
\left|
\begin{array}{*{2}{c}}
-\gamma^{(j)} & \frac{U_g}{2K}\\
\frac{U_g}{2K}& (-\gamma^{(j)}+s)\\
\end{array}
\right| = {\gamma^{(j)}}^2 -s\gamma^{(j)}-\frac{U_g^2}{4K^2} =0
\end{equation}

Which is the same as the Howie-Whelan equation \ref{HaW3} with $\xi_g = K^2/U_g$,
but now we know that the $\gamma^{(j)}$ are the Eigenvalues of a matrix problem.

Solution:
\begin{eqnarray}
\gamma^{(j)}  &=& \frac{1}{2} \left[ s-(-1)^j \sqrt{(U_g/K)^2 +s^2}   \right]\\
&=& \frac{1}{2} \left[ s-(-1)^j \sqrt{(1/\xi^2 +s^2}   \right]\\
&=& \frac{1}{2\xi_g} \left[ w-(-1)^j \sqrt{(1+w^2}   \right]\\
\end{eqnarray}

We made the substitution $w =s\xi_g$, in which the parameter $w$ characterizes the tilt out of the exact Bragg condition ($w=0$). The excitation error $s$ is zero in the exact Bragg condition, isn't it.

 The separation is $\Delta k_{z,min} = \gamma^{(1)}-\gamma^{(2)}= \frac{U_g}{2K}=\frac{1}{\xi_g}$.

By use of the eigenvalues $\gamma^{(j)}$, the linear systems of equations can be solved for the $C_g^{(j)}$. For the amplitude $\epsilon^{(j)} C_g^{(j)} = C_0^{j} C_g^{(j)}$ of the four Bloch waves with the vector $\vec{k}_0^{(j)}+\vec{g}$ we obtain:
\begin{eqnarray}
C_0^{j} C_0^{(j)} &=& \frac{1}{2}\left[ 1 + (-1)^j \frac{w}{\sqrt{1+w^2}} \right]\\
C_0^{j} C_g^{(j)} &=& -\frac{1}{2}\left[  \frac{(-1)^j}{\sqrt{1+w^2}} \right]
\end{eqnarray}

We put this into the equation for the scattered wave and substitute the thickness $t$ for the $z$ component of the vector $\vec{r}$:
\begin{eqnarray}
\psi_0(t)&=&\sum_{j=1}^2 C_0^{j} C_0^{(j)}e^{2\pi i k_z^{(j)}t}\\
\psi_g(t)&=&\sum_{j=1}^2 C_0^{j} C_g^{(j)}e^{2\pi i k_z^{(j)}t}\, e^{2\pi i k_z gx}
\end{eqnarray}
and we find (omitting common phase factors):
\begin{eqnarray}
\psi_0(t)= \cos (\pi \sqrt{1+w^2}\frac{t}{\xi_g}) & - \frac{iw}{\sqrt{1+w^2}} \sin (\pi \sqrt{1+w^2}\frac{t}{\xi_g})\\
\psi_g(t)=  &  \frac{i}{\sqrt{1+w^2}} \sin (\pi \sqrt{1+w^2}\frac{t}{\xi_g}) \, e^{2\pi i k_z gx}
\end{eqnarray}

The intensities of the transmission $T$ and reflection $R$ become:
\begin{equation}
\underbrace{\psi_g \psi_g^*}_R = \underbrace{1-\psi_0 \psi_0^*}_{1-T} =  \frac{1}{1+w^2} \sin^2 (\pi \sqrt{\sqrt{1+w^2} ( \frac{1}{\xi_g}})_{s_{eff}} \cdot t)
\end{equation}

The solution is the Pendellösung of two coupled oscillators (in mechanics: two pendulums connected with a spring).

Even in exact Bragg condition, the intensity oscillates between primary beam and Bragg reflected beam with increasing film thickness. Look at the plot below, the Pendellösung is shown without absorption.

Normally one would want to add an absorption term to reduce the intensity with thickness. This absorption term is better named a damping term and stems form the inelatic scattering to random angles instead of the  considered (here two)  Bragg angles

These oscillations of the intensities are commonly called ``rocking curve``.

```python
# ------ Input ------
xi_g = 4  # extiction distance (in terms of relative thickness )
omega = 0 # tilt from Bragg condition
damping = 0.
# --------------------

t = np.linspace(0,8,401)

plt.figure()
plt.plot(t, (1-np.sin(np.pi * np.sqrt(np.sqrt(1.+ omega**2)*1/xi_g)*t)**2) * 1/np.exp(t*damping), label='incident beam')
plt.plot(t, np.sin(np.pi * np.sqrt(np.sqrt(1.+ omega**2)*1/xi_g)*t)**2 * 1/np.exp(t*damping), label='reflected beam')

plt.legend();
```

    <IPython.core.display.Javascript object>

## Summary of Bethe Theory
The solution is the Pendellösung of two coupled oscillators.

The periodicity is the extinction length $\xi$, which tells at which thickness a beam is completely vanished.

Considering some absorption (well it's not a real absorption, but inelastic scattering) then we see that the amplitudes decrease slowly.

## Using Bethe Theory for Thickness Determination

We will do this in a lab and it will be your homework.

- The accurate thickness of the sample is an important but hard to obtain parameter, but it influences the contrast in all imaging modes.
- Be aware that with different techniques you perform different thickness measurement. In any high resolution image and diffraction experiment, you always look at the thickness of the crystalline part of the sample, omitting the contribution of contamination and amorphous surface layer (from sample preparation).
- In the Analytic Section of this class we learn how to the thickness from the  whole sample.

We can observe the above rocking curve in convergent beam electron diffraction patterns (CBED).

But we have to ensure that:
- Excitation error is as small as possible
- We are in two beam condition

### Experimental Considerations
- Choose a convergence angle $\alpha$ so that $\alpha < \theta_B$, to avoid overlapping of disks in the ZOLZ.
- The 000 disk usually contains concentric diffuse fringes, the Kossel-Möllenstedt fringes

- If you move the specimen, then you will see that the number of this fringes changes. In fact the number of each fringes increases by one every time the thickness increases by one extinction length.
- The foil thickness can be measured precisely at the point where you do your other analysis.

Please be aware that dynamic effects also occur for the HOLZ lines in a CBED pattern.

In practice to simplify the interpretation, we don't use zone axis conditions, but tilt to two--beam conditions with only one strongly excited Bragg beam.
- The CBED disks contain then parallel rather than concentric intensity oscillations as shown in the earlier figure.
- In fact, this intensity oscillations are equivalent to the rocking curve intensity oscillations discussed earlier.
- It helps, of you use an energy filter for this method.

### Thickness Determination
Because the oscillations are symmetric in the hkl disk we concentrate the analysis on this disk.
- The middle of the hkl disk is bright and originates from the exact Bragg condition ($\vec{s}=0$).
- We measure the distance between the midle (bright fringe) of the $hkl$ disk and the dark lines.

You obtain a deviation $s_i$ for each fringe from the equation:
\begin{equation}
s_i=\lambda\frac{\Delta\theta_i}{e\theta_B d^2}
\end{equation}
The Bragg angle $\theta_B$ is known from the separation of two disks and the lattice spacing $d$ is known from the sample or can be calculated through the camera length.

If the extinction distance $\xi_g$ is known you can calculate the foil thickness $t$ with:
\begin{equation}
\frac{1}{t^2} = \frac{s_1^2}{n_k^2}+\frac{1}{\xi_g^2n_k^2}
\end{equation}
where $n_k$ is an integer.

### Data Analysis
- assign $n=1$ to the first fringe $s_1$
- assign $n=2$ to the second fringe $s_2$
  and so on for all other fringes
- plot $(s_1/n_k)^2$ versus $(1/n_k)^2$.
- if you get a straight line, then you are finished and you have $k=i+j$, where $j$ is the largest integer $<(t/\xi_g)$.
- if not repeat the same thing with $n=2$ for $s_1$, $n=3$ for $s_2$, etc.
- repeat this increase by one till you get a straight line
- the slope of the line is $1/\xi_g^2$
- the extrapolated value for $1/n_k^2$ is $1/t^2$.

The whole procedure is summarized in the figure below.

## More about Bloch Waves
We can replace the  exponential functions by trigonometric functions and get:
\begin{eqnarray}
A^{(1)} &=& \cos \frac{\beta}{2}\\
A^{(2)} &=& \sin \frac{\beta}{2}
\end{eqnarray}

Some of the Bloch waves are located (have their maxima) between the atoms. These Bloch waves channel and are more or less undisturbed.

Another set is located on the atomic rows and will cause much more inelastic scattering than the other, also they will travel much faster.

The second set is especially important for Z-contrast image, where a small convergent beam is located at the atomic rows. You might consider the atoms like little lenses which keep the beam focused on the column.

## Navigation

-   **Next: [Mutli-Slice Theory](CH2_D02-Multi_Slice.ipynb)**
-   **Chapter 2: [Diffraction](CH2_00-Diffraction.ipynb)**
-   **List of Content: [Front](../_MSE672_Intro_TEM.ipynb)**

```python

```


text:
```
import numpy.matlib
import numpy as np

```

```
import random
A=-1
B= 1
C= 17
i= 0
n= 0
for i in range(10000):
  a= random.uniform(A,B)
  b= random.uniform(A,B)
  if a**2+b**2<1:
    n=n+1
print(str(4*n/10000))
```

    3.14

# Example 3
Given i.i.d $\{\alpha_i : i \in 1, 2, ..., N\}$, we use
$$\bar \alpha_N = \frac{1}{N} \sum_{i=1}^N \alpha_i$$
as its estimator of the mean $\mathbb E[\alpha_1]$ and use
$$\beta_N = \frac{1}{N} \sum_{i=1}^N (\alpha_i - \bar \alpha_N)^2$$
as the estimator of $Var(\alpha_1)$. Suppose $\alpha_1 \in L^4$, then
- Prove $\beta_N$ is biased.
- Prove that $\beta_N$ is consistent in $L^2$.
- Can you propose an unbiased estimator?

__Proof:__
because $\{\alpha_i\}$ is  iid
\begin{equation}
\begin{aligned}
Bias（\beta_N）= &= \mathbb E[\beta_N] - Var(\alpha_1) \\
\end{aligned}
\end{equation}
so
\begin{equation}
\begin{aligned}
Bias(\beta_N) &= \mathbb E[\beta_N] - Var(\alpha_1) \\
&= \mathbb E[\frac{1}{N} \sum_{i=1}^N (\alpha_i - \bar \alpha_N)^2] - Var(\alpha_1) \\
&= \frac{1}{N} \mathbb E[\sum_{i=1}^N (\alpha_i^2- 2 \alpha_i \bar \alpha_N + \bar \alpha_N^2)] - Var(\alpha_1) \\
&= \frac{1}{N} \mathbb E[\sum_{i=1}^N \alpha_i^2 - 2N \cdot \bar \alpha_N^2 + N \cdot \bar \alpha_N^2] - Var(\alpha_1) \\
&= \frac{1}{N} \sum_{i=1}^N \mathbb E[\alpha_i^2] - \frac{1}{N}\mathbb E[N \cdot \bar \alpha_N^2] - Var(\alpha_1) \\
&= \mathbb E[\alpha_1^2] - \mathbb E[\bar \alpha_N^2]- Var(\alpha_1) \\
&= \mathbb E[\alpha_1]^2 - \mathbb E[\bar \alpha_N^2] .......................(1)\\
\end{aligned}
\end{equation}
because
\begin{equation}
\begin{aligned}
\mathbb E[\bar \alpha_N] &= \mathbb E[\frac{1}{N} \sum_{i=1}^N \alpha_i] \\
&= \frac{1}{N} \sum_{i=1}^N \mathbb E[\alpha_i] \\
&= \frac{1}{N} \cdot N \cdot \mathbb E[\alpha_1] \\
&= \mathbb E[\alpha_1], \\
\end{aligned}
\end{equation}
so
\begin{equation}
\begin{aligned}
(1)=E[\bar \alpha_N]^2-E[\bar \alpha_N^2]=-Var(\bar \alpha_N)=-Var(\frac{1}{N} \sum_{i=1}^N \alpha_i)= -\frac{1}{N^2} \sum_{i=1}^N Var(\alpha_i)=- \frac{1}{N^2} \cdot N \cdot Var(\alpha_1) = -\frac{1}{N} Var(\alpha_1)\\
\end{aligned}
\end{equation}
βN  is biased

 - To show βN is consistent in L^2, we need to show that
$$\mathbb E[(\beta_N - Var(\alpha_1))^2] \to 0.$$

- find the unbiased estimator:
we already know if $$\beta_N = \frac{1}{N} \sum_{i=1}^N (\alpha_i - \bar \alpha_N)^2$$
\begin{equation}
\begin{aligned}
Bias（\beta_N）= \mathbb E[\beta_N] - Var(\alpha_1) =-\frac{1}{N} Var(\alpha_1)\\
\end{aligned}
\end{equation}
so
\begin{equation}
\begin{aligned}
\mathbb E[\beta_N] = \frac{N-1}{N} Var(\alpha_1)\\
\mathbb E[\beta_2]= \frac{N}{N-1}E[\beta_N] \\
\beta_2= \frac{1}{N-1} \sum_{i=1}^N (\alpha_i - \bar \alpha_N)^2
\end{aligned}
\end{equation}

```
D=0
m=0
pi_list = []
for i in range(100):
  A=-1
  B= 1
  C= 17
  i= 0
  n= 0
  for z in range(10000):
      a= random.uniform(A,B)
      b= random.uniform(A,B)
      if a**2+b**2<1:
       n=n+1
  pi_list.append(4*n/10000)
  D=4*n/10000+D
d= D/100
for i in range(100):
  m = (pi_list[i] - d)**2+m
print(str(m/100))
```

    0.00029377369600000047

```
import matplotlib.pyplot as plt
y = []
for t in range (6):
  D=0
  m=0
  pi_list = []
  for i in range(100):
    A=-1
    B= 1
    C= 17
    i= 0
    n= 0
    for z in range((2**(5+t))):
        a= random.uniform(A,B)
        b= random.uniform(A,B)
        if a**2+b**2<1:
         n=n+1
    pi_list.append(4*n/(2**(5+t)))
    D=4*n/(2**(5+t))+D
  d= D/100
  for i in range(100):
    m = (pi_list[i] - d)**2+m
  y.append(m/100)
x = [n for n in range(5,11)]
print(y)
plt.loglog(x,y)
```

```
import random
A= 0
B= 1
C= 17
i= 0
n= 0
for i in range(100000000):
  a= random.uniform(A,B)
  if a<0.01:
    n=n+100
  if a>0.01:
    n=n+1
print(str(n/100000000))
```

    1.98990892

```
import matplotlib.pyplot as plt
y = []
for t in range (6):
  D=0
  m=0
  i_list = []
  for i in range(100):
    A= 0
    B= 1
    C= 17
    i= 0
    n= 0
    for z in range((2**(5+t))):
        a= random.uniform(A,B)
        if a<0.01:
          n=n+100
        if a>0.01:
          n=n+1
    i_list.append(n/(2**(5+t)))
    D=n/(2**(5+t))+D
  d= D/100
  for i in range(100):
    m = (i_list[i] - d)**2+m
  y.append(m/100)
x = [n for n in range(5,11)]
print(y)
plt.loglog(x,y)
```

ext: ipynb

max_stars_repo_path: hw5/omc_integral and mcpi.ipynb

max_issues_repo_path: hw5/omc_integral and mcpi.ipynb

max_forks_repo_path: hw5/omc_integral and mcpi.ipynb


text:
# Gaussian Process (GP) smoothing

This example deals with the case when we want to **smooth** the observed data points $(x_i, y_i)$ of some 1-dimensional function $y=f(x)$, by finding the new values $(x_i, y'_i)$ such that the new data is more "smooth" (see more on the definition of smoothness through allocation of variance in the model description below) when moving along the $x$ axis.

It is important to note that we are **not** dealing with the problem of interpolating the function $y=f(x)$ at the unknown values of $x$. Such problem would be called "regression" not "smoothing", and will be considered in other examples.

If we assume the functional dependency between $x$ and $y$ is **linear** then, by making the independence and normality assumptions about the noise, we can infer a straight line that approximates the dependency between the variables, i.e. perform a linear regression. We can also fit more complex functional dependencies (like quadratic, cubic, etc), if we know the functional form of the dependency in advance.

However, the **functional form** of $y=f(x)$ is **not always known in advance**, and it might be hard to choose which one to fit, given the data. For example, you wouldn't necessarily know which function to use, given the following observed data. Assume you haven't seen the formula that generated it:

```python
%pylab inline
figsize(12, 6);
```

    Populating the interactive namespace from numpy and matplotlib

```python
import numpy as np
import scipy.stats as stats

x = np.linspace(0, 50, 100)
y = (np.exp(1.0 + np.power(x, 0.5) - np.exp(x/15.0)) +
     np.random.normal(scale=1.0, size=x.shape))

plot(x, y);
xlabel("x");
ylabel("y");
title("Observed Data");
```

### Let's try a linear regression first

As humans, we see that there is a non-linear dependency with some noise, and we would like to capture that dependency. If we perform a linear regression, we see that the "smoothed" data is less than satisfactory:

```python
plot(x, y);
xlabel("x");
ylabel("y");

lin = stats.linregress(x, y)
plot(x, lin.intercept + lin.slope * x);
title("Linear Smoothing");
```

### Linear regression model recap

The linear regression assumes there is a linear dependency between the input $x$ and output $y$, sprinkled with some noise around it so that for each observed data point we have:

$$ y_i = a + b\, x_i + \epsilon_i $$

where the observation errors at each data point satisfy:

$$ \epsilon_i \sim N(0, \sigma^2) $$

with the same $\sigma$, and the errors are independent:

$$ cov(\epsilon_i, \epsilon_j) = 0 \: \text{ for } i \neq j $$

The parameters of this model are $a$, $b$, and $\sigma$. It turns out that, under these assumptions, the maximum likelihood estimates of $a$ and $b$ don't depend on $\sigma$. Then $\sigma$ can be estimated separately, after finding the most likely values for $a$ and $b$.

### Gaussian Process smoothing model

This model allows departure from the linear dependency by assuming that the dependency between $x$ and $y$ is a Brownian motion over the domain of $x$. This doesn't go as far as assuming a particular functional dependency between the variables. Instead, by **controlling the standard deviation of the unobserved Brownian motion** we can achieve different levels of smoothness of the recovered functional dependency at the original data points.

The particular model we are going to discuss assumes that the observed data points are **evenly spaced** across the domain of $x$, and therefore can be indexed by $i=1,\dots,N$ without the loss of generality. The model is described as follows:

\begin{equation}
\begin{aligned}
z_i & \sim \mathcal{N}(z_{i-1} + \mu, (1 - \alpha)\cdot\sigma^2) \: \text{ for } i=2,\dots,N \\
z_1 & \sim ImproperFlat(-\infty,\infty) \\
y_i & \sim \mathcal{N}(z_i, \alpha\cdot\sigma^2)
\end{aligned}
\end{equation}

where $z$ is the hidden Brownian motion, $y$ is the observed data, and the total variance $\sigma^2$ of each ovservation is split between the hidden Brownian motion and the noise in proportions of $1 - \alpha$ and $\alpha$ respectively, with parameter $0 < \alpha < 1$ specifying the degree of smoothing.

When we estimate the maximum likelihood values of the hidden process $z_i$ at each of the data points, $i=1,\dots,N$, these values provide an approximation of the functional dependency $y=f(x)$ as $\mathrm{E}\,[f(x_i)] = z_i$ at the original data points $x_i$ only. Therefore, again, the method is called smoothing and not regression.

### Let's describe the above GP-smoothing model in PyMC3

```python
import pymc3 as pm
from theano import shared
from pymc3.distributions.timeseries import GaussianRandomWalk
from scipy import optimize
```

    /anaconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
      from ._conv import register_converters as _register_converters

Let's create a model with a shared parameter for specifying different levels of smoothing. We use very wide priors for the "mu" and "tau" parameters of the hidden Brownian motion, which you can adjust according to your application.

```python
LARGE_NUMBER = 1e5

model = pm.Model()
with model:
    smoothing_param = shared(0.9)
    mu = pm.Normal("mu", sigma=LARGE_NUMBER)
    tau = pm.Exponential("tau", 1.0/LARGE_NUMBER)
    z = GaussianRandomWalk("z",
                           mu=mu,
                           tau=tau / (1.0 - smoothing_param),
                           shape=y.shape)
    obs = pm.Normal("obs",
                    mu=z,
                    tau=tau / smoothing_param,
                    observed=y)
```

    INFO (theano.gof.compilelock): Waiting for existing lock by process '16678' (I am process '16664')
    INFO (theano.gof.compilelock): To manually release the lock, delete /.theano/compiledir_Linux-4.15--generic-x86_64-with-debian-buster-sid-x86_64-3.6.8-64/lock_dir
    INFO (theano.gof.compilelock): Waiting for existing lock by process '16678' (I am process '16664')
    INFO (theano.gof.compilelock): To manually release the lock, delete /.theano/compiledir_Linux-4.15--generic-x86_64-with-debian-buster-sid-x86_64-3.6.8-64/lock_dir
    INFO (theano.gof.compilelock): Waiting for existing lock by process '16678' (I am process '16664')
    INFO (theano.gof.compilelock): To manually release the lock, delete /.theano/compiledir_Linux-4.15--generic-x86_64-with-debian-buster-sid-x86_64-3.6.8-64/lock_dir
    INFO (theano.gof.compilelock): Waiting for existing lock by process '16678' (I am process '16664')
    INFO (theano.gof.compilelock): To manually release the lock, delete /.theano/compiledir_Linux-4.15--generic-x86_64-with-debian-buster-sid-x86_64-3.6.8-64/lock_dir
    INFO (theano.gof.compilelock): Waiting for existing lock by process '16678' (I am process '16664')
    INFO (theano.gof.compilelock): To manually release the lock, delete /.theano/compiledir_Linux-4.15--generic-x86_64-with-debian-buster-sid-x86_64-3.6.8-64/lock_dir

Let's also make a helper function for inferring the most likely values of $z$:

```python
def infer_z(smoothing):
    with model:
        smoothing_param.set_value(smoothing)
        res = pm.find_MAP(vars=[z], fmin=optimize.fmin_l_bfgs_b)
        return res['z']
```

Please note that in this example, we are only looking at the MAP estimate of the unobserved variables. We are not really interested in inferring the posterior distributions. Instead, we have a control parameter $\alpha$ which lets us allocate the variance between the hidden Brownian motion and the noise. Other goals and/or different models may require sampling to obtain the posterior distributions, but for our goal a MAP estimate will suffice.

### Exploring different levels of smoothing

Let's try to allocate 50% variance to the noise, and see if the result matches our expectations.

```python
smoothing = 0.5
z_val = infer_z(smoothing)

plot(x, y);
plot(x, z_val);
title("Smoothing={}".format(smoothing));
```

It appears that the variance is split evenly between the noise and the hidden process, as expected.

Let's try gradually increasing the smoothness parameter to see if we can obtain smoother data:

```python
smoothing = 0.9
z_val = infer_z(smoothing)

plot(x, y);
plot(x, z_val);
title("Smoothing={}".format(smoothing));
```

### Smoothing "to the limits"

By increading the smoothing parameter, we can gradually make the inferred values of the hidden Brownian motion approach the average value of the data. This is because as we increase the smoothing parameter, we allow less and less of the variance to be allocated to the Brownian motion, so eventually it aproaches the process which almost doesn't change over the domain of $x$:

```python
fig, axes = subplots(2, 2)

for ax, smoothing in zip(axes.ravel(), [0.95, 0.99, 0.999, 0.9999]):

    z_val = infer_z(smoothing)

    ax.plot(x, y)
    ax.plot(x, z_val)
    ax.set_title('Smoothing={:05.4f}'.format(smoothing))
```

This example originally contributed by: Andrey Kuzmenko, http://github.com/akuz

```python

```


text:
```python
from IPython.core.display import display_html
from urllib.request import urlopen

link = 'https://gist.github.com/robblack007/eca03fa9f7586860235d/raw/ef05a2f29febc94a9c9f99ca20fd0b65e74ed451/custom.css'
display_html(urlopen(link).read(), raw=True)
```

<link href='http://fonts.googleapis.com/css?family=EB+Garamond' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Lato:100,100italic' rel='stylesheet' type='text/css'>
<link rel="stylesheet" href="//cdn.jsdelivr.net/font-hack/2.013/css/hack-extended.min.css">

<style>
	@font-face{
   		font-family: "Computer Modern";
        	src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');}

        div.cell{
	        max-width:    800px;
	        margin-left:  auto;
	        margin-right: auto;}

	h1 {
        	font-family: 'Lato', sans-serif;}

	h2 {
        	font-family: 'Lato', serif;}

	h3{
		font-family:   'Lato', serif;
        	margin-top:    12px;
		margin-bottom: 3px;}

        h4{
		font-family: 'Lato', serif;}

	h5 {
        	font-family: 'Lato', serif;}

	div.text_cell_render{
        	font-family:  'EB Garamond', 'Computer Modern', 'Helvetica Neue', Arial, Helvetica, Geneva, sans-serif;
	        line-height:  140%;
	        font-size:    14pt;
		max-width:    650px;
		margin-left:  auto;
        	margin-right: auto;}

	.CodeMirror{
        	font-family: 'Inconsolata', monospace;}

	.text_cell_render h1 {
	        font-weight:   100;
	        font-size:     48pt;
		line-height:   100%;
	        color:         #CD2305;
	        margin-bottom: 0.5em;
	        margin-top:    0.5em;
	        display:       block;}

	.text_cell_render h2 {
        	font-weight:   100;
        	font-size:     36pt;}

	.text_cell_render h3 {
        	font-weight:   100;
        	font-size:     21pt;}

	.text_cell_render h4 {
        	font-weight:   100;
        	font-size:     18pt;}

	.text_cell_render h5 {
	        font-weight:   100;
	        font-size:     16pt;
	        color:         #CD2305;
	        font-style:    italic;
	        margin-bottom: 0.5em;
	        margin-top:    0.5em;
	        display:       block;}

	.text_cell_render h6 {
	        font-weight:   100;
	        font-size:     14pt;
	        color:         #CD2305;
	        font-style:    italic;
	        margin-bottom: 0.5em;
	        margin-top:    0.5em;
	        display:       block;}

	.warning{
	    	background-color: #FFFFE0;
	    	border-color: #FFCC99;
	    	border-left: 5px solid #FFCC99;
	    	padding: 0.5em;}

	.error{
	    	background-color: #fcf2f2;
	    	border-color: #dFb5b4;
	    	border-left: 5px solid #dfb5b4;
	    	padding: 0.5em;}

</style>

# Dinámica de un robot manipulador planar

```python
from sympy import var, sin, cos, Matrix, Integer, eye, Function, Rational, exp, Symbol, I, solve, pi, trigsimp, dsolve, sinh, cosh, simplify
from sympy.physics.mechanics import mechanics_printing
mechanics_printing()
```

```python
var("m1 m2 J1 J2 L1 L2 l1 l2 t g")
```

```python
q1 = Function("q1")(t)
q2 = Function("q2")(t)
```

```python
x1 = l1*cos(q1)
y1 = l1*sin(q1)
v1 = x1.diff("t")**2 + y1.diff("t")**2
v1.trigsimp()
```

```python
x2 = L1*cos(q1) + l2*cos(q1 + q2)
y2 = L1*sin(q1) + l2*sin(q1 + q2)
v2 = x2.diff("t")**2 + y2.diff("t")**2
v2.trigsimp()
```

```python
ω1 = q1.diff("t")
ω2 = q2.diff("t")
```

```python
K1 = Rational(1, 2)*m1*v1 + Rational(1, 2)*J1*ω1**2
K1
```

```python
K2 = Rational(1, 2)*m1*v2 + Rational(1, 2)*J2*ω2**2
K2
```

```python
U1 = m1*g*y1
U1
```

```python
U2 = m2*g*y2
U2
```

```python
K = K1 + K2
K
```

```python
U = U1 + U2
U
```

```python
L = (K - U).expand().simplify()
L
```

```python
τ1 = (L.diff(q1.diff(t)).diff(t) - L.diff(q1)).simplify().expand().collect(q1.diff(t).diff(t)).collect(q2.diff(t).diff(t))
```

```python
τ2 = (L.diff(q2.diff(t)).diff(t) - L.diff(q2)).simplify().expand().collect(q1.diff(t).diff(t)).collect(q2.diff(t).diff(t))
```

```python
τ1
```

```python
τ2
```

```python
from scipy.integrate import odeint
from numpy import linspace
```

```python
def pendulo_doble(estado, tiempo):
    # Se importan funciones necesarias
    from numpy import sin, cos, matrix
    # Se desenvuelven variables del estado y tiempo
    q1, q2, q̇1, q̇2 = estado
    t = tiempo

    # Se declaran constantes del sistema
    m1, m2 = 1, 1
    l1, l2 = 1, 1
    g = 9.81

    # Se declaran constantes del control
    kp1, kp2 = -30, -60
    kv1, kv2 = -20, -20

    # Señales de control nulas
    tau1, tau2 = 0, 0

    # Posiciones a alcanzar
    qd1, qd2 = 1, 1

    # Se declaran señales de control del sistema
    #tau1 = kp1*(q1 - qd1) + kv1*q̇1
    #tau2 = kp2*(q2 - qd2) + kv2*q̇2

    # Se calculan algunos terminos comunes
    ϕ1 = m1*l1**2
    ϕ2 = m1*l1*l2
    ϕ3 = m1*l2**2

    # Se calculan las matrices de masas, Coriolis,
    # y vectores de gravedad, control, posicion y velocidad
    M = matrix([[2*ϕ1 + 2*ϕ2*cos(q2) + ϕ3, ϕ2*cos(q2) + ϕ3],
                [ϕ2*cos(q2) + ϕ3, ϕ3]])
    C = matrix([[-2*ϕ2*sin(q2)*q̇2, -ϕ2*sin(q2)*q̇2], [ϕ2*sin(q2)*q̇1, 0]])
    G = matrix([[m1*l1*cos(q1) + m2*l1*cos(q1) + m2*l2*cos(q1 + q2)], [m2*l2*cos(q1 + q2)]])
    Tau = matrix([[tau1], [tau2]])
    q = matrix([[q1], [q2]])
    q̇ = matrix([[q̇1], [q̇2]])

    # Se calcula la derivada del estado del sistema
    qp1 = q̇1
    qp2 = q̇2

    qpp = M.I*(Tau - C*q̇ - G)
    qpp1, qpp2 = qpp.tolist()

    return [qp1, qp2, qpp1[0], qpp2[0]]
```

```python
t = linspace(0, 10, 1000)
estados_simulados = odeint(func = pendulo_doble, y0  = [0, 0, 0, 0], t = t)
```

```python
q1, q2, q̇1, q̇2 = list(zip(*estados_simulados.tolist()))
```

```python
%matplotlib notebook
from matplotlib.pyplot import plot, style, figure
from mpl_toolkits.mplot3d import Axes3D
style.use("ggplot")
```

```python
fig1 = figure(figsize=(8, 8))

ax1 = fig1.gca()

p1, = ax1.plot(t, q1)
p2, = ax1.plot(t, q2)
ax1.legend([p1, p2],[r"$q_1$", r"$q_2$"])
ax1.set_ylim(-4, 4)
ax1.set_xlim(-0.1, 10);
```

    <IPython.core.display.Javascript object>

```python

```

max_stars_repo_path: Clases/.ipynb_checkpoints/Dinamica pendulo doble completa-checkpoint.ipynb


text:
# Bayesian Hyperparameter Optimization using Gaussian Processes

TODO: animation of GP optimization

TODO: intro

TODO: outline

First, let's import the packages we'll use, and set some settings.

```
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import RobustScaler

# Plot settings
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
mpl.rcParams['image.cmap'] = 'viridis'
```

We'll also be using [CatBoost](https://catboost.ai/) for prediction, so let's install and import that as well.

```
!pip install -q catboost
```

        100% |████████████████████████████████| 60.1MB 545kB/s

```
from catboost import CatBoostRegressor
```

A more complete version of the code we'll be working with here is available [on GitHub](http://github.com/brendanhasz/dsutils).  We'll also import the two relevant components from that repository which this post shows how to build:

```
!pip install git+http://github.com/brendanhasz/dsutils.git
from dsutils.optimization import GaussianProcessOptimizer
from dsutils.optimization import optimize_params
```

    Collecting git+http://github.com/brendanhasz/dsutils.git
      Cloning http://github.com/brendanhasz/dsutils.git to /tmp/pip-req-build-4ktk2eyd
    Building wheels for collected packages: dsutils
      Building wheel for dsutils (setup.py) ... done
      Stored in directory: /tmp/pip-ephem-wheel-cache-y_8w1syc/wheels/a1/ff/2a/75bdc08e9c96d4917294db5e6faf99ef3de673f37992c52278
    Successfully built dsutils
    Installing collected packages: dsutils
    Successfully installed dsutils-0.1

## The Math

TODO: gaussian process, expensive function, etc

To determine what point to sample next, we can find the hyperparameters which maximize an "acquisition function".  The acquisition function is a function that when evaluated at $\mathbf{x}$ (some combination of hyperparameters), tells us how advantagous it would be to evaluate our expensive function $f$ at that point.  That way, we can use the acquisition function to find the next combination of hyperparameters we should try!

There are a few different functions which we could use for the acquisition function, including the probability of improvement, the expected improvement, and the upper confidence bound. Here we'll use the expected improvement as our acquisition function.

The expected improvement for some combination of hyperparameters $\mathbf{x}$ is, like the name suggests, the average improvement that we expect from that combination of hyperparameters over the best set we've found so far ($\hat{\mathbf{x}}$):

$$
EI(\mathbf{x}) = \mathbb{E}[\max (0, ~ f(\mathbf{x})-f(\hat{\mathbf{x}})) ]
$$

where $f(\mathbf{x})$ is our Gaussian process' estimation as to our model's performance with hyperparameter combination $\mathbf{x}$, and $f(\hat{\mathbf{x}})$ is the best performance that we've actually acheived so far.

How do we actually compute that?  The Gaussian process is modeling the probability of our model's performance ($f$) as a function of the hyperparameters ($\mathbf{x}$), and so at any value of $\mathbf{x}$, it gives us a normal distribution for its prediction of the value of $f$:

$$
p(f(\mathbf{x})|\mathbf{x}) ~ \sim ~ \mathcal{N}(\mu(\mathbf{x}), \sigma(\mathbf{x}))
$$

So, we can integrate that probability distribution (times the improvement magnitude) above the value of our current best performance value.

TODO: diagram

This gives us how much of an improvement we can expect with hyperparameters $\mathbf{x}$ over our current best, $\hat{\mathbf{x}}$.

$$
\begin{align}
EI(\mathbf{x}) &= \int_{f(\hat{\mathbf{x}})}^\infty p(f(\mathbf{x})|\mathbf{x}) (f(\mathbf{x})-f(\hat{\mathbf{x}}))  df(\mathbf{x}) \\
&= \int_{f(\hat{\mathbf{x}})}^\infty\mathcal{N}(f(\mathbf{x})|\mu(\mathbf{x}), \sigma(\mathbf{x})) ~ (f(\mathbf{x})-f(\hat{\mathbf{x}}))  ~ df(\mathbf{x}) \\
\end{align}
$$

After "applying some tedious integration by parts" ([Jones et al., 1998](https://link.springer.com/article/10.1023/A:1008306431147)),  we have an analytical solution:

$$
EI(\mathbf{x}) = (\mu(\mathbf{x})-f(\hat{\mathbf{x}}))\Phi(z) + \sigma(\mathbf{x}) \phi(z)
$$

where $\Phi$ is the cumulative distribution function of the standard normal distribution, $\phi$ is the probability density function of the standard normal, and $z$ is the scaled difference between the Gaussian process mean and our current best performance:

$$
z = \frac{\mu(\mathbf{x})-f(\hat{\mathbf{x}})}{\sigma(\mathbf{x})}
$$

Using our Gaussian process to estimate $\mu(\mathbf{x})$ and $\sigma(\mathbf{x})$, we can then search for the $\mathbf{x}$ which we expect to give us the best expected improvement if we take our next sample there:

$$
\mathbf{x}_{t+1} = \arg \max_{\mathbf{x}_{t+1}} EI(\mathbf{x}_{t+1})
$$

For more information,
[Brochu et al., 2010](https://arxiv.org/abs/1012.2599) is a great tutorial on Bayesian optimization, which includes an intro to Gaussian processes and info about several different types of acquisition functions.  But enough math - on to the code!

## Bayesian Optimization

Bayesian optimization isn't specific to finding hyperparameters - it lets you optimize *any* expensive function.  That includes, say, the parameters of a simulation which takes a long time, or the configuration of a scientific research study, or the appearance of a website during an A/B test.  Before we write code to optimize our hyperparameters, we'll build a system which can optimize any expensive function using Bayesian optimization.

### The code

To optimize our expensive function, let's build a python class which we can use to keep track of the hyperparameter combinations we've sampled, the resulting error, and then suggest the next combination of hyperparameters that we should try based on that information.

When initializing the optimizer, the user should define the lower and upper bounds of the parameters to optimize (so our optimizer can search some confined space!).  Also, the user should be able to specify whether they want to maximize or minimize the function (we'll do that via a `minimize` arg).

For example, to create an optimizer which finds the maximum of a function with 2 parameters, the first a float between 0 and 1, and the second an integer between 10 and 100, the user should be able to do:

```
gpo = GaussianProcessOptimizer([0, 10], [1, 100],
                               minimize=False)
```

So to start, let's build the class definition and the constructor (which defines how one initializes the object, as above):

```
class GaussianProcessOptimizer():
    """Bayesian function optimizer which uses a Gaussian process to model the
    expensive function, and expected improvement as the acquisition function.
    """

    def __init__(self, lb, ub, minimize=True):
        """Gaussian process-based optimizer

        Parameters
        ----------
        lb : list
            Lower bound for each parameter.
        ub : list
            Upper bound for each parameter.
        minimize : bool
            Whether to minimize (True) or maximize (False).
        """

        # Store parameters
        self.lb = lb
        self.ub = ub
        self.minimize = minimize
        self.db = [ub[i]-lb[i] for i in range(self.num_dims)]
        self.bounds = [(lb[i], ub[i]) for i in range(self.num_dims)]

        # Gaussian process to use for estimating the function
        self.gp = GaussianProcessRegressor(
            kernel=RationalQuadratic()+WhiteKernel(noise_level=1e-3),
            alpha=0.0,
            n_restarts_optimizer=10,
        )

        # Keep track of points so far, and the best
        self.x = []
        self.y = []
        self.opt_x = None
        if minimize:
            self.opt_y = np.inf
        else:
            self.opt_y = -np.inf
```

The class will also need a function to add points to the history, and the resulting function value corresponding to that point.  This function should allow multiple y values per x value, because when we use it for estimating the cross-validated scores, we'll have a different score for each fold!  Also, the function should keep track of the point corresponding to the best y value which has been sampled so far.  We'll add this function to the `GaussianProcessOptimizer` class:

```
def add_point(self, x, y):
    """Add a point to the history of sampled points."""

    # Append to sample record
    if isinstance(y, list): #repeated x values
        for ty in y:
            self.x.append(x)
            self.y.append(ty)
        ty = np.array(y).mean()
    else:
        self.x.append(x)
        self.y.append(y)
        ty = y

    # Store best point so far
    if self.minimize:
        if ty < self.opt_y:
            self.opt_y = ty
            self.opt_x = x
    else:
        if ty > self.opt_y:
            self.opt_y = ty
            self.opt_x = x
```

The meat of the optimizer is its ability to fit a Gaussian process to the sampled points, so we'll add a private function to do that fitting.  We'll normalize the y-values, which will cause our Gaussian process to assume that, in areas of parameter space which we have not explored yet, the y-values are roughly equal to the mean of the y-values for the points which we have sampled.  The way in which we normalize y will need to be kept track of, so that we can *un*-normalize the predictions of the Gaussian process.

We'll also normalize the x-values so that we have to worry less about the exact value of the length-scale parameter of the Gaussian process (see my [previous post](https://brendanhasz.github.io/2018/10/10/hmm-vs-gp.html) on Gaussian processes to see why that can be a pain).  Let's add to the `GaussianProcessOptimizer` class a method to fit the Gaussian process:

```
def _fit_gp(self, step=None):
    """Fit the Gaussian process to data."""

    # Normalize y
    self._y_mean = y.mean()
    self._y_std = y.std()
    y = (y-self._y_mean)/self._y_std

    # Normalize x
    for iD in range(self.num_dims):
        x[:, iD] = (x[:, iD]-self.lb[iD]) / self.db[iD]

    # Fit the Gaussian process
    self.gp = self.gp.fit(x, y)
```

Our optimizer will also need to be able use the Gaussian process to predict the y-values (e.g. the cross-validated performance) for a given x-value (e.g. the hyperparameter values).  We need to normalize the new x values in the same way we did when fitting the Gaussian process (above), and un-normalize the predicted y-values as discussed above.  Let's another private method to the class which performs that prediction:

```
def _pred_gp(self, x, return_std=False):
    """Predict y with the Gaussian process."""

    # Normalize x
    for iD in range(self.num_dims):
        x[:, iD] = (x[:, iD] - self.lb[iD]) / self.db[iD]

    # Predict y
    y, y_std = self.gp.predict(x, return_std=True)

    # Convert y back to true scale
    y = y*self._y_std+self._y_mean
    y_std = y_std*self._y_std

    # Return std dev if requested
    if return_std:
        return y, y_std
    else:
        return y
```

The last key peice of the optimizer is the acquisition function.  Let's add one last private method to the `GaussianProcessOptimizer` class which uses the Gaussian process' estimate and uncertainty as to the y-values to compute the expected improvement:

```
def _expected_improvement(self, x):
    """Compute the expected improvement at x."""

    # Predict performance at x
    mu, sigma = self._pred_gp(x.reshape(-1, self.num_dims),
                              return_std=True)

    # Compute and return expected improvement
    flip = np.power(-1, self.minimize)
    z = flip*(mu-self.opt_y)/sigma
    return flip*(mu-self.opt_y)*norm.cdf(z) + sigma*norm.pdf(z)
```

Now that our optimizer is able to fit the Gaussian process to our data, we need it to be able to suggest points to sample.  The easiest method of suggesting points is just generating a random point that's within the parameter bounds:

```
def random_point(self, get_dict=False):
    """Get a random point within the bounds."""
    return np.random.uniform(self.lb, self.ub)
```

But what we really want is for the optimizer to suggest a point which is either likely to have a good y-value, or a point for which we're highly uncertain as to the y-values in that area of parameter space.  We'll add to the `GaussianProcessOptimizer` class a method which first fits a Gaussian process to the points sampled so far, uses that fit Gaussian process to compute the expected improvement, and then searches for the x-values which result in the best expected improvement:

```
def next_point(self):
    """Get the point with the highest expected improvement."""

    # Fit the Gaussian process to samples so far
    self._fit_gp()

    # Find x with greatest expected improvement
    x = self.random_point()
    best_score = np.inf
    n_restarts = 10
    for iR in range(n_restarts):

        # Maximize expected improvement
        res = minimize(lambda x: -self._expected_improvement(x),
                       self.random_point(),
                       method='L-BFGS-B',
                       bounds=self.bounds)

        # Keep x if it's the best so far
        if res.fun < best_score:
            best_score = res.fun
            x = res.x

    # Return x with highest expected improvement
    return x
```

After sampling as many points as we have time to sample, we'll want the optimizer to spit out the point which results in the best y-value.  There's two ways we could do this.  The simplest way would be for the optimizer to simply search the history of sampled points to find the point with the highest y-value.  However, a different method would be to use the Gaussian process to find the point with the best *predicted* y-value.  We'll add to the `GaussianProcessOptimizer` class a function which takes a `expected` argument that allows the user to choose between the two options.  When `expected=True` (the default), the opimizer will search for the point which has the best y-value, as predicted by the Gaussian process.  When `expected=False`, the optimizer object will simply return the point in the sample history which had the highest y-value.

```
def best_point(self, expected=True):
    """Get the best point (expected or actual)"""

    # Return best point which was actually sampled
    if not expected:
        return self.x[self.y.index(min(self.y))]

    # Fit the Gaussian process to samples so far
    self._fit_gp()

    # Find x with greatest expected score
    flip = np.power(-1, self.minimize)
    rx = x.reshape(-1,self.num_dims)
    score_func = lambda x: flip*self._pred_gp(rx)
    x = self.random_point()
    best_score = np.inf
    n_restarts = 10
    for iR in range(n_restarts):

        # Maximize expected improvement
        res = minimize(score_func,
                       self.random_point(),
                       method='L-BFGS-B',
                       bounds=self.bounds)

        # Keep x if it's the best so far
        if res.fun < best_score:
            best_score = res.fun
            x = res.x

    # Return x with highest expected value
    return x
```

Finally, we want a way to visualize the Gaussian process' estimate of the expensive function's value.  We'll add to the optimizer class a method which plots the sampled points along with the Gaussian process' estimate of the function value in either one or two dimensions:

```
def plot_surface(self, x_dim=0, y_dim=None, res=100):
    """Plot the estimated surface of the function"""

    # Fit the Gaussian process to points so far
    self._fit_gp()

    # 1D plot
    if y_dim is None:

        # Predict y as a fn of x
        res = 100 #resolution
        x_pred = np.ones((res, self.num_dims))
        x_pred *= (np.array(self.bounds)
                   .mean(axis=1)
                   .reshape(-1, self.num_dims))
        x_pred[:,x_dim] = np.linspace(self.lb[x_dim],
                                      self.ub[x_dim], res)
        y_pred, y_err = self._pred_gp(x_pred, return_std=True)

        # Plot the Gaussian process' estimate of the function
        plot_err(x_pred[:, x_dim], y_pred, y_err)

        # Plot the sampled points
        for iP in range(len(self.x)):
            plt.plot(self.x[iP][x_dim],
                     self.y[iP], '.', color='0.6')

    # 2D plot
    else:

        # Predict y as a fn of x
        res = 100 #resolution
        x_pred = np.ones((res*res, self.num_dims))
        x_pred *= (np.array(self.bounds)
                   .mean(axis=1)
                   .reshape(-1, self.num_dims))
        xp, yp = np.meshgrid(
            np.linspace(self.lb[x_dim], self.ub[x_dim], res),
            np.linspace(self.lb[y_dim], self.ub[y_dim], res))
        x_pred[:,x_dim] = xp.reshape(-1)
        x_pred[:,y_dim] = yp.reshape(-1)
        y_pred = self._pred_gp(x_pred, return_std=False)

        # Plot the Gaussian process
        plt.imshow(y_pred.reshape((res, res)), aspect='auto',
                   interpolation='bicubic', origin='lower',
                   extent=(self.lb[x_dim], self.ub[x_dim],
                           self.lb[y_dim], self.ub[y_dim]))
        plt.colorbar()

        # Plot the sampled points
        for iP in range(len(self.x)):
            plt.plot(self.x[iP][x_dim],
                     self.x[iP][y_dim],
                     '.', color='0.6')
```

### Minimizing a 1D Function

Now we can find the minimum (or maximum) of any expensive function with our Bayesian optimizer.  Supose we have some expensive function (the function below isn't expensive, of course, but is just for demonstration purposes):

```
# Some 1D function
func = lambda x: -np.sin(x)/x + 0.1*np.random.randn()
```

We can use the optimizer class we just created to find the minimum of that function in as few steps as possible.  The optimizer lets us sample random points in the range:

```
# Create an optimizer object
# (using the class we just created!)
gpo = GaussianProcessOptimizer([-10], [10])

# Randomly sample points
for i in range(10):

    # Get a new random point in the range
    new_params = gpo.random_point()

    # Compute the function value at that point
    # (this is usually the expensive part)
    value = func(new_params)

    # Store the outcome
    gpo.add_point(new_params, value)

    # Plot the function estimate
    gpo.plot_surface()
```

```
#gpo.plot_surface()
# TODO: animation using step arg
```

But more importantly, it provides suggestions as to where to sample next, in order to most effectively minimize the function:

```
# Find the minimum w/ as few steps as possible
for i in range(10):

    # Get a new random point in the range
    new_params = gpo.next_point()

    # Compute the function value at that point
    # (this is usually the expensive part)
    value = func(new_params)

    # Store the outcome
    gpo.add_point(new_params, value)

    # Plot the function estimate
    gpo.plot_surface()
```

    /usr/local/lib/python3.6/dist-packages/sklearn/gaussian_process/gpr.py:480: ConvergenceWarning: fmin_l_bfgs_b terminated abnormally with the  state: {'grad': array([-2.50286965e-05,  1.68961200e-04, -1.21799232e-04]), 'task': b'ABNORMAL_TERMINATION_IN_LNSRCH', 'funcalls': 90, 'nit': 32, 'warnflag': 2}
      ConvergenceWarning)

```
#gpo.plot_surface()
# TODO: animation using step arg
```

Notice how when sampling the points the optimizer suggests, the samples concentrate around the minimum of the function, allowing us to better-explore promising areas of parameter space.

### Minimizing a 2D Function

Our optimizer would only be moderately useful if it could only optimize one parameter at a time!  Luckily Gaussian processes can handle multiple dimensions, and therefore so can our optimizer.  For example, suppose we have an expensive-to-evaluate two-dimensional function (again, this isn't expensive but just for demonstration purposes):

```
# Some 2D function
func = lambda x: (-np.sin(x[0])/x[0]
                  - np.sin(x[1])/x[1]
                  + 0.1*np.random.randn())
```

We can create a 2D optimizer and make some random samples in exactly the same way as before:

```
# Create an optimizer object
gpo = GaussianProcessOptimizer([-10, -10], [10, 10], dtype=[float, float])

# Perform a few random samples
for i in range(40):

    # Get a new random point in the range
    new_params = gpo.random_point()

    # Compute the function value at that point
    value = func(new_params)

    # Store the outcome
    gpo.add_point(new_params, value)

    # Plot the surface
    #gpo.plot_surface(x_dim=0, y_dim=1)
```

```
gpo.plot_surface(x_dim=0, y_dim=1)
plt.grid(False)
plt.show()
```

And we can also make directed samples in the same way as before:

```
# Sample directed by the optimizer
for i in range(10):

    # Point suggested by optimizer
    new_params = gpo.next_point()

    # Compute the function value at that point
    value = func(new_params)

    # Store the outcome
    gpo.add_point(new_params, value)

    # Plot the surface
    #gpo.plot_surface(x_dim=0, y_dim=1)
```

```
gpo.plot_surface(x_dim=0, y_dim=1)
plt.grid(False)
plt.show()
```

Notice how the optimizer now focuses on promising areas, occasionally exploring the more unexplored areas of parameter-space.  Also keep in mind that the `GaussianProcessOptimizer` isn't limited to just 1 or 2 dimensions - in theory it can handle as many dimensions as you can throw at it.  But of course the curse of dimensionality will catch up with you fairly quickly and you'll need an inordinate number of samples to perform effective optimization.

## Hyperparameter Optimization

Now that we have a Bayesian optimizer, we can create a function to find the hyperparameters of a machine learning model which optimize the cross-validated performance.  All this function needs is the x and y data, the predictive model (in the form of an sklearn Estimator), and the hyperparameter bounds.  The user should also be able to specify what metric to optimize (e.g. mean squared error or accuracy), in addition to whether that metric should be maximized or minimized.  This is because some metrics are better when lower, like mean squared error, while others are better when higher, like the coefficient of determination.  Users should also be able to specify how many samples to be taken, and how many random samples to take before switching to Bayesian optimization.  This function won't be a method of `GaussianProcessOptimizer`, but rather a stand-alone function which uses that class.

```
def optimize_params(X, y, model, bounds,
                    metric=make_scorer(mean_squared_error),
                    minimize=True,
                    n_splits=3,
                    shuffle=True,
                    max_evals=50,
                    n_random=5):
    """Optimize model parameters using cross-fold validation.

    Parameters
    ----------
    X : pandas DataFrame
        Independent variable values (features)
    y : pandas Series
        Dependent variable values (target)
    model : sklearn Estimator
        Predictive model to optimize
    bounds : dict
        Parameter bounds.
    metric : sklearn scorer
        Metric to use for evaluation.
    minimize : bool
        Whether to minimize ``metric``.
        If true, minimize; if false, maximize.
    n_splits : int
        Number of cross-validation folds.
    shuffle : bool
        Whether to shuffle samples before splitting into CV folds.
    max_evals : int
        Max number of cross-validation evaluations to perform.
    n_random : int
        Number of evaluations to use random parameter combinations
        before switching to Bayesian global optimization.

    Returns
    -------
    opt_params : dict
        Optimal parameters.
    optimizer : dsutils.optimization.GaussianProcessOptimizer
        Optimizer used to select the points.  Contains the history
        of all points which were sampled.
    """

    # Collect info about parameters to optimize
    Np = len(bounds) #number of parameters
    step_params = [e for e in bounds]
    steps = [e.split('__')[0] for e in step_params]
    params = [e.split('__')[1] for e in step_params]
    lb = [bounds[e][0] for e in step_params]
    ub = [bounds[e][1] for e in step_params]

    # Initialize the Gaussian process optimizer
    gpo = GaussianProcessOptimizer(lb, ub, minimize=minimize)

    # Create a cross-fold generator
    kf = KFold(n_splits=n_splits, shuffle=shuffle)

    # Search for optimal parameters
    for i in range(max_evals):

        # Get next set of parameters to try
        if i < n_random:
            new_params = gpo.random_point()
        else:
            new_params = gpo.next_point()

        # Modify model to use new parameters
        for iP in range(Np):
            tP = {params[iP]: new_params[iP]}
            model.named_steps[steps[iP]].set_params(**tP)

        # Compute and store cross-validated metric
        scores = cross_val_score(model, X, y, cv=kf,
                                 scoring=scorer, n_jobs=n_jobs)

        # Store parameters and scores
        gpo.add_point(new_params, scores)

    # Return optimal parameters and the optimizer object used
    opt_params = dict(zip(step_params, gpo.best_point()))
    return opt_params, gpo
```

### Optimizing One Hyperparameter

Let's try the hyperparameter optimizer out on some real data.

TODO: test on some real data, optimize a hyperparam in 1D, plot score vs hyperparam value

```
# Load the diabetes dataset
# TODO: well, maybe could use one that's a bit bigger...
from sklearn.datasets import load_diabetes
diabetes = load_diabetes()
X_data = pd.DataFrame(diabetes['data'])
y_data = pd.Series(diabetes['target'], X_data.index)
y_data = (y_data-y_data.mean())/y_data.std()
```

```
# CatBoost pipeline
catboost = Pipeline([
    ('scaler', RobustScaler()),
    ('regressor', CatBoostRegressor(verbose=False,
                                    loss_function='RMSE',
                                    learning_rate=0.1))
])

# Parameter bounds
bounds = {
    'regressor__learning_rate': [0.01, 0.5, float],
}

# Find the optimal parameters
opt_params, gpo = optimize_params(
    X_data, y_data, catboost, bounds,
    n_splits=3, shuffle=False,
    max_evals=2, n_random=2)
```

```
gpo.plot_surface()
```

### Optimizing Multiple Hyperparameters

TODO: test on some real data, optimize the hyperparams in 2D, plot animation for 2D case

```
# CatBoost pipeline
catboost = Pipeline([
    ('scaler', RobustScaler()),
    ('regressor', CatBoostRegressor(verbose=False,
                                    loss_function='RMSE',
                                    learning_rate=0.1,
                                    l2_leaf_reg=2.0))
])

# Parameter bounds
bounds = {
    'regressor__learning_rate': [0.01, 0.5, float],
}
```

## Conclusion

TODO: Other packages for Bayesian optimization

* [HyperOpt](http://hyperopt.github.io/hyperopt)
* [skopt](https://scikit-optimize.github.io)
* [MOE](https://github.com/Yelp/MOE)
* [HyperparameterHunter](https://github.com/HunterMcGushion/hyperparameter_hunter/blob/master/README.md)
* [Spearmint](https://github.com/HIPS/Spearmint)
* [BayesOpt](https://github.com/rmcantin/bayesopt)
* [SMAC](https://automl.github.io/SMAC3/stable)
* [HPOlib2](https://automl.github.io/HPOlib2/stable)

TODO: Also mention that tree parzen estimators have been used instead of gaussian processes

max_stars_repo_name: brendanhasz/dsutils

max_issues_repo_name: brendanhasz/dsutils

max_forks_repo_name: brendanhasz/dsutils


text:
# Quantum phase estimation

Your task in this notebook is to implement the quantum Fourier transform on a set of 3 qubits in Qiskit, and then use it to estimate the eigenvalue of a simple Hamiltonian.

### Part A: Implementing the QFT

```python
import numpy as np

from qiskit import QuantumRegister, ClassicalRegister
from qiskit import QuantumCircuit
from qiskit import execute, BasicAer

import qiskit.tools.visualization as qvis

import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning)
```

**Task 1.** In the lecture slides I showed the QFT in general; using this template, write out the circuit for the 3-qubit case. Don't forget the SWAP at the end!

**Task 2.** Recall that the rotations have the form
\begin{equation}
 R_k = \begin{pmatrix}
       1 & 0 \\ 0 & e^{2\pi i / 2^k}
       \end{pmatrix}
\end{equation}
Evaluate the unitary operations for the rotations present in Task 1. Do these gates look familiar?

** Task 3.** Now it's time to implement the circuit. I've set up a 3-qubit quantum register, so all you have to do is implement the gates. For reference, [this page](https://qiskit.org/documentation/terra/summary_of_quantum_operations.html#controlled-operations-on-qubits) contains information on the gates you can implement natively in Qiskit. Hint: look at the `cu1` gates.

```python
# Empty register
q = QuantumRegister(3)

# This object will be our circuit
qft = QuantumCircuit(q)

# Apply gates to the first qubit

# Apply gates to the second qubit

# Apply  remaining gates
```

### Part B: Eigenvalue estimation

In this part, we're going to implement the phase estimation algorithm to get the eigenvalue of the following unitary:
\begin{equation}
 U = \begin{pmatrix}
       1 & 0 \\
       0 & e^{5\pi i / 4} \\
       \end{pmatrix}
\end{equation}

** Task 4. ** What is the eigenvalue $\varphi$ of this matrix as expressed in the phase estimation problem, i.e. $e^{2\pi i \varphi}$? What is its expansion in the notation $0.\varphi_1 \cdots \varphi_n$? How many qubits do you need to represent this exactly?

Now let's set up the two quantum registers: one for the qubits to represent phase, and 1 for the eigenvector since it is a 2-dimensional vector. We'll also need some classical bits to hold the measurement outcomes of the phase register.

```python
# Change this to the number of qubits needed in your upper register
n_phase_qubits = -1

ph_reg = QuantumRegister(n_phase_qubits)
eig_reg = QuantumRegister(1)
meas_reg = ClassicalRegister(n_phase_qubits)

qpe = QuantumCircuit(ph_reg, eig_reg, meas_reg)
```

** Task 5. ** Apply the phase estimation algorithm to your registers. Recall that you'll also need to initialize the eigenvector register into the proper eigenvector.

```python
# Initialize the eigenvector in the lower register
#

# Perform the phase estimation; we're applying a controlled U varying amounts of times
#
```

** Task 6. ** Apply the inverse QFT to your phase register.

Note: unfortunately, Qiskit does not yet support a simple circuit reversal operation; so you'll have to manually list the inverses of each gate. Don't forget to take the adjoint!

```python
# Apply the inverse QFT to the phase register of the qpe circuit
```

** Final task. ** Let's simulate the circuit and measure to get the vector corresponding to our phase!

```python
qpe.measure(ph_reg, meas_reg)

backend = BasicAer.get_backend('qasm_simulator') # the device to run on
result = execute(qpe, backend, shots=1000).result()
counts = result.get_counts(qpe)

eigenvalue = 0

# Your code here: Extract the counts and compute the decimal value
# Recall that Qiskit measurement results come out 'backwards', so make sure
# to put the output computational basis state in the reverse order first

print(f"Phase estimation routine returned the eigenvalue {eigenvalue}.")
```

max_stars_repo_name: a-capra/Intro-QC-TRIUMF

max_issues_repo_name: a-capra/Intro-QC-TRIUMF

max_forks_repo_name: a-capra/Intro-QC-TRIUMF


text:
# Intro

This Notebook documents a few investigations for modelling point and disk loads using analytic solutions to the equations of elastic equilibrium. These solutions are explained in detail in Farrell 1972 and also Chapter 3 of the book Contact Mechanics and given briefly here.

## Disk Load

For a disk load of unit weight (force), the solution for vertical displacement from Farrell 72 is:

\begin{equation}
u(0,r) =
    \begin{cases}
    -\frac{\sigma}{\pi^2\mu\eta \alpha} E\big(\frac{r}{
    \alpha} \big) & r < \alpha \\
    -\frac{\sigma r}{\pi^2\mu\eta \alpha^2} \Big\lbrack E\big(\frac{\alpha}{
    r}\big) - \big( 1 - \frac{\alpha^2}{r^2} K\big(\frac{\alpha}{r}\big) \Big\rbrack & r > \alpha
    \end{cases}
\end{equation}

where $E(x)$ and $K(x)$ are the complete elliptic integrals, which are special functions commonly tabulated (and included in scipy). As for constants, $\mu$ is the shear modulus, $\lambda$ is Lame's first parameter, $\eta = \mu + \lambda$ and $\alpha$ is the disk radius.

This scenario is also derived in Contact Mechanics, in terms of Young's Modulus and Poisson's Ratio, but I haven't gone through it in detail. The relationship between YM, Poisson's Ratio and Lame Parameters amongst more is here: ( or see SubSurf wiki if you don't have a magnifying glass to hand http://subsurfwiki.org/wiki/1st_Lam%C3%A9_parameter)

## Point Load

For a point load, the solution from Farrel 72 is now obtainable for all z, not only z=0. It is:

\begin{equation}
u(z,r) = -\frac{1}{4 \pi \mu R} \big( \frac{\sigma}{\eta} + \frac{z^2}{R^2} \big)
\end{equation}

where $R = \sqrt{r^2 + z^2}$.

## Argus et al

Argus et al 2014 and 2017 use a more complicated formulation for a radially-stratified spherical Earth. They may even take into account graviational effects. However, we could expect the simple analytic solution detailed above to agree broadly with theirs. In their paper, they give two specific examples: loading from two 0.62 GT loads, one with a radius of 7 km and the other with a radius of 14 km.

# Reproducing Argus et al 2014

I use the analytic model for a disk with the same radius and mass as Argus, and reproduce similar results to theirs. Note that to plot their result, I have picked results off their graph by hand.

```python
%matplotlib inline
```

```python
import numpy as np
from scipy import special
import matplotlib.pyplot as plt
import seaborn as sns

g = 9.81
lamda_gran = 25*10**9 # Lame parameter of granite
mew_gran = 40*10**9 # Shear modulus of granite

def calc_deformation_disk(M,lamda,mew,alpha,r):
    '''
    Calculates the deformation for a disk of mass M and radius alpha, on top of a substrate with Lame parameter lambda and shear modulus mew.
    M: mass of the disk, kg
    lamda: Lame parameter of substrate
    mew: shear modulus of substrate
    alpha: disk radius, in metres
    r: array of distances from centre of disk at which to calculate solution. In metres. eg r = np.linspace(0,50*10**3,num=1000) to go to 50km distance.

    '''
    sigma=lamda+2*mew
    nabla=lamda + mew

    defm=np.zeros_like(r)

    r_disk = r[r<=alpha]
    r_postdisk = r[r>=alpha]

    defm[r<=alpha]=-M*g* (sigma/(np.pi**2 * mew * nabla * alpha) * special.ellipe((r_disk/alpha)**2) )

    defm[r>=alpha]= -M *g* (sigma* r_postdisk / (np.pi**2 * mew * nabla * alpha**2)) * (special.ellipe((alpha/r_postdisk)**2) - (1 - (alpha/r_postdisk)**2) * special.ellipk((alpha/r_postdisk)**2))

    return defm

argus_traced_14_disp = -1*np.array([2.6, 2.3, 1.4, 0.7, 0.3, 0.15]) # These are the displacements which I picked by hand from Argus 2017
argus_traced_14_distance = [0,9,16,20,36,50] # And these are the corresponding distances in kilometres

M_Argus1 = np.pi * (14*10**3)**2 * 1 * 1000
#M_Argus2 = np.pi * (7*10**3)**2 * 4 * 1000
alpha_Argus1 = 14*10**3
#alpha_Argus2 = 7*10**3

r_argus = np.linspace(0,50*10**3,num=1000)

defm_Argus1 = calc_deformation_disk(M_Argus1,lamda_gran,mew_gran,alpha_Argus1,r_argus)
#defm_Argus2 = calc_deformation_disk(M_Argus2,lamda_gran,mew_gran,alpha_Argus2,r_argus)
sns.set_context('poster')

fig,ax = plt.subplots(1,1,figsize=(18,12))
plt.plot(r_argus/1000,defm_Argus1*1000,color='blue',label="0.62 Gt; Radius = 14km")
plt.plot(argus_traced_14_distance,argus_traced_14_disp,'--',color='blue',label="Traced 14km from Argus et al")
#plt.plot(r_argus/1000,defm_Argus2*1000,color='orange',label="0.62 Gt; Radius = 7km")
#plt.plot(point_Argusr/1000,pointloadArgus1*1000,label='Point Load Argus1')

plt.xlabel('Distance/km')
plt.ylabel('Displacement after loading/mm')

plt.title("Reproducing Argus results")

plt.legend()
plt.show()

```

# How deep are we sensitive to: case of point loading

BELOW HERE I HAVEN'T TIDIED UP; GENERALLY UNLABELLED AND TO BE LOOKED AT AT YOUR OWN RISK!

```python

#
#
#defm_PineFlat = calc_deformation(M_Pine,lamda_gran,mew_gran,alpha,r)
#
#plt.plot(r/1000,defm_granite*1000,label="Granite Substrate")
#plt.plot(r/1000,defm_sandstone*1000,label="Sandstone Substrate")
#
#
#plt.xlabel('Distance/km')
#plt.ylabel('Displacement after loading/mm')
#
#plt.title("Modelled subsidence for Lake Kaweah")
#plt.legend()

#%%
sns.set_context('notebook')

def calc_deformation_point(R,Z,lamda,mew,M):
    sigma = lamda + 2*mew
    nabla = lamda + mew
    RR = np.sqrt(Z**2 + R**2)
    disp = -M/(4*np.pi*mew*RR) * (sigma/nabla + Z**2/RR**2)
    return disp

M_kaweah = 1.968*10**(11)
z = np.linspace(1,4000,num=200)
r_2d = np.linspace(1,4000,num=200)
R,Z = np.meshgrid(r_2d,z)
disp = calc_deformation_point(R,Z,lamda_gran,mew_gran,M_kaweah)

fig2= plt.figure()
plt.pcolormesh(R/1000,Z/1000,disp*1000,vmax=0,vmin=-5)
cbar = plt.colorbar()
cbar.set_label('Vertical displacement / mm')

plt.gca().invert_yaxis()
plt.xlabel('Radial Distance r / km')
plt.ylabel('z / km')

#%% Put the point load on the original plot

sns.set_context('poster')

pointload=disp[2,:]
point_r = R[2,:]

#ax.plot(point_r/1000,pointload*1000,label='Point Load')
fig.legend()

```

# How Important is Substrate

```python
import numpy as np
from scipy import special
import matplotlib.pyplot as plt
import seaborn as sns

g = 9.81

M_kaweah = 1.968*10**(11)
lamda_gran = 25*10**9
lamda_sandstone = 15*10**9
lamda_uncons = 1*10**9

mew_gran = 40*10**9
mew_sandstone=30*10**9
mew_uncons = 1*10**9
r = np.linspace(0,10*10**3,num=100)
alpha = 1.25*10**3

def calc_deformation(M,lamda,mew,alpha,r):
    sigma=lamda+2*mew
    nabla=lamda + mew

    defm=np.zeros_like(r)

    r_disk = r[r<=alpha]
    r_postdisk = r[r>=alpha]

    defm[r<=alpha]=-M*g* (sigma/(np.pi**2 * mew * nabla * alpha) * special.ellipe((r_disk/alpha)**2) )

    defm[r>=alpha]= -M *g* (sigma* r_postdisk / (np.pi**2 * mew * nabla * alpha**2)) * (special.ellipe((alpha/r_postdisk)**2) - (1 - (alpha/r_postdisk)**2) * special.ellipk((alpha/r_postdisk)**2))

    return defm

#%%
sns.set_context('poster')

fig,ax = plt.subplots(1,1,figsize=(18,12))

defm_granite = calc_deformation(M_kaweah,lamda_gran,mew_gran,alpha,r)
defm_sandstone = calc_deformation(M_kaweah,lamda_sandstone,mew_sandstone,alpha,r)
defm_uncons = calc_deformation(M_kaweah,lamda_uncons,mew_uncons,alpha,r)

plt.plot(r/1000,defm_granite*1000,label="Granite Substrate")
plt.plot(r/1000,defm_sandstone*1000,label="Sandstone Substrate")
#plt.plot(r/1000,defm_uncons*1000,label="Unconsolidated Seds Substrate (??)")

plt.xlabel('Distance/km')
plt.ylabel('Displacement after loading/mm')

plt.title("Modelled subsidence for Lake Kaweah")

plt.legend()

```

```python
#%% Check that I am using the right Lame parameters

Vp = 8000
Vs = 4200
rho = 2600

lamda_granite_Check = rho * (Vp**2 - 2*Vs**2)
mew_granite_Check = rho * Vs**2

```

```python
from IPython.display import HTML
HTML('''
The raw code for this IPython notebook is by default hidden for easier reading.
To toggle on/off the raw code, click <a href="javascript:code_toggle()">here</a>.''')

```

The raw code for this IPython notebook is by default hidden for easier reading.
To toggle on/off the raw code, click here.

```python
def calc_deformation_2Ddisk(M,lamda,mew,alpha,R):
    ''' Calculates the deformation for a disk of mass M and radius alpha, on top of a substrate with Lame parameter lambda and shear modulus mew. R should be a meshgrid(r).'''
    sigma=lamda+2*mew
    nabla=lamda + mew

    defm=np.zeros_like(R)

    r_disk = r[r<=alpha]
    r_postdisk = r[r>=alpha]

    defm[r<=alpha]=-M*g* (sigma/(np.pi**2 * mew * nabla * alpha) * special.ellipe((r_disk/alpha)**2) )

    defm[r>=alpha]= -M *g* (sigma* r_postdisk / (np.pi**2 * mew * nabla * alpha**2)) * (special.ellipe((alpha/r_postdisk)**2) - (1 - (alpha/r_postdisk)**2) * special.ellipk((alpha/r_postdisk)**2))

    return defm

```

max_stars_repo_name: envgp/SnowInv

max_issues_repo_name: envgp/SnowInv

max_forks_repo_name: envgp/SnowInv


text:
# 15.1. Diving into symbolic computing with SymPy

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

```python
var('x y')
```

```python
x, y = symbols('x y')
```

```python
expr1 = (x + 1) ** 2
expr2 = x**2 + 2 * x + 1
```

```python
expr1 == expr2
```

```python
simplify(expr1 - expr2)
```

```python
expr1.subs(x, expr1)
```

```python
expr1.subs(x, pi)
```

```python
expr1.subs(x, S(1) / 2)
```

```python
_.evalf()
```

```python
f = lambdify(x, expr1)
```

```python
import numpy as np
f(np.linspace(-2., 2., 5))
```

max_stars_repo_name: jhgoebbert/jupyter-jsc-notebooks

max_issues_repo_name: jhgoebbert/jupyter-jsc-notebooks

max_forks_repo_name: jhgoebbert/jupyter-jsc-notebooks