text:
```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use('classic')
%matplotlib inline
```

# Class 8: Deterministic Time Series Models II

## Nonlinear First-Order Difference Equations

Recall the Solow growth model with no exogenous productivity growth (written in "per worker" terms):

\begin{align}
y_t & = Ak_t^{\alpha}\\
k_{t+1} & = i_t + (1-\delta)k_{t}\\
y_t & = c_t + i_t\\
c_t & = (1-s)y_t,
\end{align}

where $y_t$ is output per worker, $k_t$ is capital per worker, $c_t$ is consumption per worker, and $i_t$ is investment per worker. We've seen that by treating investment as an exogenous quantity, then the capital accumulation equation could be viewed as a linear first-order difference equation. But in the Solow model, investment is *not* exogenous and is equal to the savings rate times output per worker:

\begin{align}
i_t & = sy_t = sAk_t^{\alpha}
\end{align}

Therefore the equilibrium law of motion for capital per worker is:

\begin{align}
k_{t+1} & = sAk_t^{\alpha} + (1-\delta)k_{t}, \label{eqn:capital_solution}
\end{align}

which is a *nonlinear* first-order difference equation. In equilibrium capital in period $t+1$ is a *concave* (i.e., increasing at a decreasing rate) function of capital in period $t$. We can iterate on the nonlinear law of motion just like we iterated on the linear difference equation.

Let's suppose the following values for the simulation:

| $$k_0$$ | $$s$$ | $$A$$  | $$\alpha$$ | $$\delta$$  | $$T$$  |
|---------|-------|--------|------------|-------------|--------|
| 10      | 0.1   | 10     | 0.35       |  0.1        | 101    |

Where $T$ is the total number of simulation periods (i.e., $t$ will range from $0$ to $100$).

```python
# Create variables 'k0', 's', 'A', 'alpha','delta','T' to store parameter values for the simulation

# Initialize capital per worker variable 'capital' as an array of zeros with length T

# Set first value of capital equal to k0

# Iterate over t in range(T-1) to update subsequent values in the capital array

# Construct a plot of simulated capital per worker

```

It looks like capital per worker is converging toward a steady state value near $35$. You can verify that the steady state is in fact $34.55$.

So computing nonlinear difference equations is essentially the same as computing linear differrence equations. However, establishing the stability properties of nonlinear difference equations is more involved and we'll skip that. But in case you're interested, I've included an optional discussion at the bottom of this notebook.

## Systems of Difference Equations

Most macroeconomic models are dynamic *systems* of equations. They are collections of several equations that simultaneously determine several variables. Working with them can be a little intimidating, but in principle, it's similar to working with a single dynamic equation.

### Example: Solow Growth Model Without Exogenous Growth

The basic Solow growth model is a system of four equations that determine values of four variables $k_t$, $y_t$, $c_t$, and $i_t$. In the Solow model, capital per worker is *state variable* because the the value of $k_t$ was actually determined in period $t-1$. In fact, in the Solow model without exogenous population of productivity growth, capital is the *only* state variable. That means that if you know the value of capital per worker, then you can compute all of the other quantities. The nonstate variables are called *control* variables.

To make this point clear, here is the Solow growth model *solved* in terms of $k_t$:

\begin{align}
y_t & = Ak_t^{\alpha}\\
i_t & = sAk_t^{\alpha}\\
c_t & = (1-s)Ak_t^{\alpha}\\
k_{t+1} & = sAk_t^{\alpha} + (1-\delta)k_{t}.
\end{align}

To simulate all of the endogenous variables, do these two steps:

1. Simulate an array of values for $k_t$
2. Compute the other variables using the array of $k_t$ values.

Since we've already simulated capital, we can just compute simulated values for $y_t$, $c_t$, and $i_t$ right away.

```python
# Create variables 'output', 'consumption', and 'investment' equal to the respective variables in the model

# Construct a plot of simulated output per worker, consumption per worker, investment per worker, capital per worker.

# Create legend. Use ax.legend(loc='center left', bbox_to_anchor=(1, 0.5)) to put legend outside of axis

# Add grid if you want

```

Of course, if we want to compute multiple simulations, then we will want to write a function to compute the simulated values so that we don't have to repeat ourselves. The function in the following cell returns a `DataFrame` with columns containing simulated values of capital per worker, output per worker, consumption per worker, and investment per worker.

```python
# Define a function that returns a DataFrame of simulated value from the Solow model with no exogenous growth. CELL PROVIDED
def solow(s,A,alpha,delta,T,k0):

    '''Function for computing a simulation of the Solow growth model without exogenous growth.

    The model is assmed to be in "per-worker" quantities:

            y[t]   = A*k[t]^alpha
            k[t+1] = i[t] + (1-delta)*k[t]
            y[t]   = c[t] + i[t]
            c[t]   = (1-s)*y[t]

    Args:
        s (float):      Saving rate
        A (float):      TFP
        alpha (float):  Capital share in Cobb-Dougla production function
        delta (float):  Capital depreciation rate
        T (int):        Number of periods to simulate
        k0 (float):     Initial value of capital per worker

     Returns:
         Pandas DataFrame
     '''

    # Initialize capital per worker values
    capital = np.zeros(T)

    # Set first value of capital per worker equal to k0
    capital[0] = k0

    # Iterate over t in range(T-1) to update subsequent values in the capital array
    for t in range(T-1):
        capital[t+1] = s*A*capital[t]**alpha + (1-delta)*capital[t]

    # Compute the values of the other variables
    output = A*capital**alpha
    consumption = (1-s)*output
    investment = s*output

    # Put simulated data into a DataFrame
    df = pd.DataFrame({'output':output,'consumption':consumption,'investment':investment,'capital':capital})

    # Return the simulated data
    return df
```

Use the function `solow()` to simulate the Solow model with the same parameters that we've been using.

```python
# Simulate the model and store results in a variable called 'df'

# Print the first five rows of df

```

Now, use the function `solow()` to simulte the model for five different initial values of capital per worker: $k_0 = 10, 20, 30, 40, 50$. Plot the trajectories for $y_t$ together.

```python
# Create variable 'initial_ks' that store the five initial values of k[t]

# Create a figure and axis

# Iterate over k0 in initial_ks and plot

# Add title

# Create legend. Use ax.legend(loc='center left', bbox_to_anchor=(1, 0.5)) to put legend outside of axis

# Add grid if you want

```

### Solow Model with Exogenous Labor

Now consider the Solow growth model with exogenous labor growth. That is, let $L_t$ denote the quantity of labor and suppose that the population (and therefore the labor supply) grows at the constant rate $n$. This means:

\begin{align}
L_{t+1} & = (1+n)L_t,
\end{align}

where the intial value of labor $L_0$ must be given. The rest of the model in aggregate terms is:

\begin{align}
Y_t & = AK_t^{\alpha}L_t^{1-\alpha}\\
K_{t+1} & = I_t + (1-\delta)K_{t}\\
Y_t & = C_t + I_t\\
C_t & = (1-s)Y_t,
\end{align}

where the intial value of capital $K_0$ is also given.

Since capital *and* labor are both determined in the previous period by their respective laws of motion, the Solow model with exogenous labor growth has *two* state variables: $K_t$ and $L_t$. We can solve the capital law of motion for capital in terms of only the variables capital and labor and so the two equations that determine how the state of the economy evolves are:

\begin{align}
K_{t+1} & = sAK_t^{\alpha}L_t^{1-\alpha} + (1-\delta)K_{t}\\
L_{t+1} & = (1+n)L_t
\end{align}

If we iterate on these two to compute simulations of $K_t$ and $L_t$, then we can compute $Y_t$, $C_t$, and $I_t$ easily.

Let's suppose the following values for a simulation:

| $$L_0$$ | $$n$$  | $$K_0$$ | $$s$$  | $$A$$  | $$\alpha$$ | $$\delta$$ | $$T$$ |
|---------|--------|---------|--------|--------|------------|------------|-------|
| 1       | 0.01   | 10      | 0.1    | 10     | 0.35       | 0.1        | 101   |

Compute simulated values for $K_t$ and $L_t$ and use those values to compute and plot simulated series for output $Y_t$ and output *per worker* $Y_t/L_t$ side-by-side in the same figure.

```python
# Create variables 'K0', 'L0', 'n', 's', 'A', 'alpha','delta','T' to store parameter values for the simulation

# Initialize capital variable 'capital' as an array of zeros with length T

# Set first value of capital equal to K0

# Initialize the labor variable 'labor' as an array of zeros with length T

# Set first value of capital equal to L0

# Iterate over t in range(T-1) to update subsequent values in the capital and labor arrays

# Create variables 'output' and 'output_pw' that store output and output per worker

# Create figure

# Construct a plot of simulated output

# Construct a plot of simulated output per worker

```

Because of exogenous labor growth, there is long-run growth in output but not output per worker. The function in the cell below simulates the Solow model with exogenous labor growth.

```python
# Define a function that returns a DataFrame of simulated value from the Solow model with exogenous labor growth. CELL PROVIDED
def solow_w_L(s,A,alpha,delta,n,T,K0,L0):

    '''Function for computing a simulation of the Solow growth model with exogenous labor.

            Y[t]   = A*K[t]^alpha*L[t]^(1-alpha)
            K[t+1] = I[t] + (1-delta)*K[t]
            Y[t]   = C[t] + I[t]
            C[t]   = (1-s)*Y[t]
            L[t+1] = (1+n)*L[t]

    Args:
        s (float):      Saving rate
        A (float):      TFP
        alpha (float):  Capital share in Cobb-Douglas production function
        delta (float):  Capital depreciation rate
        n (float):      Labor growth rate
        T (int):        Number of periods to simulate
        k0 (float):     Initial value of capital per worker
        L0 (float):     Initial labor

     Returns:
         Pandas DataFrame
     '''

    # Initialize capital values
    capital = np.zeros(T)

    # Set first value of capital equal to K0
    capital[0] = K0

    # Initialize labor values
    labor = np.zeros(T)

    # Set first value of labor equal to L0
    labor[0] = L0

    # Iterate over t in range(T-1) to update subsequent values in the capital and labor arrays
    for t in range(T-1):
        capital[t+1] = s*A*capital[t]**alpha*labor[t]**(1-alpha) + (1-delta)*capital[t]
        labor[t+1] = (1+n)*labor[t]

    # Compute the values of the other aggregate variables
    output = A*capital**alpha*labor**(1-alpha)
    consumption = (1-s)*output
    investment = s*output

    # Compute the values of the "per worker" variables
    capital_pw = capital/labor
    output_pw = output/labor
    consumption_pw = consumption/labor
    investment_pw = investment/labor

    # Put simulated data into a DataFrame
    df = pd.DataFrame({'output':output,
                       'consumption':consumption,
                       'investment':investment,
                       'capital':capital,
                       'output_pw':output_pw,
                       'consumption_pw':consumption_pw,
                       'investment_pw':investment_pw,
                       'capital_pw':capital_pw})

    # Return the simulated data
    return df
```

Use the function `solow_w_L()` to replicate the previous simulation and plots.

```python
# Replicate the previous simulation exercise using the solow_no_exo_growth() function. CELL PROVIDED

# Simulate the model
df = solow_w_L(s,A,alpha,delta,n,T,K0,L0)

# Create a figure
fig = plt.figure(figsize=(12,4))

# Construct a plot of simulated output
ax = fig.add_subplot(1,2,1)
ax.plot(df['output'],lw=3,alpha=0.75)
ax.set_title('Ouput')
ax.grid()

# Construct a plot of simulated output per worker
ax = fig.add_subplot(1,2,2)
ax.plot(df['output_pw'],lw=3,alpha=0.75)
ax.set_title('Ouput per Worker')
ax.grid()
```

## Stability of Nonlinear First-Order Difference Equations (Optional)

Unlike the linear first-order difference equation, stability (i.e., whether the process diverges to infinity in absolute value) of nonlinear first-order difference equations is less straightforward to establish. In fact, proving or disproving *global stability* -- Will the process remain finite for *any* initial condition? -- is particularly challenging. An easier task is to establish the existience of *local stability*: that the process is stable in a *neighborhood* of a given point. The idea is to use the first-order Taylor series approximation (https://en.wikipedia.org/wiki/Taylor_series) to approximate the nonlinear equation with a linear one. An example using the Solow model will make this easier to explain.

Let $k^*$ denote the steady state of capital per worker in the Solow model with no exogenous growth. Then:

\begin{align}
k^* & = \left(\frac{sA}{\delta}\right)^{\frac{1}{1-\alpha}}
\end{align}

The first-order Taylor series approximation to the capital law of motion around steady state capital per worker is:

\begin{align}
k_{t+1} & \approx \left[ sA\left(k^*\right)^{\alpha} + (1-\delta)k^*\right] + \left[\alpha sA\left(k^*\right)^{\alpha-1} + 1-\delta\right]\left(k_t - k^*\right). \label{eqn:capital_approx}
\end{align}

Equation (\ref{eqn:capital_approx}) is a linear first-order difference equation. As long as the coefficient on $k_t$,

\begin{align}
\alpha sA\left(k^*\right)^{\alpha-1} + 1-\delta,
\end{align}

is less than 1 in absolute value, then the approximated model is stable and we say that $k_t$ is stable in a neighborhood of the steady state $k^*$. Let's compute this coefficient for the Solow model using the same parameterization we used earlier:

| $$s$$ | $$A$$  | $$\alpha$$ | $$\delta$$ |
|-------|--------|------------|------------|
| 0.1   | 10     | 0.35       |  0.1       |

```python
# Parameters

# Computer steady state

# Compute coefficient

```

So steady state capital per worker is about 35 and the coefficient on $k_t$ in the approximation is 0.935. So the Solow model as we have parameterized it is stable in the neighborhood of $k^*$. Can you find the upper bound on $\alpha$ for which the model will *not* be stable given $A = 10$, $s = 0.1$ and $\delta = 0.1$?

An obvious question is: How good is the linear approximation? It turns out that for the capital law of motion in the Solow model, the linear approximation is really good. Let's plot the exact law of motion and the approximated law of motion for $k_t \in[0,100]$.

```python
# Create array of values for k[t]

# Compute k[t+1] exactly using the law of motion for capital per worker

# Compute k[t+1] approximately using the approximation to the law of motion for capital per worker

# Create figure

# Plot exact and approximate k[t+1] against for k[t]

```

The exact law of motion and the approximated law of motion intersect at the steady state and are hardly distinguishable elsewhere. You can see the difference between them only by zooming in on the extreme values.

```python
# Create figure

# Plot exact and approximate k[t+1] against for k[t] between 0 and 10

# Plot exact and approximate k[t+1] against for k[t] between 80 and 100

```

For $k_t$ near zero, the approximated law of motion is greater than the exact one but not by much. For $k_t$ near 100, the two laws of motion are still really close to each other. The point is that we can be confident that in the Solow growth model, capital will tend to converge toward the steady state regardless of its initial value.


text:
#
Representación de funciones

```python
from sympy import *
init_printing(use_latex='mathjax')
x, y, z = symbols('x, y, z')
```

El comando básico para dibujar funciones es **plot**. Si solamente le pasamos como argumento una función nos la dibuja en el intervalo (-10,10) del eje $x$. El tamaño del eje $y$ dependerá de la función. Para variar el dominio de la gráfica se le debe añadir como argumento una tupla, que conste de la letra $x$, el límite inferior y el superior.

###
Dibuja la gráfica de la función $f(x)=x^2$ en el intervalo $(-3,4)$.

```python
plot(x**2,(x,-3,4))
```

A la gráfica de podemos hacer muchas modificaciones, añadiendo argumentos optativos. Algunos de ellos son:

* **ylim**. Para modificar el eje $y$.

* **line_color**. Para cambiar el color.

* **xlabel**. Para poner nombre al eje $x$

* **title**. Para poner un título a la gráfica.

###
Varia alguno de los parámetros mencionados.

```python

```

###
Modificar el tamaño y guardar la imagen.

Se pueden dibujar varias funciones en la misma gráfica. La manera más sencilla, cuando ambas gráficas tienen el mismo dominio, es escribir las distintas funciones separadas por comas. Sin embargo para controlar individualmente cada función es mejor utilizar otro método, que consiste en realizar "una tupla de tuplas".

###
Dibuja,  en unos mismos ejes, las gráficas de una parábola  y del coseno.

```python
plot((x**2,(x,-3,3)), (cos(x),(x,-5,5)))
```

Con lo visto anteriormente se puede "simular" la gráfica de una función definida a trozos.

###
Representa la siguiente función definida a trozos:

$$
f(x)=
\begin{cases}
2x-3 & \mathrm{ si } & -4<x<2\\
x^2 & \mathrm{ si } & -2\leq x <2 \\
x & \mathrm{ si }& 2\leq x <5
\end{cases}
$$
rrrplot((2*x-3,(x,-4,-2)),(x**2,(x,-2,2)),(x,(x,2,5)))

```python

```

ext: ipynb

max_stars_repo_path: notebooks/07.- Representacion de funciones.ipynb

max_stars_repo_name: disoftw/python-for-maths

max_issues_repo_name: disoftw/python-for-maths

max_forks_repo_name: disoftw/python-for-maths


text:
```python
from IPython.display import HTML
import numpy as np
import matplotlib.pyplot as plt
import math as m

## . . I've used python symbolic package 'sympy' to complete the tests
## . . You may need to install this package using:
## . .         conda install sympy
from sympy import plot_implicit, symbols, Eq, And
import sympy.plotting as syp

## . . Turn off warnings
import warnings
warnings.filterwarnings("ignore")

HTML('''
The raw code for this Jupyter notebook is by default hidden for easier reading.
To toggle on/off the raw code, click <a href="javascript:code_toggle()">here</a>.
<style>
.output_png {
    display: table-cell;
    text-align: center;
    vertical-align: middle;
}
</style>
''')
```

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

###
Module 3: Complex Numbers

This module focuses on some of the fundamental mathematical skills that are required for the material in this course. Starting from elementary arithmatic, we will refresh our memory of the complex plane; complex exponentials, and multiplication/division and roots; and finally complex functions and the concept of analycity.

## Motivation for Being in the Complex Plane

In many instances in geophysical digitial signal processing we apply [integral transforms](https://en.wikipedia.org/wiki/Integral_transform) that project our geophysical data from one domain to another.  The reason why we do this is that many of the problems we are trying to solve (e.g., filtering out unwanted components, signal-noise separation, etc) are easier to solve in the transform domain. Ideally, the integral transform is **invertible** in the sense that it has an **inverse transform** that allows you to bring your (filtered, separated, etc) data back to the original domain without losing any of its information content.

### Forward and inverse integral transforms

Mathematically, we may write the **integral transform** as the following:

$$(Tg)(u) = \int_{t_1}^{t_2} g(t)\, K(t,u) \,\mathrm{d}t, \tag{1}$$

where the input is a function $g(t)$ of variable $t$, and the output is another function $Tg$ of transform variable $u$. Often we symbolize the output of this transform as $\widehat{G}(u)$.  There are numerous useful integral transforms, where each is specified by a choice of the [integration kernel function](https://en.wikipedia.org/wiki/Integral_transform) $K(u,t)$ of two variables - $t$ and $u$ - as well as the limits applied - $t_1$ and $t_2$.

Some kernels have an associated inverse kernel $K^{-1}(t,u)$, which roughly speaking, delivers an **inverse transform**:

$$ g(t) =  \int_{u_1}^{u_2}  (Tg)(u) \, K^{-1}(t,u) \,\mathrm{d}u, \tag{2} $$

when applied with the appropriate limits $u_1$ and $u_2$. There are many [different integral transforms](https://en.wikipedia.org/wiki/Integral_transform#Table_of_transforms) that are used across a whole range of the physical sciences and mathematics.

### Fourier and Inverse Fourier Transform Pair

The most important set of integral transforms that you will likely experience over your career is the Forward and Inverse Fourier Transform pair (denoted $\mathcal{F}$ and $\mathcal{F}^{-1}$).  By defining $u=f$, we see that the kernel functions for these transforms, $K_{\mathcal{F}}(f,t)$ and $K^{-1}_{\mathcal{F}}(f,t)$, are given by:

$$K_{\mathcal{F}}(f,t) = e^{-i2\pi ft} \quad \mathrm{and} \quad K^{-1}_{\mathcal{F}}(f,t)=e^{i 2\pi f t}. \tag{3}$$

In this course, we will sometimes use the **angular frequency** definition $\omega = 2\pi f$ so that we may write

$$K_{\mathcal{F}}(\omega,t) = e^{-i \omega t} \quad \mathrm{and} \quad K^{-1}_{\mathcal{F}}(\omega,t)=e^{i\omega t}. \tag{4}$$

When we include the limits $t_1=u_1=-\infty$ and $t_2=u_2=\infty$ we may write:

$$\widehat{G}(\omega)\equiv (\mathcal{F}g)(\omega) = \int_{-\infty}^{\infty} g(t) \, e^{-i \omega t} \, \mathrm{d}t \quad \mathrm{and} \quad g(t) =  \int_{-\infty}^{\infty}  (\mathcal{F}g)(\omega) \, e^{i\omega t} \,\mathrm{d}\omega =  \int_{-\infty}^{\infty} \widehat{G}(\omega) \, e^{i\omega t} \,\mathrm{d}\omega. \tag{5}$$

The above Fourier Transform process looks pretty straightforward to compute.  However, there are a number of fundamental questions that need to be addressed:

(1) Does the formulation work for all $g(t)$?
   * **NO!** There are some restrictions on the types of $g(t)$ to which we can apply the Fourier Transform. For example, does it have "finite" (good!) or "infinite" (bad!) energy in the signal?

(2) Is it possible to specify some $g(t)$ for which the Fourier Transform "fails"?
   * **YES!** It is actually very easy to specify a function that would "fail" - e.g., $g(t) = e^{-t} (1-t)^{-1}$ around $t=1$.  This is quite common when we are using designer functions (i.e., "filters") to do specific filtering tasks.

(3) Are we able to analyze these types of signals?
   * **YES!** Just because a signal is not well defined on the unit circle in the complex plane does not mean we have to give up.  We have to expand our ability to analyze signals in complex plane beyond the unit circle.  Doing so leads to some course material that requires us to be comfortable with complex numbers and some of the fundamental principles of complex analysis.  Hopefully, you have encountered these topics in your previous courses and just need a short review.

# Complex Numbers: A Review

You have probably encountered the expression $i=\sqrt{-1}$ with the imaginary unit many times before - likely first with certain solutions to the quadratic equation:

$$ax^2+bx+c=0, \tag{6}$$

given by

$$x=\frac{-b \pm \sqrt{b^2-4ac}}{2a}. \tag{7}$$

When $b^2 < 4ac$ and we can find no real number equal to x.  The solutions in this case are termed **complex numbers**.  We have also assumed the existence of certain arithmetic operations on complex number that largely parallel our experience with real numbers (generally due to some sort of *mathematical symmetry*). However, we do not need to delve into such a detailed discussion for this course.

### Elementary Arithmetic

We first consider simple powers of the imaginary unit $i=\sqrt{-1}$:

* $i^0 = 1$
* $i^1 = i$
* $i^2 = \sqrt{-1}\sqrt{-1} = -1$
* $i^3 = i^2i=-i$
* $i^4 = i^2 i^2=(-1)(-1)=1$
* $i^5 = i$

Note: squaring $i$ causes an effective positive rotation of the output. Applying a square root or other fractional power will naturally cause a negative rotation of the output.

When performing the operations of addition, subtraction, multiplication and division we treat $i$ as though it were an algebraic constant.  For example, if $z_1=3-2i$ and $z_2=2-i$ then:

$$\begin{eqnarray}
{\rm Addition:}\quad z_1+z_2 = (3+2i)+(2-i) &=& 5+i \tag{9a}\\
{\rm Subtraction:}\quad  z_1-z_2=(3+2i)-(2-i) &=& 1+3i \tag{9b}\\
{\rm Multiplication:}\quad  z_1z_2=(3+2i)(2-i)  &=& 6-3i+4i-2i^2 = 6+i+2 = 8+i\tag{9c}\\
\end{eqnarray}$$

Division requires remembering a simple trick.  To compute $\frac{z_1}{z_2}=\frac{3+2i}{2-i}$ we first multiply both the numerator and denominator by the **complex conjugate** (i.e., if $z=a+ib$ then $\overline{z}=\overline{a+ib} = a-ib$ where $[a,b]\in \Re$ and the overbar indicates a complex conjugate) of the term in the denominator.  Thus,

$${\rm Division:}\quad \frac{z_1}{z_2} = \frac{3+2i}{2-i} = \frac{(3+2i)(2+i)}{(2-i)(2+i)} = \frac{4+7i}{5}. \tag{10}$$

This process is also called **rationalizing the denominator**.

We often want to extract out just the real or imaginary components of a complex number, which is written notationally using the **real and imaginary operators**, $\Re$ and $\Im$, respectively. Thus, if $z=a+ib$ we have the following:

$$\Re(z)=a, \tag{10a}$$
and
$$\Im(z)=b. \tag{10b}$$

Note that the result of either operation is a **real number**.

### QUICK PROBLEMS:

(1) Let $z=3-3i$ and $w=1+\sqrt{3}i$.  Find the result of the following operations:

   (a) $z+w$

   (b) $z-w$

   (c) $z\,w$

   (d) $z\,/\,w$

Let's test our answers using the **sympy** package (symbolic python):

```python
# . . Since the variables are not numbers by symbols
# . . I've used python symbolic package 'sympy' to complete the tests
import sympy as sm
z = 3-           3j;
w = 1+sm.sqrt(3)*1j;

# . . Now evaluate as print statements
print('Q1a: z+w=',z+w)
print('Q1b: z-w=',z-w)
print('Q1c: z*w=',z*w)
print('Q1d: z/w=',z/w)
```

    Q1a: z+w= 4.0 - 3.0*I + 1.0*sqrt(3)*I
    Q1b: z-w= 2.0 - 3.0*I - 1.0*sqrt(3)*I
    Q1c: z*w= (1 + 1.0*sqrt(3)*I)*(3.0 - 3.0*I)
    Q1d: z/w= (3.0 - 3.0*I)/(1 + 1.0*sqrt(3)*I)

(2) Let $z=a+ib$ and $w=c+id$, where $[a,b,c,d] \in \mathcal{R}$ (i.e. are real numbers).  Show that:

   (a) $\overline{z+w} = \overline{z}+\overline{w}$

   (b) $\overline{z-w}=\overline{z}-\overline{w}$

   (c) $\overline{z\,w}=\overline{z} \,\overline{w}$

   (d) $\overline{z\,/\,w} = \overline{z} \, / \, \overline{w}$

Let's again test our answers using python:

```python
# . . Let's define all of the required symbols
a,b,c,d,z,w =symbols('a b c d z w')

z = a+1j*b
w = c+1j*d

# . . NOTE: Python uses the '==' operator to test equality
print('Test 2a')
print('z1=conj(z+w)      =',np.conj(z+w))
print('z2=conj(z)+conj(w)=',np.conj(z)+np.conj(w))
print ('Does z1=z2?',np.conj(z+w)==np.conj(z)+np.conj(w))
print('')
print('Test 2b')
print('z1=conj(z-w)      =',np.conj(z-w))
print('z2=conj(z)-conj(w)=',np.conj(z)-np.conj(w))
print ('Does z1=z2?',np.conj(z-w)==np.conj(z)-np.conj(w))
print('')
print('Test 2c')
print('z1=conj(z*w)      =',np.conj(z*w))
print('z2=conj(z)*conj(w)=',np.conj(z)*np.conj(w))
print ('Does z1=z2?',np.conj(z*w)==np.conj(z)*np.conj(w))
print('')
print('Test 2d')
print('z1=conj(z/w)      =',np.conj(z/w))
print('z2=conj(z)/conj(w)=',np.conj(z)/np.conj(w))
print ('Does z1=z2?',np.conj(z/w)==np.conj(z)/np.conj(w))
```

    Test 2a
    z1=conj(z+w)      = conjugate(a) - 1.0*I*conjugate(b) + conjugate(c) - 1.0*I*conjugate(d)
    z2=conj(z)+conj(w)= conjugate(a) - 1.0*I*conjugate(b) + conjugate(c) - 1.0*I*conjugate(d)
    Does z1=z2? True

    Test 2b
    z1=conj(z-w)      = conjugate(a) - 1.0*I*conjugate(b) - conjugate(c) + 1.0*I*conjugate(d)
    z2=conj(z)-conj(w)= conjugate(a) - 1.0*I*conjugate(b) - conjugate(c) + 1.0*I*conjugate(d)
    Does z1=z2? True

    Test 2c
    z1=conj(z*w)      = (conjugate(a) - 1.0*I*conjugate(b))*(conjugate(c) - 1.0*I*conjugate(d))
    z2=conj(z)*conj(w)= (conjugate(a) - 1.0*I*conjugate(b))*(conjugate(c) - 1.0*I*conjugate(d))
    Does z1=z2? True

    Test 2d
    z1=conj(z/w)      = (conjugate(a) - 1.0*I*conjugate(b))/(conjugate(c) - 1.0*I*conjugate(d))
    z2=conj(z)/conj(w)= (conjugate(a) - 1.0*I*conjugate(b))/(conjugate(c) - 1.0*I*conjugate(d))
    Does z1=z2? True

# The Complex Plane

We are familiar with the real line as a geometric device for visualizing real numbers.  Since the complex number $z=x+iy$ is actually two dimensional, it is not surprising that many consider complex numbers as points on a plane.  For example, the point $z=3+4i$ can be represented on the complex plane by:

```python
# . . Express Complex Number
z = 3+4j
a = np.array([z])

# . . Plot unit circle
t = np.linspace(0,2*np.pi,101)
plt.plot(np.cos(t),np.sin(t))

# . . Plot data
plt.plot(a.real, a.imag,'ro')

# . . Make box bigger
plt.axes().set_aspect('equal')
m = max(max(abs(a.real)),max(abs(a.imag)))
plt.xlim(-6,6),plt.ylim(-6,6)

# . . Draw horizontal and vertical axes
plt.plot([0,0],[-6,6],'k',linewidth=0.5)
plt.plot([-6,6],[0,0],'k',linewidth=0.5)

# . . Annotate plot with labels
plt.text(z.real+0.25, z.imag+0.25,'$z$',fontsize=15)
plt.show()
```

**Figure 1. Illustration of point $z=3+4i$ in the complex plane. Note that the blue circle is the unit circle defined by equation $|z|=1$, which I have included for reference.**

## Polar representation of complex numbers

This point can be thought as having Cartesian coordinates $[x,iy]=[3,4]$ or as the vector drawn from the origin to this point.  The complex number $z=3+4i$ can also be defined by **polar coordinate** variables, $r$ and $\theta$, which represent the vector magnitude from the origin of the complex plane and the angle made with respect to the $x$-axis, respectively. Recall that for $z=x+iy$ we have the following forward and inverse coordinate transformations

**Forward (Cartesian to Polar)**

$$r = \sqrt{x^2+y^2}, \tag{11a}$$

$$\theta = {\rm atan}\left(\frac{y}{x}\right). \tag{11b}$$

Note: this is sometimes called atan2 in programming languages and requires two inputs - atan2(y,x) - to put the angle in the correct quadrant.  Using an atan function (if defined) will give you the angle in the first quadrant, which is probably not what you are looking for!

**Inverse (Polar to Cartesian)**

$$x = r \, {\rm cos}\, \theta, \tag{11c}$$

$$y = r \, {\rm sin} \,\theta. \tag{11d}$$

Thus, for the example $z=3+4i$ we have $r=\sqrt{3^2+4^2}=5$ and $\theta=\mathrm{atan}(4/3) = 0.93$ radians.  As in analytic geometry $\theta$ is not unique and we can add to $\theta$ any multiple of $2\pi$ radians:

$$\theta = 0.93+2\pi n,\quad n=0,\pm 1,\pm 2, ... \tag{11e}$$

We often call $r$ the **modulus** of $z$ and denote it by $|z|$.  The angle $\theta$ is called the **argument** of $z$ and is denoted $arg(z)=\theta$ or sometimes just $\angle(z)=\theta$.  It would come as no surprise that complex numbers can be represented by vectors, and that the geometric description of addition and subtraction of vectors is similar to that you would have experienced in other courses.

To help you visualize the unit circle, here is a diagram with various key angles:

**Figure 2. Illustrating various key angles in unit circle.  The values at these angles would be good to memorize for this course.**

### QUICK PROBLEMS

(3) Express the following numbers in the form $re^{i\theta}$ with $r>0$ and $0\le\theta<2\pi$.

* $5i$
* $-7+7i$
* $-4-4\sqrt{3}i$
* $-1$

(4) Express the following numbers in the form $z=x+iy$

* $2e^{i3\pi}$
* $e^{i13\pi/4}$
* $7e^{-i\pi/6}$
* $e^{-i\pi/2}$

## Regions in the Complex Plane

In this course we will sometimes discuss **regions within the complex plane**.  For example:
* (1) a disk contained within the unit circle centered at $c=1+2i$
* (2) the complement of the disk in (1)
* (3) the right half plane
* (4) the upper half plane

Let's again use the sympy library to plot out these regions (though it is difficult to visually see if the sets are open or closed).

```python
## . . Define x and iy symbolically
x,iy = symbols("x iy")

## . . Plot the regions using the Sympy plot_implicit function that can take inequalities
plot_implicit((x-1)**2+(iy-2)**2<1,(x,-5,5),(iy,-3.5,3.5),title='$|z-c|<r$')
plt.show()
```

**Figure 3. Plot showing a disk contained within the unit circle centered at $c=1+2i$.**

```python
plot_implicit((x-1)**2+(iy-2)**2>1,(x,-5,5),(iy,-3.5,3.5),title='Complement of $|z-c|<r$')
plt.show()
```

**Figure 4. Plot showing the complement of the disk shown in Figure 3.**

```python
plot_implicit(x >0,(x,-5,5),(iy,-3.5,3.5),title='Right-Half Plane')
plt.show()
```

**Figure 5. Plot showing the right-half plane.**

```python
plot_implicit(iy>0,(x,-5,5),(iy,-3.5,3.5),title='Upper-half Plane')
plt.show()
```

**Figure 6. Plot showing the upper-half plane.**

Note that there is a difference between an **open set** and a **closed set**.  An open set is a 2D region with the boundary points removed.  For example, the inequality $|z|<r$ would be an open set while $|z|\le r$ would be a closed set.

### QUICK PROBLEMS:

(5) What would the following regions look like:

   * An annulus defined by $r < | z-c | < R$

   * A infinite open sector defined by $\alpha < arg(z) < \beta$

   * A finite open sector defined by $\alpha < arg(z) < \beta$ and $r<|z|<R$

(6) Describe the following regions:

   * $|3z-6|\le 9$

   * $\left|\frac{z-1}{1+z}\right|<1$

   * $|z| \le |1-z|$

# The Complex Exponential

In previous courses you would have seen that functions of $e^x$, $\mathrm{sin} \, x $ and $\mathrm{cos} \,x $ can be expanded in power series:

$$e^x = 1 + \frac{x}{1!} + \frac{x^2}{2!}+ \frac{x^3}{3!} + \frac{x^4}{4!}+ \mathcal{O}\left(x^5\right) \tag{12}$$

$$\mathrm{sin} \, x = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \frac{x^9}{9!} + \mathcal{O}\left(x^{11}\right)\tag{13}$$

$$\mathrm{cos} \,x = 1 - \frac{x^2}{2!} + \frac{x^4}{4!}- \frac{x^6}{6!}+ \frac{x^8}{8!} + \mathcal{O}\left(x^{10}\right) \tag{14}$$

where, e.g., $\mathcal{O}(x^5)$ means "terms of order $x^5$ and greater", and $!$ is the factorial operator. If we now set $x=i\theta$ and substitute this into $e^x$ we get:

$$
\begin{eqnarray}
e^{i\theta} &=& 1 + \frac{{i\theta}}{1!} + \frac{({i\theta})^2}{2!}+ \frac{({i\theta})^3}{3!}+ \frac{({i\theta})^4}{4!}+ \frac{({i\theta})^5}{5!}+ \mathcal{O}\left(x^{6}\right) \\
\, &=&1 + \frac{{i\theta}}{1!} - \frac{{\theta}^2}{2!}- \frac{i{\theta}^3}{3!}+ \frac{{\theta}^4}{4!}+ \frac{i{\theta}^5}{5!} + \mathcal{O}\left(x^{6}\right)\\
\, &=&\left(  1 - \frac{{\theta}^2}{2!} + \frac{{\theta}^4}{4!} +\mathcal{O}\left(x^{6}\right)\right) + i\left(\frac{{\theta}}{1!} - \frac{{\theta}^3}{3!} + \frac{{\theta}^5}{5!} + \mathcal{O}\left(x^{7}\right) \right) \tag{15}\\
\, &=& \mathrm{cos}\,\theta + i \mathrm{sin}\,\theta. \\
\end{eqnarray}
$$

This is **Euler's formula** that defines $e^{i\theta}$ for any given real-valued $\theta$ as a vector of unit length making angle $\theta$ with the real axis.  Alternatively, the Euler Formulae for the cos $\theta$ and sin $\theta$ functions in terms of complex exponentials are given by:

$$
{\rm cos}\,\theta = \frac{{\rm e}^{i\theta} + {\rm e}^{-i\theta}}{2} \quad \mathrm{and} \quad {\rm sin}\,\theta = \frac{{\rm e}^{i\theta} - {\rm e}^{-i\theta}}{2i}. \tag{16}
$$

Note that multiplying $e^{i\theta}$ by $r$ gives $re^{i\theta}$ which is a complex number having modulus $r$ and argument $\theta$.   Thus, $z=re^{i\theta}$ is convenient when we wish to display the **polar coordinates** of $z$ and the notation $z=x+iy$ is convenient when we want to display the **Cartesian components** of $z$.  For example, all of the following are equivalent representations:

$$1+i = \sqrt{2}e^{i\pi/4} = -\sqrt{2} e^{i5\pi/4} = -\sqrt{2} e^{-i3\pi/4}. \tag{17}$$

# Geometric aspects of Multiplication and Division

Based on the rules from the study of exponents involving real numbers, it is fairly straightforward to suppose that the following operations involving imaginary exponents will hold: $e^{ia}e^{ib}$,  $e^{ia}/e^{ib}$ and  $(e^{ia})^b$ where $[a,b]\in\Re$.

To illustrate this, let's define $z=re^{i\theta}$ and $w=se^{i\gamma}$ where $[r,s,\theta,\gamma]\in\Re$. We can then specify the following multiplication and division operations:

**Multiplication**:  We obtain the following:

$$  z\,w = r\, e^{i\theta} s\,e^{i\gamma} = r\,s\, e^{i\theta}\,e^{i\gamma} = r\,s\, e^{i(\theta+\gamma)}. \tag{18}$$

The last argument tells us that the modulus is $|z\,w|=|z|\,|w|$ and the argument is $arg(zw) = arg(z)+arg(w)$.

**Division:** We obtain the following:

$$\frac{z}{w} = \frac{r\, e^{i\theta}}{s\,e^{i\gamma}} = \frac{r}{s}\, e^{i(\theta-\gamma)}, \tag{19}$$

where the modulus is given by $|z\,/\,w|= |z|\,/\,|w|$ and the argument by $arg(z\,/\,w)=arg(z)-arg(w)$.

### QUICK PROBLEMS:

(7) Prove **de Moivre's Theorem**:

$$ \left(\cos\,\theta+i\sin\,\theta\right)^{n}=\cos\,n\theta+i\sin\,n\theta, \tag{20}$$

where $n$ is an integer.  (Hint: it's straightforward with Euler's Identity. You can also do this more formally through an inductive proof that does not rely on Euler's Identity.)

(8) By squaring $e^{i\theta}=\mathrm{cos}\, \theta + i\,\mathrm{sin} \theta$, find formulas for $\mathrm{cos}\,2\theta$ and $\mathrm{sin}\,2\theta$ in terms of $\mathrm{cos}\,\theta$ and $\mathrm{sin}\,\theta$.

(9) Find formulas for $\mathrm{cos}\,3\theta$ and $\mathrm{sin}\,3\theta$ in terms of $\mathrm{cos}\,\theta$ and $\mathrm{sin}\,\theta$.

# Roots of a Complex Number

Suppose we are asked for the cube root of $8i$.  We would seek seek $z$ such that $z^3=8i$. First, we recognize that the **principle value** of $i=e^{i\pi/2}$. Recognizing that we can have $2\pi n$ multiples of this (where $n$ is integer) allows us to write:

$$8i = 8 e^{i(\pi/2+2\pi n)},  \quad n=0,\pm 1,\pm 2, ... \tag{21}$$

Thus, we can enforce the cube root on both sides:

$$\begin{eqnarray}
z^3 &=& 8 e^{i(\pi/2+2\pi n)}\tag{22a} \\
z   &=& \left(8 e^{i(\pi/2+2\pi n)} \right)^\frac{1}{3} \tag{22b}\\
z   &=& 8^\frac{1}{3} \left(e^{i(\pi/2+2\pi n)} \right)^\frac{1}{3} \tag{22c}\\
z   &=& 2 e^{i(\pi/6+2\pi n/3)}, \quad n=0,\pm 1,\pm 2, ... \tag{22d}\\
\end{eqnarray}$$

It looks like there are an infinite number of different values!  However, let's look at this in complex plane using an [Argand Diagram](http://mathworld.wolfram.com/ArgandDiagram.html).  In the example below, I have defined an *argand(a)* function that will plot the location of a given set of points $z$ on the complex plane and then draw the vector from the origin. In this example, I am plotting the points in the above example from $n=0$ to $n=\pm 50$.

```python
def argand(a):
    ## . . Loop over and plot all points given in the passed array a
    for x in range(len(a)):
        plt.plot([0,a[x].real],[0,a[x].imag],'ro-',label='python')
    limit=np.max(np.ceil(np.absolute(a))) # set limits for axis
    plt.ylabel('Imaginary')
    plt.xlabel('Real')
    plt.axis('square')

    # . . Draw horizontal and vertical axes
    plt.plot([0,0],[-6,6],'k',linewidth=0.5)
    plt.plot([-6,6],[0,0],'k',linewidth=0.5)

## . . Create array with different values of n from -50<n<50
nn = np.arange(-50,50,1)
a = 2*(np.cos(np.pi*(1/6+2/3*nn)) + 1j*np.sin(np.pi*(1/6+2/3*nn)))
argand(a)
plt.xlim((-2.5,2.5))
plt.ylim((-2.5,2.5))
plt.show()
```

**Figure 7. Plotting the three cube roots of $z=8i$ using the defined Argand function.**

In fact, as illustrated in Figure 7, there are only three distinct answers!  This should not be surprising since $z^3=8i$ is a [cubic equation](https://en.wikipedia.org/wiki/Cubic_function) and therefore we should expect exactly three solutions.

### QUICK PROBLEMS

(10) Solve the following equations and show the answers on an Argand diagram:

   * $z^2+4=0$

   * $z^3-8=0$

   * $z^4=81(1+i)/\sqrt{2}$

   * $z^3+i=0$

   * $z^6=-64$

# Complex Functions

Complex analysis is concerned with complex functions that are differentiatable in some defined domain.  However, we should define what we mean by **complex function**.  A function $f$ on a complex domain $S$ assigns to every $z \in S$ a complex number $w$ called the value of $f$ at $z$:

$$ w = f(z). \tag{23}$$

For example, $w = f(z) = z^2+3z$ is a complex function defined for all of $z$, such that the domain $S$ is the whole complex plane.  Generally speaking, $w$ is itself a complex number, $w = u+iv$, where both $u$ and $v$ can be dependent on $z=x+iy$.  Thus, we write:

$$w = f(z) = u(x,y) + iv(x,y) \tag{24}$$

where $u(x,y)$ and $v(x,y)$ are functions of two real variables.

**EXAMPLE**

**Q1:** If $w=f(z)=z^2$, what are the associated $u(x,y)$ and $v(x,y)$?

**A1:** Because $z=x+iy$ we may write the following:

$$w = z^2 = (x+iy)^2 = x^2+2ixy+(iy)^2 = x^2-y^2 + 2ixy$$

Therefore, by inspection we see that $u(x,y) = x^2-y^2$ and $v(x,y) = 2xy$.

## Limits and Continuity

A function $f(z)$ is said to have the **limit** $l$ as $z$ approaches some $z_0$:

$$ \lim_{z\rightarrow z_0} f(z) = l, \tag{25}$$

if and only if $f$ is defined in a neighbourhood of $z_0$ and if the values of $f$ are "close" to $l$ for all $z$ "close" to $z_0$. More precisely, in the disk $|z-z_0|<\delta$ we have $|f(z)-l|<\epsilon$.  This is similar to how limits are defined in the calculus of 1D real variables; however, one big difference is that $z$ may approach $z_0$ **from any direction** (i.e., in $\theta$, not just from the left (-) and right (+) along the real axis).

A function $f(z)$ is said to be **continuous** at $z=z_0$ if and only if $f(z_0)$ is defined and

$$\lim_{z\rightarrow z_0} f(z) = f(z_0). \tag{26}$$

$f(z)$ is said to be **continuous in a domain** if it is continuous at each point in the domain.

**EXAMPLE**

**Q1:** Is $f(z)=\sum_{n=1}^{10} z^n$ continuous at $z=0$?

**A1:** Let's look at the following limit:

$$\lim_{z\rightarrow 0}\sum_{n=1}^{10} z^n  = \lim_{r\rightarrow 0}\sum_{n=1}^{10} (re^{i\theta })^n = \lim_{r\rightarrow 0}\sum_{n=1}^{10} r^ne^{in\theta }=0.$$

Thus, because the limit approaches $l=0$, the function is continuous at $z=0$.

### QUICK PROBLEMS

(11) Find out (and give a reason) whether $f(z)$ is continuous at $z=0$ if $f(0)=0$ and for $z\neq0$ the function is equal to:
   * $\frac{\Im{(z)}}{(1-|z|)}$

   * $\frac{\Im{(z^2)}}{|z|^2}$

(Hint: Use the polar notation, assert $r\rightarrow 0$ and see if this is valid for all $\theta$).

# Derivative of complex function

**Provided that the limit exists**, the **derivative** of a complex function $f$ at a point $z_0$ is written $f^\prime(z_0)$ and is defined by

$$f^\prime(z_0) = \lim_{\Delta z\rightarrow 0} \frac{f(z_0+\Delta z)-f(z_0)}{\Delta z}. \tag{28} $$

Then $f(z)$ is said to be **differentiable** at $z_0$.  If we write $\Delta z = z-z_0$ then we can change the limit argument and denominator as follows:

$$f^\prime(z_0) = \lim_{z\rightarrow z_0} \frac{f(z_0+\Delta z)-f(z_0)}{z-z_0}. \tag{29}$$

Note that differentiability of $f$ at $z_0$ requires along whatever path $z$ approaches $z_0$ the expressions always approaches a certain value and that all these values are equal.

**EXAMPLES**

**Q1**: Is the function $f(z)=z^2$ differentiable for all $z$? If yes, what is its derivative?.

**A1:** Yes, the function $f(z)=z^2$ is differentiable for all $z$ and has the derivative $f^\prime(z)=2z$ because:

$$ f^\prime(z) = \lim_{\Delta z\rightarrow 0} \frac{ (z+\Delta z)^2 -z^2 }{\Delta z} = \lim_{\Delta z\rightarrow 0} \frac{ z^2+2z\Delta z +\Delta z^2 -z^2 }{\Delta z}= \lim_{\Delta z\rightarrow 0} \left(2z +\Delta z\right) = 2z. \tag{30}$$

**Q2**: Is the function $f(z)=\overline{z}$ differentiable for all $z$? If yes, what is its derivative?.

**A2:** No, the function $f(z)=\overline{z}$ is not differentiable because:

$$ \frac{f(z+\Delta z)-f(z)}{\Delta z} = \frac{\overline{z+\Delta z}-\overline{z}}{\Delta z} =\frac{\overline{\Delta z}}{\Delta z} = \frac{\Delta x - i\Delta y}{\Delta x + i \Delta y}. \tag{31}$$

If $\Delta y=0$ and we approach $z_0$ along $\Delta x$ then this expression is equal to $1$; however, if $\Delta x=0$ and we approach $z_0$ along $\Delta y$ then the expression is $-1$.  Thus, the limit does not exist and $\overline{z}$ is not differentiable.

### QUICK PROBLEMS

(12) Are the following functions differentiable? If so, what is the derivative?:

   * $f(z) = \left(z^3 + i\right)^2 $

   * $f(z) =\frac{i}{(1-z)^2}$

## Analytic Functions

Complex analysis is concerned with the theory and application of [analytic functions](https://en.wikipedia.org/wiki/Analytic_function) - i.e., functions that are differentiable in some domain such that we can "do calculus".  An **analytic function** is said to be analytic in a domain $D$ if $f(z)$ is defined and differentiable at all points in $D$. The function $f(z)$ is said to be **analytic** at a point $z=z_0$ in D if f(z) is analytic in a neighbourhood of $z_0$.

Importantly, non-negative integer powers $\{1, z, z^2, z^3, ...\}$ are analytic in the complex plane and so are [polynomials](https://en.wikipedia.org/wiki/Polynomial) of the form:

$$f(z) = c_0 + c_1z + c_2 z^2 + c_3z^3 + ... + c_nz^n = \sum_{n=0}^{N} c_n z^n. \tag{32} $$

The quotient of two polynomials $g(z)$ and $h(z)$:

$$f(z) = \frac{g(z)}{h(z)} \tag{33}$$

is called a [rational function](https://en.wikipedia.org/wiki/Rational_function).  Rational function $f(z)$ is analytic except at points where $h(z)=0$ [also called the poles of $f(z)$].  We also can assume that many common factors of of $g$ and $h$ cancel out.  Rational functions will be quite important in the context of designing filters in the later sections of this course (Modules 12 and 13).

# Cauchy-Riemann Equations

Finally, we are going to conclude this overview of complex numbers with a brief discussion of the [Cauchy-Riemann equations](https://en.wikipedia.org/wiki/Cauchy%E2%80%93Riemann_equations), which are arguably the most important equations in [complex analysis](https://en.wikipedia.org/wiki/Complex_analysis).  These equations provide a test criterion for the analycity of a complex function $w=f(z)= u(x,y)+iv(x,y)$.  Roughly, $f$ is analytic in a domain $D$ if and only if the first partial derivatives of $u$ and $v$ satisfy the two Cauchy-Riemann equations:

$$\begin{eqnarray}
\frac{\partial u}{\partial x} &=& \frac{\partial v}{\partial y}, \tag{34a}\\
\frac{\partial u}{\partial y} &=&-\frac{\partial v}{\partial x}, \tag{34b}\\
\end{eqnarray}$$

everywhere in domain $D$.

**Q1**: Show that $f(z)=z^2$ is analytic for all $z$.

**A1:** First, we observe that $z^2 =x^2-y^2+i2xy$ and thus $u = x^2-y^2$ and $v=2xy$.  We then apply the Cauchy-Riemann test:

$$\begin{eqnarray}
\frac{\partial u}{\partial x} = 2x, &\quad \, \quad & \frac{\partial u}{\partial y}=-2y, \tag{35a}\\
\frac{\partial v}{\partial x} = 2y, &\quad\,\quad & \frac{\partial v}{\partial y}=2x.
\tag{35b}\\
\end{eqnarray}$$

Thus, $f(z)=z^2$ satifies the Cauchy-Riemann equations and is analytic for all of $z$.

### QUICK PROBLEMS

(13) Determine whether the following functions are analytic:

   * $f(z) = z^4$

   * $f(z) = e^{2x}\left(\mathrm{cos}\, y + i  \,\mathrm{sin} \, y\right)$


text:
# Neural Nets multiple
Version 0.1, in `nnn`

# Note:
This is the source for `nnn.py`, the foundation library of these neural net experiments. It contains a simple library, and a quantity, never sufficient, of test code, guarded by `if __name__ = '__main__'`. It makes various noises as it runs, but should not fail any asserts. See the bottom of the file for the procedure to produce the importable library file `nnn.py`.

Should do [Working efficiently with jupyter lab](https://florianwilhelm.info/2018/11/working_efficiently_with_jupyter_lab/)

```python
from collections.abc import MutableSequence
from itertools import cycle
import numpy as np
from types import FunctionType
```

A network built of components which:
1. accept an ordered set of reals (we'll use `numpy.array`, and  call them vectors) at the input port and produce another at the output port - this is forward propagation. ${\displaystyle f\colon \mathbf {R} ^{n}\to \mathbf {R} ^{m}}$
1. accept an ordered set of reals at the output port, representing the gradient of the loss function at the output, and produce the gradient of the loss function at the input port - this is back propagation, aka backprop. ${\displaystyle b\colon \mathbf {R} ^{m}\to \mathbf {R} ^{n}}$
1. from the gradient of the loss function at the output, calculate the partial of the loss function w.r.t the internal parameters ${\displaystyle \frac{\partial E}{\partial w} }$
1. accept a scalar $\eta$ to control the adjustment of internal parameters. _Or is this effected by scaling the loss gradient before passing? YES_
1. update internal parameters ${\displaystyle w \leftarrow w - \eta \frac{\partial E}{\partial w} }$

```python
class Layer(MutableSequence):
    def __init__(self):
        pass

    def __call__(self, x):
        """Compute response to input"""
        raise NotImplementedError

    def backprop(self, output_delE):
        """Use output error gradient to adjust internal parameters, return gradient of error at input"""
        raise NotImplementedError

    def state_vector(self):
        """Provide the layer's learnable state as a vector"""
        raise NotImplementedError

    def set_state_from_vector(self, sv):
        """Set the layer's learnable state from a vector"""
        raise NotImplementedError

    # __delitem__, __getitem__, __len__, __setitem__, insert

    def __delitem__(self, ix):
        raise NotImplementedError

    def __getitem__(self, ix):
        raise NotImplementedError

    def __len__(self):
        raise NotImplementedError

    def __setitem__(self, ix, value):
        raise NotImplementedError

    def insert(self, ix, value):
        raise NotImplementedError
```

A network built of a cascade of layers:

```python
class Network:
    version_string = '0.6'

    def __init__(self):
        self.layers = []
        self.eta = 0.1 #FIXME

    def extend(self, net):
        self.layers.append(net)
        return self

    def __call__(self, input):
        v = input
        for net in self.layers:
            v = net(v)
        return v

    def learn(self, facts, eta=None):
        eta = eta or self.eta
        for x, ideal in facts:
            y = self(x)
            e = y - ideal
            egrad = e * eta / np.atleast_2d(e).shape[-2] #TESTME: normalization
            for net in reversed(self.layers):
                egrad = net.backprop(egrad)
        #loss = float(e.dot(e.T))/2.0
        loss = np.einsum('...ij,...ij', e, e) / (2.0 * np.atleast_2d(e).shape[-2])
        self.eta = eta
        return loss

    def losses(self, facts):
        return [np.einsum('...ij,...ij', e, e) / (2.0 * np.atleast_2d(e).shape[-2]) \
                for e in (self(x) - ideal for x, ideal in facts)]

    def state_vector(self):
        """Provide the network's learnable state as a vector"""
        return np.concatenate([layer.state_vector() for layer in self.layers])

    def set_state_from_vector(self, sv):
        """Set the layer's learnable state from a vector"""
        i = 0
        for layer in self.layers:
            lsvlen = len(layer.state_vector())
            layer.set_state_from_vector(sv[i:i+lsvlen])
            i += lsvlen
        return self
```

___

## Useful Layers

### Identify

```python
class IdentityLayer(Layer):
    def __call__(self, x):
        return x

    def backprop(self, output_delE):
        return output_delE

    def state_vector(self):
        return np.array([])

    def set_state_from_vector(self, sv):
        pass
```

### Affine
A layer that does an [affine transformation](https://mathworld.wolfram.com/AffineTransformation.html) aka affinity, which is the classic fully-connected layer with output offsets.

$$ \mathbf{M} \mathbf{x} + \mathbf{b} = \mathbf{y} $$
where
$$
\mathbf{x} = \sum_{j=1}^{n} x_j \mathbf{\hat{x}}_j \\
\mathbf{b} = \sum_{i=1}^{m} b_i \mathbf{\hat{y}}_i \\
\mathbf{y} = \sum_{i=1}^{m} y_i \mathbf{\hat{y}}_i
$$
and $\mathbf{M}$ can be written
$$
\begin{bmatrix}
    m_{1,1} & \dots & m_{1,n} \\
    \vdots & \ddots & \vdots \\
    m_{m,1} & \dots & m_{m,n}
\end{bmatrix} \\
$$

#### Error gradient back-propagation
$$
\begin{align}
 \frac{\partial loss}{\partial\mathbf{x}}
  &= \frac{\partial loss}{\partial\mathbf{y}} \frac{\partial\mathbf{y}}{\partial\mathbf{x}} \\
  &= \mathbf{M}^\mathsf{T}\frac{\partial loss}{\partial\mathbf{y}}
\end{align}
$$

#### Parameter adjustment
$$
 \frac{\partial loss}{\partial\mathbf{M}}
 = \frac{\partial loss}{\partial\mathbf{y}} \frac{\partial\mathbf{y}}{\partial\mathbf{M}}
 = \frac{\partial loss}{\partial\mathbf{y}} \mathbf{x} \\
 \frac{\partial loss}{\partial\mathbf{b}}
 = \frac{\partial loss}{\partial\mathbf{y}} \frac{\partial\mathbf{y}}{\partial\mathbf{b}}
 = \frac{\partial loss}{\partial\mathbf{y}}
$$

#### Adapting to `numpy`

In `numpy` it is more convenient to use row vectors, particularly for calculating the transform on multiple inputs in one operation. We use the identity $ \mathbf{M} \mathbf{x} = (\mathbf{x} \mathbf{M}^\mathsf{T})^\mathsf{T}.$ To avoid cluttering names, we will use `M` in the code below to hold $\mathbf{M}^\mathsf{T}$.

```python
class AffineLayer(Layer):
    """An affine transformation, which is the classic fully-connected layer with offsets.

    The layer has n inputs and m outputs, which numbers must be supplied
    upon creation. The inputs and outputs are marshalled in numpy arrays, 1-D
    in the case of a single calculation, and 2-D when calculating the outputs
    of multiple inputs in one call.
    If called with 1-D array having shape == (n,), e.g numpy.arange(n), it will
    return a 1-D numpy array of shape (m,).
    If called with a 2-D numpy array, input shall have shape (k,n) and will return
    a 2-D numpy array of shape (k,m), suitable as input to a subsequent layer
    that has input width m.
    """
    def __init__(self, n, m, multiples=1):
        self.M = np.empty((multiples, n, m))
        self.b = np.empty((multiples, 1, m))
        self.randomize()

    def randomize(self):
        self.M[:] = np.random.randn(*self.M.shape)
        self.b[:] = np.random.randn(*self.b.shape)

    def __call__(self, x):#input[a,b,c] is net a, batch member b, dimension c
        self.input = x
        self.output = x @ self.M + self.b
        return self._squeezeit(self.output, ref=self.input)

    def backprop(self, output_delE):
        input_delE = output_delE @ self.M.transpose(0,2,1)
        self.M -= np.einsum('...ki,...kj->...ji', np.atleast_2d(output_delE), np.atleast_2d(self.input))
        self.b -= np.sum(np.atleast_2d(output_delE), 0)
        #return input_delE.squeeze(tuple(range(input_delE.ndim - output_delE.ndim)))
        return self._squeezeit(input_delE, ref=output_delE)

    def state_vector(self):
        return np.concatenate((self.M.ravel(), self.b.ravel()))

    def set_state_from_vector(self, sv):
        """Set the layer's learnable state from a vector"""
        l_M = len(self.M.ravel())
        l_b = len(self.b.ravel())
        self.M[:] = sv[:l_M].reshape(self.M.shape)
        self.b[:] = sv[l_M : l_M + l_b].reshape(self.b.shape)

    def _squeezeit(self, v, ref):
        if self.M.shape[0] == 1:
            return v.squeeze(tuple(range(v.ndim - ref.ndim)))
        return v

    def __getitem__(self, ix):
        rv = self.__class__(*self.M.shape[-2:])
        rv.M = self.M[ix].reshape((1,) + self.M[ix].shape)
        rv.b = self.b[ix].reshape((1,) + self.b[ix].shape)
        return rv

    def __setitem__(self, ix, a):
        if not np.shares_memory(self.M[ix], a.M):
            self.M[ix,:] = a.M
        if not np.shares_memory(self.b[ix], a.b):
            self.b[ix,:] = a.b

    def __len__(self):
        return self.M.shape[0]
```

### Map
Maps a scalar function on the inputs, for e.g. activation layers.

```python
class was_MapLayer(Layer):
    """Map a scalar function on the input taken element-wise"""
    def __init__(self, fun, dfundx):
        self.vfun = np.vectorize(fun)
        self.vdfundx = np.vectorize(dfundx)

    def __call__(self, x):
        self.input = x
        return self.vfun(x)

    def backprop(self, output_delE):
        input_delE = self.vdfundx(self.input) * output_delE
        return input_delE

    def state_vector(self):
        return np.array([])

    def set_state_from_vector(self, sv):
        pass
```

```python
class MapLayer(Layer):
    """Map a scalar function on the input taken element-wise"""
    def __init__(self, f, df):
        if hasattr(f, '__call__') and hasattr(df, '__call__'):
            self.f, self.df = np.vectorize(f), np.vectorize(df)
            self.single = True
        elif hasattr(f, '__getitem__') and hasattr(f[0], '__call__') and \
                hasattr(df, '__getitem__') and hasattr(df[0], '__call__'):
            self.f = [np.vectorize(fun) for fun in f]
            self.df = [np.vectorize(fun) for fun in df]
            self.single = False
        else:
            raise ValueError("MapLayer requires a function or sequence of functions")

    def __call__(self, x):
        self.input = x
        if self.single:
            return self.f(x)
        return np.array([f(v) for f,v in zip(cycle(self.f),x)])

    def backprop(self, output_delE):
        if self.single:
            return self.df(self.input) * output_delE
        return np.array([df(v) for df,v in zip(cycle(self.df), output_delE)])

    def state_vector(self):
        return np.array([])

    def set_state_from_vector(self, sv):
        pass
```

---

# Tests
*Incomplete* \
Also `unittest` the `.py` version with a separate test script, see `test-nn_v3.py`.

Make a few test arrays:

```python
if __name__ == '__main__':
    one_wide = np.atleast_2d(np.arange(1*4)).reshape(-1,1)
    print(f"one_wide is:\n{one_wide}")
    two_wide = np.arange(2*4).reshape(-1,2)
    print(f"two_wide is:\n{two_wide}")
    three_wide = np.arange(3*4).reshape(-1,3)
    print(f"three_wide is:\n{three_wide}\n")
```

    one_wide is:
    [[0]
     [1]
     [2]
     [3]]
    two_wide is:
    [[0 1]
     [2 3]
     [4 5]
     [6 7]]
    three_wide is:
    [[ 0  1  2]
     [ 3  4  5]
     [ 6  7  8]
     [ 9 10 11]]

## Tooling for Testing

```python
if __name__ == '__main__':
    import sympy

    class VC():
        def grad(f, x, eps=1e-6):
            epsihat = np.eye(x.shape[-1]) * eps
            yp = np.apply_along_axis(f, 1, x + epsihat)
            ym = np.apply_along_axis(f, 1, x - epsihat)
            return (yp - ym)/(2 * eps)

        def tensor_grad(f, x, eps=1e-6):
            return np.apply_along_axis(lambda v: VC.grad(f, v, eps), 1, x)

    def closenuf(a, b, tol=0.001):
        return np.allclose(a, b, rtol=tol)

    def arangep(n, starting_index=0):
        sympy.sieve.extend_to_no(starting_index + n)
        return np.array(sympy.sieve._list[starting_index:starting_index + n])
```

```python
VC.grad(lambda x:x**2, three_wide[1])
```

    array([[ 6.,  0.,  0.],
           [ 0.,  8.,  0.],
           [ 0.,  0., 10.]])

```python
VC.grad(lambda x: x @ np.array([[0,1], [-1, 0]]), two_wide[:2])
```

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

```python
VC.tensor_grad(lambda x:x**2, three_wide)
```

    array([[[ 0.        ,  0.        ,  0.        ],
            [ 0.        ,  2.        ,  0.        ],
            [ 0.        ,  0.        ,  4.        ]],

           [[ 6.        ,  0.        ,  0.        ],
            [ 0.        ,  8.        ,  0.        ],
            [ 0.        ,  0.        , 10.        ]],

           [[12.        ,  0.        ,  0.        ],
            [ 0.        , 14.        ,  0.        ],
            [ 0.        ,  0.        , 15.99999999]],

           [[17.99999999,  0.        ,  0.        ],
            [ 0.        , 19.99999999,  0.        ],
            [ 0.        ,  0.        , 21.99999998]]])

---

Input to a layer can be a single (row) vector, or a vertical stack of row vectors,
a 2-d array that resembles a matrix. We need to test each layer class with both single and stacked input.

## Identity layer

```python
if __name__ == '__main__':
    iL = IdentityLayer()

    # It's transparent from input to output
    assert np.equal(iL(np.arange(5)), np.arange(5)).all()
    assert (iL(three_wide) == three_wide).all()

    # It back-propagates the loss gradient without alteration
    assert np.equal(iL.backprop(np.arange(7)), np.arange(7)).all()
    assert (iL.backprop(three_wide) == three_wide).all()

    # It works for stacked input
    # (see above)
```

## Map layer

### Test single vector input behavior

```python
if __name__ == '__main__':
    mL = MapLayer(lambda x:x**2, lambda d:2*d)

    # It applies the forward transformation
    assert np.equal(mL(np.array([-2,1,3])), np.array([4,1,9])).all()

    # It back-propagages the loss gradient
    x = np.array([1,2,2])
    y = mL(x)

    # for loss function, use L2-distance from some ideal
    # (divided by 2, for convenient gradient = error)
    ideal = np.array([2,3,5])
    loss = lambda v: (v - ideal).dot(v - ideal) / 2.0
    loss_at_y = loss(y)
    print(f"x = {x}, y = {y}, loss at y = {loss_at_y}")

    # find numerical gradient of loss function at y, the layer output
    grad_y = VC.grad(loss, y)
    print(f"∇𝑙𝑜𝑠𝑠(𝑦) = {grad_y}")

    # find the numerical gradient of the loss w.r.t. the input of the layer
    grad_x = VC.grad(lambda x:loss(mL(x)), x)
    print(f"∇𝑙𝑜𝑠𝑠(𝑥) = {grad_x}")

    # The backprop method does the same
    _ = mL(x) # Make sure the last x is in the right place
    in_delE = mL.backprop(grad_y)
    print(f"backprop({grad_y}) = {in_delE}")
    assert closenuf(in_delE, grad_x)

    # The backprop operation did not change the behavior
    assert np.equal(mL(x), y).all()
```

    x = [1 2 2], y = [1 4 4], loss at y = 1.5
    ∇𝑙𝑜𝑠𝑠(𝑦) = [-1.  1. -1.]
    ∇𝑙𝑜𝑠𝑠(𝑥) = [-2.  4. -4.]
    backprop([-1.  1. -1.]) = [-2.  4. -4.]

### Test stacked-vectors input:

```python
if __name__ == '__main__':
    mL = MapLayer(lambda x:x**2, lambda d:2*d)

    two_wide_sq = np.array([[ 0,  1],
                            [ 4,  9],
                            [16, 25],
                            [36, 49]])
    # It applies the forward transformation
    assert np.equal(mL(two_wide), two_wide_sq).all()

    # It back-propagages the loss gradient
    x = two_wide
    y = mL(x)

    # for loss function, use L2-distance from some ideal
    # (divided by 2, for convenient gradient = error)
    ideal = two_wide * 2 + 11
    #print(y - ideal)
    #loss = lambda v: (v - ideal).dot(v - ideal) / 2.0
    loss = lambda v: np.einsum('ij,ij', v-ideal, v-ideal) / (2 * v.shape[0])
    loss_at_y = loss(y)
    print(f"x =\n{x}\ny =\n{y}, loss = {loss_at_y}\n")

    # find numerical gradient of loss function at y, the layer output
    grad_y = VC.tensor_grad(loss, y)
    print(f"∇𝑙𝑜𝑠𝑠(𝑦) = {grad_y}\n")

    # find the numerical gradient of the loss w.r.t. the input of the layer
    grad_x = VC.tensor_grad(lambda x:loss(mL(x)), x)
    print(f"∇𝑙𝑜𝑠𝑠(𝑥) = {grad_x}\n")

    # The backprop method does the same
    _ = mL(x) # Make sure the last x is in the right place
    in_delE = mL.backprop(grad_y)
    print(f"backprop({grad_y}) =\n{in_delE}")
    assert closenuf(in_delE, grad_x)

    # The backprop operation did not change the behavior
    assert np.equal(mL(x), y).all()
```

    x =
    [[0 1]
     [2 3]
     [4 5]
     [6 7]]
    y =
    [[ 0  1]
     [ 4  9]
     [16 25]
     [36 49]], loss = 152.5

    ∇𝑙𝑜𝑠𝑠(𝑦) = [[-34.00000003 -36.00000002]
     [-26.00000002 -19.99999998]
     [ -1.99999999  12.00000001]
     [ 37.9999999   59.99999985]]

    ∇𝑙𝑜𝑠𝑠(𝑥) = [[   0.          -71.99999999]
     [-103.99999996 -120.00000001]
     [ -16.          120.00000002]
     [ 456.0000001   840.00000027]]

    backprop([[-34.00000003 -36.00000002]
     [-26.00000002 -19.99999998]
     [ -1.99999999  12.00000001]
     [ 37.9999999   59.99999985]]) =
    [[  -0.          -72.00000005]
     [-104.00000008 -119.99999987]
     [ -15.99999996  120.00000012]
     [ 455.99999885  839.99999788]]

### Test stacked-vectors input with multiple functions:

```python
if __name__ == '__main__':
    # A list of lambdas is applied cyclically to the argument sequence
    mL = MapLayer([lambda x,i=i:x**i for i in range(2, 4)],
                  [lambda x,i=i:i*x**(i-1) for i in range(2, 4)])

    x = np.ones(5) * 10
    y = np.array([100, 1000, 100, 1000, 100])
    dy = np.array([20, 300, 20, 300, 20])
    assert np.equal(mL(x), y).all()
    assert np.equal(mL.backprop(x), dy).all()

    # When applied to a multi, it cycles through application correctly
    x = np.arange(2*3*4).reshape(2,3,4)
    y = np.array(
        [[[    0,     1,     4,     9],
          [   16,    25,    36,    49],
          [   64,    81,   100,   121]],

         [[ 1728,  2197,  2744,  3375],
          [ 4096,  4913,  5832,  6859],
          [ 8000,  9261, 10648, 12167]]])
    assert np.equal(mL(x), y).all()
```

## Affine layer

### Test single vector input behavior

Test, for single input-vector operations:
* input and output widths
* state vector setting and getting
* forward calculation

```python
if __name__ == '__main__':
    # Affine
    a = AffineLayer(2,3)

    # The input and output widths are correct
    assert a(np.arange(2)).shape == (3,)

    # Its internal state can be set
    a.set_state_from_vector(np.arange(9))
    # and read back
    assert (a.state_vector() == np.arange(9)).all()
    # NOTE: The two assertions below are commented out because they depend
    # on white-box knowledge, and are duplicative of other tests
    #assert np.equal(a.M, np.array([[0, 1, 2],
    #                               [3, 4, 5]])).all()
    #assert np.equal(a.b, np.array([6, 7, 8])).all()

    # Its internal state observed using numerical gradient is correct
    x = np.random.rand(2)
    y = a(x)
    dydx = VC.grad(a, x)
    b = y - x.dot(dydx)
    #print(dydx, b)
    #print(dydx, np.arange(6).reshape(2,-1))
    assert closenuf(dydx, np.arange(6).reshape(2, -1))
    #print(b, np.arange(6, 9))
    assert closenuf(b, np.arange(6, 9))

    # It performs a single-input forward calculation correctly
    x = np.array([2, 1])
    y = a(x)
    #print(f"a.M is:\n{a.M}\na.b is {a.b}\nx is: {x}\ny is: {y}\n")
    assert (y == np.array([9, 13, 17])).all()

    # It performs a different single-input forward calculation correctly
    a.set_state_from_vector(np.array([ 2,  3,  5,  7, 11, 13, 17, 19, 23]))
    x = np.array([[29, 31]])
    y = a(x)
    assert (y == np.array([[292, 447, 571]])).all()
```

Test, for single input-vector operations:
* back-propagation of the loss gradient
* learning (change in forward function) from the back-prop operation

```python
if __name__ == '__main__':
    # Affine
    a = AffineLayer(2,3)
    a.set_state_from_vector(np.arange(9))

    # Doing a single-input-vector calculation
    x = np.array([2, 1])
    y = a(x)
    assert np.equal(y, np.array([9, 13, 17])).all()

    # It back-propagages the loss gradient
    ideal = np.array([11,12,10])
    loss = lambda v: (v - ideal).dot(v - ideal) / 2.0
    loss_at_y = loss(y)
    print(f"x = {x}, y = {y}, loss = {loss_at_y}")
    grad_y = VC.grad(loss, y)
    print(f"∇𝑙𝑜𝑠𝑠(𝑦) = {grad_y}")
    grad_x = VC.grad(lambda x:loss(a(x)), x)
    print(f"∇𝑙𝑜𝑠𝑠(𝑥) = {grad_x}")

    # Back-propagate the loss gradient from layer output to input
    _ = a(x) # Make sure the last x is in the right place
    out_delE = grad_y * 0.1 # Backprop one-tenth of the loss gradient
    in_delE = a.backprop(out_delE)
    print(f"backprop({out_delE}) = {in_delE}")

    # The loss gradient back-propagated to the layer input is correct
    assert closenuf(in_delE / 0.1, grad_x)

    # And how did the learning affect the layer?
    print(f"Now a({x}) = {a(x)}, loss = {loss(a(x))}")
    print(f"state_vector is {a.state_vector()}")
    # FIXME: Check the change is correct
```

    x = [2 1], y = [ 9. 13. 17.], loss = 27.0
    ∇𝑙𝑜𝑠𝑠(𝑦) = [-2.          1.          7.00000001]
    ∇𝑙𝑜𝑠𝑠(𝑥) = [15.00000002 33.00000001]
    backprop([-0.2  0.1  0.7]) = [1.5 3.3]
    Now a([2 1]) = [10.2 12.4 12.8], loss = 4.319999988809314
    state_vector is [0.4 0.8 0.6 3.2 3.9 4.3 6.2 6.9 7.3]

### Test batch operations

Test, for batch operations:
* input and output widths
* forward calculation

```python
if __name__ == '__main__':
    # Affine
    a = AffineLayer(2,3)
    a.set_state_from_vector(np.arange(9))

    # The input and output widths for the forward calculation are correct
    x = two_wide
    y = a(two_wide)
    assert y.shape[0] == x.shape[0]
    assert y.shape[1] == 3

    # The input and output widths for the backprop calculation are correct
    bp = a.backprop(three_wide * 0.001)
    assert bp.shape[0] == three_wide.shape[0]
    assert bp.shape[1] == x.shape[1]

    # The forward calculation is correct (in at least two instances)
    a.set_state_from_vector(np.arange(9))
    x = np.array([[0, 1],
                  [2, 3],
                  [4, 5],
                  [6, 7]])
    assert (a(x) == np.array([[ 9, 11, 13],
                              [15, 21, 27],
                              [21, 31, 41],
                              [27, 41, 55]])).all()
    #print(f"a.M is:\n{a.M}\na.b is {a.b}\nx is: {x}\ny is: {y}")
    a.set_state_from_vector(np.array([ 2,  3,  5,  7, 11, 13, 17, 19, 23]))
    y = a(x)
    #print(f"x is: {x}\ny is: {y}")
    assert (y == np.array([[ 24,  30,  36],
                           [ 42,  58,  72],
                           [ 60,  86, 108],
                           [ 78, 114, 144]])).all()
```

Test, for batch operations:
* back-propagation of the loss gradient
* learning (change in forward function) from the back-prop operation

```python
if __name__ == '__main__':
    # Affine
    a = AffineLayer(2,3)
    a.set_state_from_vector(np.arange(9))
    x = np.array([[0, 1],
                  [2, 3],
                  [4, 5],
                  [6, 7]])
    y = a(x)

    # It back-propagages the loss gradient

    # for loss function, use L2-distance from some ideal
    # (divided by 2, for convenient gradient = error)
    ideal = x @ arangep(2*3).reshape(2,3) + arangep(3,6) # A known, different parameter setting
    print(f"y - ideal =\n{y - ideal}")
    #loss = lambda v: (v - ideal).dot(v - ideal) / 2.0
    loss = lambda v: np.einsum('ij,ij', v-ideal, v-ideal) / (2 * v.shape[0])
    loss_at_y = loss(y)
    print(f"x =\n{x}\nideal =\n{ideal}\ny =\n{y}, loss = {loss_at_y}\n")

    # find numerical gradient of loss function at y, the layer output
    grad_y = VC.tensor_grad(loss, y)
    print(f"∇𝑙𝑜𝑠𝑠(𝑦) =\n{grad_y}")

    # find the numerical gradient of the loss w.r.t. the input of the layer
    grad_x = VC.tensor_grad(lambda x:loss(a(x)), x)
    print(f"∇𝑙𝑜𝑠𝑠(𝑥) =\n{grad_x}")

    # Back-propagate the loss gradient from layer output to input
    _ = a(x) # Make sure the last x is in the right place
    out_delE = grad_y * 0.01 # Backprop one percent of the loss gradient
    in_delE = a.backprop(out_delE)
    print(f"backprop({out_delE}) = {in_delE}")

    # The loss gradient back-propagated to the layer input is correct
    #assert closenuf(in_delE / 0.1, grad_x)

    # And how did the learning affect the layer?
    print(f"Now a({x}) = {a(x)}, loss = {loss(a(x))}")
    print(f"state_vector is {a.state_vector()}")
    # FIXME: Check the change is correct
```

    y - ideal =
    [[-15. -19. -23.]
     [-27. -37. -45.]
     [-39. -55. -67.]
     [-51. -73. -89.]]
    x =
    [[0 1]
     [2 3]
     [4 5]
     [6 7]]
    ideal =
    [[ 24  30  36]
     [ 42  58  72]
     [ 60  86 108]
     [ 78 114 144]]
    y =
    [[ 9. 11. 13.]
     [15. 21. 27.]
     [21. 31. 41.]
     [27. 41. 55.]], loss = 3765.5

    ∇𝑙𝑜𝑠𝑠(𝑦) =
    [[ -55.99999895  -81.33333358 -102.66666595]
     [ -47.99999988  -67.99999983  -83.99999933]
     [ -40.00000035  -54.66666653  -65.33333317]
     [ -31.99999992  -41.33333346  -46.66666655]]
    ∇𝑙𝑜𝑠𝑠(𝑥) =
    [[ -286.6666664  -1006.6666664 ]
     [ -235.9999994   -836.00000062]
     [ -185.33333332  -665.33333393]
     [ -134.66666633  -494.66666701]]
    backprop([[-0.55999999 -0.81333334 -1.02666666]
     [-0.48       -0.68       -0.83999999]
     [-0.4        -0.54666667 -0.65333333]
     [-0.32       -0.41333333 -0.46666667]]) = [[ -2.86666665 -10.06666661]
     [ -2.35999998  -8.35999996]
     [ -1.85333333  -6.65333333]
     [ -1.34666667  -4.94666666]]
    Now a([[0 1]
     [2 3]
     [4 5]
     [6 7]]) = [[ 16.99999999  21.93333333  26.06666661]
     [ 44.44        60.94666666  74.41333313]
     [ 71.88000001  99.95999999 122.75999966]
     [ 99.32000002 138.97333332 171.10666618]], loss = 325.0071954815595
    state_vector is [ 4.48000001  7.02666667  9.09333331  9.24       12.48       15.07999996
      7.75999999  9.45333333 10.98666665]

### Test batch operations when the affine layer has only one input

Test, for batch operations:
* input and output widths
* forward calculation

```python
if __name__ == '__main__':
    # Affine
    a = AffineLayer(1,3)
    a.set_state_from_vector(np.arange(6))

    # The input and output widths for the forward calculation are correct
    x = one_wide
    y = a(one_wide)
    assert y.shape[0] == x.shape[0]
    assert y.shape[1] == 3

    # The input and output widths for the backprop calculation are correct
    bp = a.backprop(three_wide * 0.001)
    assert bp.shape[0] == three_wide.shape[0]
    assert bp.shape[1] == x.shape[1]

    # The forward calculation is correct (in at least two instances)
    a.set_state_from_vector(np.arange(6))
    x = np.array([[0],
                  [1],
                  [2],
                  [3]])
    assert (y == np.array([[ 3.,  4.,  5.],
                           [ 3.,  5.,  7.],
                           [ 3.,  6.,  9.],
                           [ 3.,  7., 11.]])).all()
    #print(f"a.M is:\n{a.M}\na.b is {a.b}\nx is: {x}\ny is: {y}")
    a.set_state_from_vector(np.array([ 2,  3,  5,  7, 11, 13]))
    y = a(x)
    #print(f"x is: {x}\ny is: {y}")
    assert (a(x) == np.array([[ 7, 11, 13],
                              [ 9, 14, 18],
                              [11, 17, 23],
                              [13, 20, 28]])).all()
```

Test, for batch operations:
* back-propagation of the loss gradient
* learning (change in forward function) from the back-prop operation

```python
if __name__ == '__main__':
    # Affine
    a = AffineLayer(1,3)
    a.set_state_from_vector(np.arange(6))
    x = np.array([[0],
                  [1],
                  [2],
                  [3]])
    y = a(x)
    #print(f"x =\n{x}\ny =\n{y}")

    # It back-propagages the loss gradient

    # for loss function, use L2-distance from some ideal
    # (divided by 2, for convenient gradient = error)
    ideal = x @ arangep(1*3).reshape(1,3) + arangep(3,6) # A known, different parameter setting
    print(f"y - ideal =\n{y - ideal}")
    #loss = lambda v: (v - ideal).dot(v - ideal) / 2.0
    loss = lambda v: np.einsum('ij,ij', v-ideal, v-ideal) / (2 * v.shape[0])
    loss_at_y = loss(y)
    print(f"x =\n{x}\nideal =\n{ideal}\ny =\n{y}, loss = {loss_at_y}\n")

    # find numerical gradient of loss function at y, the layer output
    grad_y = VC.tensor_grad(loss, y)
    print(f"∇𝑙𝑜𝑠𝑠(𝑦) =\n{grad_y}")

    # find the numerical gradient of the loss w.r.t. the input of the layer
    grad_x = VC.tensor_grad(lambda x:loss(a(x)), x)
    print(f"∇𝑙𝑜𝑠𝑠(𝑥) =\n{grad_x}")

    # Back-propagate the loss gradient from layer output to input
    _ = a(x) # Make sure the last x is in the right place
    out_delE = grad_y * 0.01 # Backprop one percent of the loss gradient
    in_delE = a.backprop(out_delE)
    print(f"backprop({out_delE}) = {in_delE}")

    # The loss gradient back-propagated to the layer input is correct
    #assert closenuf(in_delE / 0.1, grad_x)

    # And how did the learning affect the layer?
    print(f"Now a({x}) = {a(x)}, loss = {loss(a(x))}")
    print(f"state_vector is {a.state_vector()}")
    # FIXME: Check the change is correct
```

    y - ideal =
    [[-14. -15. -18.]
     [-16. -17. -21.]
     [-18. -19. -24.]
     [-20. -21. -27.]]
    x =
    [[0]
     [1]
     [2]
     [3]]
    ideal =
    [[17 19 23]
     [19 22 28]
     [21 25 33]
     [23 28 38]]
    y =
    [[ 3.  4.  5.]
     [ 3.  5.  7.]
     [ 3.  6.  9.]
     [ 3.  7. 11.]], loss = 570.25

    ∇𝑙𝑜𝑠𝑠(𝑦) =
    [[-22.66666667 -26.00000005 -34.00000003]
     [-22.66666667 -24.66666677 -31.33333342]
     [-22.66666667 -23.33333344 -28.66666659]
     [-22.66666667 -22.         -25.99999993]]
    ∇𝑙𝑜𝑠𝑠(𝑥) =
    [[-94.0000001 ]
     [-87.3333334 ]
     [-80.66666669]
     [-74.0000001 ]]
    backprop([[-0.22666667 -0.26       -0.34      ]
     [-0.22666667 -0.24666667 -0.31333333]
     [-0.22666667 -0.23333333 -0.28666667]
     [-0.22666667 -0.22       -0.26      ]]) = [[-0.94      ]
     [-0.87333334]
     [-0.80666667]
     [-0.74      ]]
    Now a([[0]
     [1]
     [2]
     [3]]) = [[ 3.90666667  4.96        6.2       ]
     [ 5.26666667  7.33333334  9.86666666]
     [ 6.62666667  9.70666668 13.53333333]
     [ 7.98666667 12.08000001 17.19999999]], loss = 389.28084440511896
    state_vector is [1.36       2.37333334 3.66666666 3.90666667 4.96       6.2       ]

### Test multinet batch operations

Test, for batch operations:
* input and output widths
* state vector setting and getting
* forward calculation

```python
if __name__ == '__main__':
    # Affine
    a = AffineLayer(2,3, multiples=4)
    a.set_state_from_vector(np.arange(9*4))

    # The input and output widths for the forward calculation are correct
    x = two_wide
    y = a(two_wide)
    assert y.shape[-2] == x.shape[-2]
    assert y.shape[-1] == 3

    # The input and output widths for the backprop calculation are correct
    bp = a.backprop(three_wide * 0.001)
    assert bp.shape[-2] == three_wide.shape[-2]
    assert bp.shape[-1] == x.shape[-1]

    # The forward calculation is correct (in at least two instances)
    a.set_state_from_vector(np.arange(9*4))
    x = np.array([[0, 1],
                  [2, 3],
                  [4, 5],
                  [6, 7]])
    assert (a(x) == np.array(
      [[[ 27.,  29.,  31.],
        [ 33.,  39.,  45.],
        [ 39.,  49.,  59.],
        [ 45.,  59.,  73.]],

       [[ 36.,  38.,  40.],
        [ 66.,  72.,  78.],
        [ 96., 106., 116.],
        [126., 140., 154.]],

       [[ 45.,  47.,  49.],
        [ 99., 105., 111.],
        [153., 163., 173.],
        [207., 221., 235.]],

       [[ 54.,  56.,  58.],
        [132., 138., 144.],
        [210., 220., 230.],
        [288., 302., 316.]]])).all()
    #print(f"a.M is:\n{a.M}\na.b is {a.b}\nx is: {x}\ny is: {y}")
    a.set_state_from_vector(arangep(9*4))
    y = a(x)
    #print(f"x is: {x}\ny is: {y}")
    assert (y == np.array(
      [[[ 104.,  112.,  116.],
        [ 122.,  140.,  152.],
        [ 140.,  168.,  188.],
        [ 158.,  196.,  224.]],

       [[ 136.,  140.,  150.],
        [ 228.,  240.,  270.],
        [ 320.,  340.,  390.],
        [ 412.,  440.,  510.]],

       [[ 180.,  190.,  198.],
        [ 368.,  394.,  414.],
        [ 556.,  598.,  630.],
        [ 744.,  802.,  846.]],

       [[ 218.,  232.,  240.],
        [ 510.,  540.,  564.],
        [ 802.,  848.,  888.],
        [1094., 1156., 1212.]]])).all()
```

Test, for batch operations:
* back-propagation of the loss gradient
* learning (change in forward function) from the back-prop operation

```python
if __name__ == '__main__':
    # Affine
    a = AffineLayer(2,3, multiples=2)
    a.set_state_from_vector(np.arange(9*2))
    x = np.array([[0, 1],
                  [2, 3],
                  [4, 5],
                  [6, 7]])
    y = a(x)

    # It back-propagages the loss gradient

    # for loss function, use L2-distance from some ideal
    # (divided by 2, for convenient gradient = error)
    ideal = x @ arangep(2*2*3).reshape(2,2,3) + arangep(2*3, 2*2*3).reshape(2,1,3) # A known, different parameter setting
    print(f"y - ideal =\n{y - ideal}")
    #loss = lambda v: (v - ideal).dot(v - ideal) / 2.0
    loss = lambda v: np.einsum('...ij,...ij', y-ideal, y-ideal) / (2 * np.atleast_2d(y-ideal).shape[-2])
    loss_at_y = loss(y)
    print(f"x =\n{x}\nideal =\n{ideal}\ny =\n{y}, loss = {loss_at_y}\n")

    # find numerical gradient of loss function at y, the layer output
    grad_y = VC.tensor_grad(loss, y)
    print(f"∇𝑙𝑜𝑠𝑠(𝑦) =\n{grad_y}")

    # find the numerical gradient of the loss w.r.t. the input of the layer
    grad_x = VC.tensor_grad(lambda x:loss(a(x)), x)
    print(f"∇𝑙𝑜𝑠𝑠(𝑥) =\n{grad_x}")

    # Back-propagate the loss gradient from layer output to input
    _ = a(x) # Make sure the last x is in the right place
    out_delE = grad_y * 0.01 # Backprop one percent of the loss gradient
    in_delE = a.backprop(out_delE)
    print(f"backprop({out_delE}) = {in_delE}")

    # The loss gradient back-propagated to the layer input is correct
    #assert closenuf(in_delE / 0.1, grad_x)

    # And how did the learning affect the layer?
    print(f"Now a({x}) = {a(x)}, loss = {loss(a(x))}")
    print(f"state_vector is {a.state_vector()}")
    # FIXME: Check the change is correct
```

```python
%debug
```

    > <ipython-input-26-83d7fd585499>(15)<module>()
         13     # for loss function, use L2-distance from some ideal
         14     # (divided by 2, for convenient gradient = error)
    ---> 15     ideal = x @ arangep(2*2*3).reshape(2,2,3) + arangep(2*2*3, 2*(2*3 + 3)).reshape(2,1,3) # A known, different parameter setting
         16     print(f"y - ideal =\n{y - ideal}")
         17     #loss = lambda v: (v - ideal).dot(v - ideal) / 2.0

    ipdb>  arangep(2*2*3, 2*(2*3 + 3))

    array([ 67,  71,  73,  79,  83,  89,  97, 101, 103, 107, 109, 113])

    ipdb>  2*(2*3 + 3) - 2*2*3

    6

    ipdb>  2*(2*3 + 3), 2*2*3

    (18, 12)

    ipdb>  range(7,2)

    range(7, 2)

    ipdb>  list(range(7,2))

    *** Error in argument: '(range(7,2))'

    ipdb>  exit

### Test multinet batch operations when the affine layer has only one input

Test, for batch operations:
* input and output widths
* forward calculation

```python
if __name__ == '__main__':
    # Affine
    a = AffineLayer(1,3)
    a.set_state_from_vector(np.arange(6))

    # The input and output widths for the forward calculation are correct
    x = one_wide
    y = a(one_wide)
    assert y.shape[0] == x.shape[0]
    assert y.shape[1] == 3

    # The input and output widths for the backprop calculation are correct
    bp = a.backprop(three_wide * 0.001)
    assert bp.shape[0] == three_wide.shape[0]
    assert bp.shape[1] == x.shape[1]

    # The forward calculation is correct (in at least two instances)
    a.set_state_from_vector(np.arange(6))
    x = np.array([[0],
                  [1],
                  [2],
                  [3]])
    assert (y == np.array([[ 3.,  4.,  5.],
                           [ 3.,  5.,  7.],
                           [ 3.,  6.,  9.],
                           [ 3.,  7., 11.]])).all()
    #print(f"a.M is:\n{a.M}\na.b is {a.b}\nx is: {x}\ny is: {y}")
    a.set_state_from_vector(np.array([ 2,  3,  5,  7, 11, 13]))
    y = a(x)
    #print(f"x is: {x}\ny is: {y}")
    assert (a(x) == np.array([[ 7, 11, 13],
                              [ 9, 14, 18],
                              [11, 17, 23],
                              [13, 20, 28]])).all()
```

Test, for batch operations:
* back-propagation of the loss gradient
* learning (change in forward function) from the back-prop operation

```python
if __name__ == '__main__':
    # Affine
    a = AffineLayer(1,3)
    a.set_state_from_vector(np.arange(6))
    x = np.array([[0],
                  [1],
                  [2],
                  [3]])
    y = a(x)
    #print(f"x =\n{x}\ny =\n{y}")

    # It back-propagages the loss gradient

    # for loss function, use L2-distance from some ideal
    # (divided by 2, for convenient gradient = error)
    ideal = x @ arangep(1*3).reshape(1,3) + arangep(3,6) # A known, different parameter setting
    print(f"y - ideal =\n{y - ideal}")
    #loss = lambda v: (v - ideal).dot(v - ideal) / 2.0
    loss = lambda v: np.einsum('ij,ij', v-ideal, v-ideal) / (2 * v.shape[0])
    loss_at_y = loss(y)
    print(f"x =\n{x}\nideal =\n{ideal}\ny =\n{y}, loss = {loss_at_y}\n")

    # find numerical gradient of loss function at y, the layer output
    grad_y = VC.tensor_grad(loss, y)
    print(f"∇𝑙𝑜𝑠𝑠(𝑦) =\n{grad_y}")

    # find the numerical gradient of the loss w.r.t. the input of the layer
    grad_x = VC.tensor_grad(lambda x:loss(a(x)), x)
    print(f"∇𝑙𝑜𝑠𝑠(𝑥) =\n{grad_x}")

    # Back-propagate the loss gradient from layer output to input
    _ = a(x) # Make sure the last x is in the right place
    out_delE = grad_y * 0.01 # Backprop one percent of the loss gradient
    in_delE = a.backprop(out_delE)
    print(f"backprop({out_delE}) = {in_delE}")

    # The loss gradient back-propagated to the layer input is correct
    #assert closenuf(in_delE / 0.1, grad_x)

    # And how did the learning affect the layer?
    print(f"Now a({x}) = {a(x)}, loss = {loss(a(x))}")
    print(f"state_vector is {a.state_vector()}")
    # FIXME: Check the change is correct
```

## Network

### Network assembly

The simplest, the empty network, does identity:

```python
if __name__ == '__main__':
    net = Network()
    assert all(x == net(x) for x in [0, 42, 'cows in trouble'])
    assert all((x == net(x)).all() for x in [np.arange(7), np.arange(3*4*5).reshape(3,4,5)])
```

A stack of maps composes the operations:

```python
if __name__ == '__main__':
    net = Network()
    net.extend(MapLayer(lambda x: x**3, lambda d: 3*d**2))
    assert all(net(x) == x**3 for x in [0, 42, -3.14])
    net.extend(MapLayer(lambda x: 7-x, lambda d: -1))
    assert all(net(x) == 7-x**3 for x in [0, 42, -3.14])

    # It operates on each element of an input vector separately
    assert (net(np.arange(4)) == 7 - np.arange(4) ** 3).all()
```

A composition of affine transformations

_[to do someday]_

### Network Learning

Test simple batch learning of a single affine layer

```python
if __name__ == '__main__':
    from pprint import pprint

    net = Network()
    a = AffineLayer(2,3)
    a.set_state_from_vector(np.arange(9)) # A well-known initial state
    net.extend(a)
    print(f"\nNet has state {net.state_vector()}")

    x = np.array([[0, 1],
                  [2, 3],
                  [4, 5],
                  [6, 7]])

    # The net wraps the layer
    y = a(x)
    assert (net(x) == y).all()

    # Make the training batch.
    # We use a separate affine layer, initialized differently, to determine the ideal
    t_a = AffineLayer(2,3)
    t_a.set_state_from_vector(arangep(9)) # A known different initial state (of primes)
    ideal = t_a(x)

    fact = (x, ideal)
    print(f"fact is:")
    pprint(fact, indent=1)
    print(f"net(x) =\n{net(x)}")

    net.eta = 0.01
    for i in range(10):
        print(f"net.learn([fact]) = {net.learn([fact])}")
    print(f"net(x) =\n{net(x)}")

    # A simple fact yielder. Since it delivers multiple facts in succession,
    # it is a "facts", aka "batch cluster"
    def fact_ory(fact, n):
        for i in range(n):
            yield fact

    def facts_printer(facts):
        # after Network.learn
        for fact in facts:
            print(f"fact: ")
            pprint(fact)
            x, ideal = fact
            print(f"from which we get:\n\tx={x}\n\tideal={ideal}\n")

    #facts_printer(fact_ory(fact, 5))

    #print(f"list(fact_ory(facts[0], 3)) =\n{list(fact_ory(facts[0], 3))}\n")
    print(f"net.learn(fact_ory(fact,10)) = {net.learn(fact_ory(fact,10))}")
    print(f"net(x) =\n{net(x)}")
    for i in range(1000):
        loss = net.learn(fact_ory(fact,10))
        if loss < 1e-7:
            break
    print(f"did {(i+1)*10} more learnings of fact. Now loss is {loss}")
    print(f"net(x) =\n{net(x)}")

    print(f"net.state_vector() = {net.state_vector()}")

    # The network has learned the target transform
    assert closenuf(net(x), fact[1])

if False: # This section of the test is misconceived. Skip it
    # Save prior results and learn again, with different batch clustering
    prev_run_loss = loss
    prev_y = net(x)
    net.set_state_from_vector(np.arange(9)) # A well-known initial state
    print(f"\nReset net to state {net.state_vector()}")

    # Try multiple batches in each call to Network.learn
    def multibatch_fact_ory(fact, n):
        for i in range(n//2):
            yield fact, fact
    facts_printer(multibatch_fact_ory(fact, 5))

    for i in range(1000):
        loss = net.learn(multibatch_fact_ory(fact,10))
        if loss < 1e-7:
            break
    print(f"did {(i+1)*10} learnings of fact. Now loss is {loss}")
    print(f"net(x) =\n{net(x)}")

    # The results should match exactly
    assert loss == prev_run_loss
    assert (net(x) == prev_y).all()
```

### Test Network.losses

```python
if __name__ == '__main__':
    # Make a network. Leave it with the default identity behavior.
    net = Network()

    x = np.array([[0, 1],
                  [2, 3],
                  [4, 5],
                  [6, 7]])
    ideal = net(x)
    facts = [(x, ideal), (x, ideal-np.array([1,-1])), (x, 2*x)]
    assert (net.losses(facts) == [0, 1, 17.5])

    # Add some layers
    net.extend(AffineLayer(2,3)).extend(MapLayer(np.sin, np.cos)).extend(AffineLayer(3,2))
    # Place it in a known state for test repeatability
    net.set_state_from_vector(np.arange(len(net.state_vector())))
    ideal = net(x)
    facts = [(x, ideal), (x, ideal-np.array([1,-1]))]
    #print(net.losses(facts))
    assert (net.losses(facts) == [0, 1])
```

---

# Publishing

To produce an importable `nnn.py`:
1. Save this notebook
1. Uncomment the `jupyter nbconvert` line below
1. Execute it.
1. Comment out the convert again
1. Save the notebook again in that form

```python
###!jupyter nbconvert --to script nnn.ipynb
```

```python

```

max_stars_repo_name: pramasoul/aix

max_issues_repo_name: pramasoul/aix

max_forks_repo_name: pramasoul/aix


text:
# Sympy

```python
import sympy as sp

# See: http://docs.sympy.org/latest/tutorial/printing.html
sp.init_printing()
```

## Live shell

http://live.sympy.org/

## Make symbols

Make one symbol:

```python
x = sp.symbols("x")
```

Make several symbols at once:

```python
x, y, z = sp.symbols("x y z")
```

## Substitute

See: http://docs.sympy.org/latest/tutorial/basic_operations.html#substitution

```python
x = sp.symbols("x")
```

```python
expr = sp.cos(x) + 1
expr.subs(x, 0)
```

## Simplify

See: http://docs.sympy.org/latest/tutorial/simplification.html#simplify

```python
x = sp.symbols("x")
```

```python
sp.simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))
```

```python
sp.simplify(sp.exp(sp.I * sp.pi) + sp.exp(sp.I * sp.pi))
```

## Factor

See: http://docs.sympy.org/latest/tutorial/simplification.html#factor

```python
x = sp.symbols("x")
```

```python
sp.factor(x**3 - x**2 + x - 1)
```

## Expand

See: http://docs.sympy.org/latest/tutorial/simplification.html#expand

```python
x = sp.symbols("x")
```

```python
sp.expand((x - 1)*(x**2 + 1))
```

## Solve

See: http://docs.sympy.org/latest/tutorial/solvers.html

```python
x = sp.symbols("x")
```

```python
eq = sp.Eq(x**2, 1)
```

```python
sp.solveset(eq, x)
```

## Derivatives

See: http://docs.sympy.org/latest/tutorial/calculus.html#derivatives

```python
x, y, z = sp.symbols("x y z")
```

### Create an unevaluated derivative

```python
sp.Derivative(sp.cos(x), x)
```

### Evaluate an unevaluated derivative

```python
diff = sp.Derivative(sp.cos(x), x)
diff
```

```python
diff.doit()
```

### Directly compute an integral

```python
sp.diff(sp.cos(x), x)
```

### Print the equation

```python
expr = sp.exp(x*y*z)
diff = sp.Derivative(expr, x, y, y, z, z, z, z)
sp.Eq(diff, diff.doit())
```

### First derivatives

```python
diff = sp.Derivative(sp.cos(x), x)
sp.Eq(diff, diff.doit())
```

```python
diff = sp.Derivative(3*sp.cos(x)**2, x)
sp.Eq(diff, diff.doit())
```

```python
diff = sp.Derivative(sp.exp(x**2), x)
sp.Eq(diff, diff.doit())
```

### Second derivatives

```python
diff = sp.Derivative(x**4, x, 2)
sp.Eq(diff, diff.doit())
```

or

```python
diff = sp.Derivative(x**4, x, x)
sp.Eq(diff, diff.doit())
```

### Third derivatives

```python
diff = sp.Derivative(x**4, x, 3)
sp.Eq(diff, diff.doit())
```

or

```python
diff = sp.Derivative(x**4, x, x, x)
sp.Eq(diff, diff.doit())
```

### Derivatives with respect to several variables at once

```python
diff = sp.Derivative(sp.exp(x*y), x, y)
sp.Eq(diff, diff.doit())
```

### Multiple derivatives with respect to several variables at once

```python
diff = sp.Derivative(sp.exp(x*y*z), x, y, y, z, z, z, z)
sp.Eq(diff, diff.doit())
```

## Integrals

See: http://docs.sympy.org/latest/tutorial/calculus.html#integrals

```python
x, y, z = sp.symbols("x y z")
```

### Create an unevaluated derivative

```python
sp.Integral(sp.cos(x), x)
```

### Evaluate an unevaluated derivative

```python
integ = sp.Integral(sp.cos(x), x)
integ
```

```python
integ.doit()
```

### Directly compute an integral

```python
sp.integrate(sp.cos(x), x)
```

### Print the equation

```python
integ = sp.Integral(sp.cos(x), x)
sp.Eq(integ, integ.doit())
```

### Create an indefinite integral (i.e. an antiderivative or primitive)

```python
integ = sp.Integral(sp.cos(x), x)
sp.Eq(integ, integ.doit())
```

### Create a definite integral

`sp.oo` means infinity.

```python
integ = sp.Integral(sp.cos(x), (x, -sp.oo, sp.oo))
sp.Eq(integ, integ.doit())
```

```python
integ = sp.Integral(sp.cos(x), (x, -sp.pi, sp.pi))
sp.Eq(integ, integ.doit())
```

```python
integ = sp.Integral(sp.exp(-x), (x, 0, sp.oo))
sp.Eq(integ, integ.doit())
```

### Multiple integrals

```python
integ = sp.Integral(sp.cos(x), (x, -sp.oo, sp.oo), (x, -sp.oo, sp.oo))
sp.Eq(integ, integ.doit())
```

### Multiple variables integrals

```python
integ = sp.Integral(sp.cos(x**2 + y**2), (x, -sp.oo, sp.oo), (y, -sp.oo, sp.oo))
sp.Eq(integ, integ.doit())
```

## Limits

See: http://docs.sympy.org/latest/tutorial/calculus.html#limits

## Series expansion

See: http://docs.sympy.org/latest/tutorial/calculus.html#series-expansion

## Finite differences

See: http://docs.sympy.org/latest/tutorial/calculus.html#finite-differences

max_stars_repo_name: jdhp-docs/python-notebooks

max_issues_repo_name: jdhp-docs/python-notebooks

max_forks_repo_name: jdhp-docs/python-notebooks