Write a program that meets the following requirements:

It constructs a NumPy matrix A to represent the following mathematical matrix: \[ A = \begin{bmatrix} 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 \\ 10 & 11 & 12 & 13 & 14 & 15 & 16 & 17 & 18 & 19 \end{bmatrix}. \]
It defines a function left_up_sum(A: np.ndarray) -> np.ndarray that adds the component (element) to the left and the component above (wrapping, if necessary) to each component. The function should pass through the array once, row-by-row, and return a new array. The function should be able to handle any size of matrix.
It defines a function left_up_sums(A: np.ndarray, n: int) -> np.ndarray that executes left_up_sum() n times and returns a new array.
It calls left_up_sums() on A and prints the returned array for the following values of n: 0, 1, 4.

The inner product of two real $n$-vectors $\vec{x}$ and $\vec{y}$ is defined as \[ \langle \vec{x}, \vec{y} \rangle = \sum_{i=0}^{n-1} x_i y_i. \] The result is a scalar. The np.inner() and np.dot() functions can be used in NumPy to find the inner product of two vectors of the same size. In this problem, we will write our own function that computes the real inner product even if they are of different sizes. Write a program that meets the following requirements:

It defines a function
```
inner_flat_trunc(x: np.ndarray, y: np.ndarray) -> float
```
that returns the truncated inner product of vectors a and b even if the sizes of the vectors do not match by using a truncated version of the one that is too long. The function should handle any shape of input arrays by using the flatten() method before truncating and taking the inner product. If both input arrays do not have dtype attribute np.dtype('float') or np.dtype('int'), the function should raise a TypeError exception.
It calls inner_flat_trunc() on the following arrays:
1. A pair of arrays from the lists:
  [-1.1, 3, 2.9, -1, -9.2, 0.1] and [1.3, 0.2, 8.3]
2. An array of the integers from $0$ through $13$ and an array of the integers from $3$ through $12$
3. An array of $21$ linearly spaced elements from $0$ through $10$ and an array of $11$ linearly spaced elements from $5$ through $25$.
4. A pair of arrays of elements from the lists [True, False, True] and [0, 1, 2] (handle the exception in the main script so it runs without raising the exception)

Consider the following mathematical matrices and vectors: $$\begin{align} A = \begin{bmatrix} 2 & 1 & 9 & 0 \\ 0 & -1 & -2 & 3 \\ -3 & 0 & 8 & -4 \end{bmatrix} \quad B = \begin{bmatrix} 0 & 9 & -1 \\ 1 & 0 & 3 \\ 0 & -1 & 1 \end{bmatrix} \quad \vec{x} = \begin{bmatrix} 1 \\ 0 \\ -1 \end{bmatrix} \quad \vec{y} = \begin{bmatrix} 3 & 0 & -1 \end{bmatrix}. \end{align}$$

Write a program that meets the following requirements:

It defines NumPy arrays to represent $A$, $B$, (column vector) $\vec{x}$, and (row vector) $\vec{y}$.
It computes and prints the following quantities:
1. $BA$
2. $A^\top B - 6 J_{4\times 3}$, where $J_{4\times 3}$ is the $4\times 3$ matrix of all $1$ components
3. $B\vec{x} + \vec{y}^\top$
4. $\vec{x}\vec{y} + B$
5. $\vec{y}\vec{x}$
6. $\vec{y}B^{-1}\vec{x}$
7. $C B$, where $C$ is the $3\times 3$ submatrix of the first three columns of $A$

The following program meets the requirements:

""""Solution to Chapter 3 problem 3H"""
import numpy as np

# %% [markdown]
## Introduction
# This program defines matrices A and B and vectors x and y. It then
# computes several quantities that require matrix operations.
# %% [markdown]
## a. Define Arrays
# %%
A = np.array(
    [
        [2, 1, 9, 0],
        [0, -1, -2, 3],
        [-3, 0, 8, -4],
    ]
)
B = np.array(
    [
        [0, 9, -1],
        [1, 0, 3],
        [0, -1, 1],
    ]
)
x = np.array([[1], [0], [-1]])
y = np.array([[3, 0, -1]])
# %% [markdown]
## b. Compute and Print Quantitites

### i. $B A$
# %%
print("B A = \n", B @ A)
# %% [markdown]
### ii. $A^\top B - 6 J_{4\times 3}$
# %%
print("A.T B - 6 J_4x3 = \n", A.T @ B - 6 * np.ones((4, 3)))
# %% [markdown]
### iii. $B \bm{x} + \bm{y}^\top$
# %%
print("B x + y.T = \n", B @ x + y.T)
# %% [markdown]
### iv. $\bm{x} \bm{y} + B$
# %%
print("x y + B = \n", x @ y + B)
# %% [markdown]
### v. $\bm{y} \bm{x}$
# %%
print("y x = \n", y @ x)
# %% [markdown]
### vi. $\bm{y} B^{-1} \bm{x}$
# %%
print("y B^-1 x = \n", y @ np.linalg.inv(B) @ x)
# %% [markdown]
### vii. $C B$
# %%
C = A[:, :3]  # First 3 columns of A (view)
print("C B = \n", C @ B)

# %% tags=["active-py"]
import sys

sys.path.append("../")
import engcom.engcom as engcom

pub = engcom.Publication(title="Problem 3H", author="Rico Picone")
pub.write(to="md")

This program prints the following to the console:

B A = 
 [[  3  -9 -26  31]
 [ -7   1  33 -12]
 [ -3   1  10  -7]]
A.T B - 6 J_4x3 = 
 [[ -6.  15. -11.]
 [ -7.   3. -10.]
 [ -8.  67. -13.]
 [ -3.  -2.  -1.]]
B x + y.T = 
 [[ 4]
 [-2]
 [-2]]
x y + B = 
 [[ 3  9 -2]
 [ 1  0  3]
 [-3 -1  2]]
y x = 
 [[4]]
y B^-1 x = 
 [[10.]]
C B = 
 [[  1   9  10]
 [ -1   2  -5]
 [  0 -35  11]]

Consider the array:

a = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8], [9, 10, 11]])  # 4x3

Write a program that performs and prints the results of the following operations on a without using for loops:

Adds 1 to all elements
Adds 1 to the last column
Flattens a to a vector
Reshapes a into a $3\times 4$ matrix
Adds the vector [1, 2, 3] to each row
Adds the vector [1, 2, 3, 4] to each column
Reshapes a to a column vector
Reshapes a to a row vector

The following program meets the requirements:

""""Solution to Chapter 3 problem DI"""
import numpy as np

# %% [markdown]
## Introduction
# This program defines a NumPy array and prints the results of
# several operations.
# %% [markdown]
## Define Array
# %%
a = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8], [9, 10, 11]])  # 4x3
# %% [markdown]
## a. Adds 1 to all elements
print("Adds 1 to all elements:\n", a + 1)
# %% [markdown]
## b. Adds 1 to the last column
b = a.copy()
b[:, -1] = b[:, -1] + 1
print("Adds 1 to the last column:\n", b)
# %% [markdown]
## c. Flattens a to a vector
print("Flattens a to a vector\n", a.flatten())
# %% [markdown]
## d. Reshapes a into a 3×4 matrix
print("Reshapes a into a 3×4 matrix:\n", a.reshape((3, 4)))
# %% [markdown]
## e. Adds the vector [1, 2, 3] to each row
c = np.array([1, 2, 3])
cr = c.reshape((1, 3))  # Reshape for broadcasting
print("Adds the vector [1, 2, 3] to each row:\n", a + cr)
# %% [markdown]
## f. Adds the vector [1, 2, 3, 4] to each column
d = np.array([1, 2, 3, 4])
dr = d.reshape((4, 1))  # Reshape for broadcasting
print("Adds the vector [1, 2, 3, 4] to each column:\n", a + dr)
# %% [markdown]
## g. Reshapes a to a column vector
print("Reshapes a to a column vector:\n", a.reshape((a.size, 1)))
# %% [markdown]
## h. Reshapes a to a row vector
print("Reshapes a to a row vector:\n", a.reshape((1, a.size)))

This program prints the following to the console:

Adds 1 to all elements:
 [[ 1  2  3]
 [ 4  5  6]
 [ 7  8  9]
 [10 11 12]]
Adds 1 to the last column:
 [[ 0  1  3]
 [ 3  4  6]
 [ 6  7  9]
 [ 9 10 12]]
Flattens a to a vector
 [ 0  1  2  3  4  5  6  7  8  9 10 11]
Reshapes a into a 3×4 matrix:
 [[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
Adds the vector [1, 2, 3] to each row:
 [[ 1  3  5]
 [ 4  6  8]
 [ 7  9 11]
 [10 12 14]]
Adds the vector [1, 2, 3, 4] to each column:
 [[ 1  2  3]
 [ 5  6  7]
 [ 9 10 11]
 [13 14 15]]
Reshapes a to a column vector:
 [[ 0]
 [ 1]
 [ 2]
 [ 3]
 [ 4]
 [ 5]
 [ 6]
 [ 7]
 [ 8]
 [ 9]
 [10]
 [11]]
Reshapes a to a row vector:
 [[ 0  1  2  3  4  5  6  7  8  9 10 11]]

Write vectorized Python functions that operate element-wise on array arguments for the following mathematical functions:

$f(x) = x^2 + 3 x + 9$
$g(x) = 1 + \sin^2 x$
$h(x, y) = e^{-3x} + \ln y$
$F(x, y) = \left\lfloor x/y \right\rfloor$
$G(x, y) = \begin{cases} x^2 + y^2 & \text{if $x > y$} \\ 2 x & \text{otherwise} \end{cases}$

The following program meets the requirements:

""""Solution to Chapter 3 problem QX"""
import numpy as np

# %% [markdown]
## Introduction
# This program defines several mathematical functions as vectorized
# functions that can handle NumPy array inputs.
# %% [markdown]
## a. $f(x) = x^2 + 3 x + 9$
# %%
def f(x: np.ndarray) -> np.ndarray:
    return x**2 + 3 * x + 9


# %% [markdown]
## b. $g(x) = 1 + \sin^2 x$
# %%
def g(x: np.ndarray) -> np.ndarray:
    return 1 + np.sin(x) ** 2


# %% [markdown]
## c. $h(x, y) = e^{-3 x} + \ln y$
# %%
def h(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    return np.exp(-3 * x) + np.log(y)


# %% [markdown]
## d. $F(x, y) = \left\lfloor x/y \right\rfloor$
# %%
def F(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    return np.floor(x / y)


# %% [markdown]
## e. $G(x, y) = \begin{cases} x^2 + y^2 & \text{if $x > y$} \\ 2 x & \text{otherwise} \end{cases}$
# %%
def G(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    return np.where(x > y, x**2 + y**2, 2 * x)


# %% [markdown]
## Call Functions and Print
# %%
functions_args = (
    (f, 1),
    (g, 1),
    (h, 2),
    (F, 2),
    (G, 2),
)  # (fun, nargs)
x = np.array([1, 5, 10, 20, 30])
y = np.array([2, 7, 5, 10, 30])
print(f"x = {x}\ny = {y}")

for function_args in functions_args:
    if function_args[1] == 1:
        printable = np.array2string(function_args[0](x), precision=3)
        print(f"{function_args[0].__name__}(x) =", printable)
    elif function_args[1] == 2:
        printable = np.array2string(function_args[0](x, y), precision=3)
        print(f"{function_args[0].__name__}(x, y) =", printable)

This program prints the following to the console:

x = [ 1  5 10 20 30]
y = [ 2  7  5 10 30]
f(x) = [ 13  49 139 469 999]
g(x) = [1.708 1.92  1.296 1.833 1.976]
h(x, y) = [0.743 1.946 1.609 2.303 3.401]
F(x, y) = [0. 0. 2. 2. 1.]
G(x, y) = [  2  10 125 500  60]

Write a program that graphs each of the following functions over the specified domain:

$f(x) = \tanh(4 \sin x)$ for $x \in [-5, 8]$
$g(x) = \sin\sqrt{x}$ for $x \in [0, 100]$
$h(x) = \begin{cases} 0 & \text{if $x < 0$} \\ e^{-x} \sin(2 \pi x) & \text{otherwise} \end{cases}$ for $x \in [-2, 6]$

""""Solution to Chapter 3 problem DN"""
import numpy as np
import matplotlib.pyplot as plt

Introduction

This program defines several mathematical functions as vectorized functions that can handle NumPy array inputs and plots them over the given domain using Matplotlib. # Define Mathematical Functions Define $f(x) = x^2 + 3 x + 9$:

def f(x: np.ndarray) -> np.ndarray:
    return np.tanh(4 * np.sin(x))

Define $g(x) = 1 + \sin^2 x$:

def g(x: np.ndarray) -> np.ndarray:
    return np.sin(np.sqrt(x))

Define $h(x, y) = e^{-3 x} + \ln y$:

def h(x: np.ndarray) -> np.ndarray:
    return np.where(x >= 0, np.exp(-x) * np.sin(2 * np.pi * x), 0)

Plotting

Define a plotting function:

def plotter(fig, fun, limits, labels):
    x = np.linspace(limits[0], limits[1], 201)
    fig.gca().plot(x, fun(x))
    fig.gca().set_xlabel(labels[0])
    fig.gca().set_ylabel(labels[1])
    return fig

Plot $f(x)$:

fig, ax = plt.subplots()
plotter(fig, fun=f, limits=(-5, 8), labels=("$x$", "$f(x)$"))

Plot $g(x)$:

fig, ax = plt.subplots()
plotter(fig, fun=g, limits=(0, 100), labels=("$x$", "$g(x)$"))

Plot $h(x)$:

fig, ax = plt.subplots()
plotter(fig, fun=h, limits=(-2, 6), labels=("$x$", "$h(x)$"))

plt.show()

Write a program that loads and plots ideal gas data with the engcom.data.ideal_gas() function in the following way:

The data it loads is over the volume domain: $V \in [0.1, 2.1]$ m$^3$
The data it loads has $3$ temperatures: $V = 300, 400, 500$ K
It plots in a single graphic $P$ versus $V$ for each of the three temperatures
Each data point should be marked with a dot $\bullet$
Sequential data points should be connected by straight lines
Each plot should be labeled with its corresponding temperature, either next to the plot or in a legend

""""Solution to Chapter 3 problem WF"""
import numpy as np
import matplotlib.pyplot as plt
import engcom.data

Introduction

This program loads and plots ideal gas data with the ideal_gas() function from the engcom.data module. # Load the Data

d = engcom.data.ideal_gas(
    V=np.linspace(0.1, 2.1, 16),  # (m^3) Volume values
    T=np.array([300, 400, 500]),  # (K) Temperature values
)

Plot

Now d is a dictionary with the following key–value pairs:

"volume": $V_{16\times 1}$ (m$^3$)
"temperature": $T_{1\times 3}$ (K)
"pressure": $P_{16\times 3}$ (Pa)

We would like to plot $P$ versus $V$ for each of the $3$ temperatures $T$; that is, plot a sequence of pairs $(P_i, V_i)$ for each $T_j$. The following code loops through the temperatures and plots to the same axes object:

fig, ax = plt.subplots()
for j, Tj in enumerate(d["temperature"].flatten()):
    x = d["volume"]  # (m^3)
    y = d["pressure"][:, j] / 1e6  # (MPa)
    ax.plot(x, y, marker="o", label=f"$T = {Tj}$ K")

Finally, we label the axes and display the figure with the following code:

ax.set_xlabel("Volume (m$^3$)")
ax.set_ylabel("Pressure (MPa)")
ax.legend()  # Show labels in legend
plt.show()

Use the data from problem 3.7 to write a program that meets the following requirements:

It loads the pressure-volume-temperature data from problem 3.7.
It estimates the work $W$ done by the gas for each of the three values of temperatures via the integral equation \[ W = -\int_{0.1}^{2.1} P(V)\ dV. \] Note: An integral can be estimated from discrete data via the trapezoidal rule, which can be executed with NumPy’s np.trapz() function.
It generates a bar chart comparing the three values of work (one for each temperature).

""""Solution to Chapter 3 problem Y1"""
import numpy as np
import matplotlib.pyplot as plt
import engcom.data

Introduction

This program loads ideal gas data with the engcom.data.ideal_gas() function. It computes the work done by the gas over the given volume. Finally, it charts the work for each temperature. # Load the Data Load the pressure-volume-temperature data as follows:

d = engcom.data.ideal_gas(
    V=np.linspace(0.1, 2.1, 16),  # (m^3) Volume values
    T=np.array([300, 400, 500]),  # (K) Temperature values
)

Define and Compute the Work

Define a function to compute the work and apply the function.

def W(P: np.ndarray, V: np.ndarray, axis: int = -1) -> float:
    """Use trapezoidal integration to estimate the work from P(V)"""
    return np.trapz(P, V, axis=0)
work = W(P=d["pressure"], V=d["volume"])

Plot

Create a bar chart of work for each temperature

x = np.arange(len(work))
labels = np.char.add(d["temperature"].flatten().astype(dtype=str), " K")
fig, ax = plt.subplots()
ax.bar(x, work / 1e6)
ax.set_xticks(x, labels=labels)
ax.set_ylabel("Work (MJ)")
plt.show()

Write a program to bin data and create histogram charts that meets the following requirements:

It defines a function
```
binner(A: np.ndarray, nbins: int) -> (np.ndarray, np.ndarray)
```
that accepts an array A of data and returns an array for the frequency of the data in each bin and an array of the bin edges. Consider the following details:
1. The bin edges should include the left edge and not the right edge, except the rightmost, which should include the right edge (“left” and “right” here mean lesser and greater).
2. The bins should be of equal width.
3. Give a default value (e.g., 10) for the nbins argument.
4. Do not use the (nice) functions np.histogram() or plt.hist() for this exercise.
It defines a function histogram(A: np.ndarray, nbins: int) that calls binner() and plt.bar() to generate a histogram chart.
It loads all of the thermal conductivity data from the engcom.data module with the engcom.data.thermal_conductivity() function.
It generates $3$ histograms, one for each of the following material categories (key): "Metals", "Liquids", and "Gases". Be sure to properly label the axes.

Import packages:

import numpy as np
import matplotlib.pyplot as plt
import engcom.data

We begin by writing the binner() function:

def binner(A: np.ndarray, nbins: int = 10) -> (np.ndarray, np.ndarray):
    """Accepts an array A of data and returns an array for the frequency
    of the data in each bin and an array of the bin edges
    """
    Af = np.array(A).flatten()  # Ensure a vector
    minA = Af.min()
    maxA = Af.max()
    bin_edges = np.linspace(minA, maxA, nbins + 1)
    bw = bin_edges[1] - bin_edges[0]  # Bin width
    A_from_min = Af - minA  # Distance from min
    Abin = np.floor(A_from_min / bw)  # Bin index of each element of A
    Abin[Abin == nbins] = (
        nbins - 1
    )  # This moves element at maxA into the proper (last) bin
    nonempty_bin_index, nonempty_bin_frequency = np.unique(
        Abin, return_counts=True
    )  # Identifies nonempty bin indices and the frequency in each bin
    bin_frequency = np.zeros(nbins)  # Initialize frequencies to zero
    bin_frequency[
        nonempty_bin_index.astype(int)
    ] = nonempty_bin_frequency
    return bin_frequency, bin_edges

Now the histogram function can be written:

def histogram(A: np.ndarray, nbins: int = 10, xlabel=None, ylabel=None):
    """Create a histogram figure."""
    bin_frequency, bin_edges = binner(A, nbins=nbins)
    bw = bin_edges[1] - bin_edges[0]  # Bin width
    fig, ax = plt.subplots()
    ax.bar(
        x=bin_edges[0:-1],
        height=bin_frequency,
        width=0.95 * bw,
        align="edge",  # Left edge at x coordinate (left bin edge)
    )
    if xlabel is not None:
        ax.set_xlabel(xlabel)
    if ylabel is not None:
        ax.set_ylabel(ylabel)
    return fig

Now we load the thermal conductivity data. The material labels (keys) aren’t necessary for the histograms, so we just extract the values. The data is loaded as follows:

data = engcom.data.thermal_conductivity()
tc_metals = list(data["Metals"].values())
tc_liquids = list(data["Liquids"].values())
tc_gases = list(data["Gases"].values())

Now call the histogram() function to generate a histogram chart for each data set.

fig = histogram(
    tc_metals,
    nbins=10,
    xlabel="Thermal conductivity (W/(m$\cdot$K))",
    ylabel="Frequency (metals)",
)

fig = histogram(
    tc_liquids,
    nbins=10,
    xlabel="Thermal conductivity (W/(m$\cdot$K))",
    ylabel="Frequency (liquids)",
)

fig = histogram(
    tc_gases,
    nbins=10,
    xlabel="Thermal conductivity (W/(m$\cdot$K))",
    ylabel="Frequency (gases)",
)
plt.show()

You will now create life. John Conway’s Game of Life is a cellular automata game that explores the notion of life. In this problem, you will write a program for the game, which is played on a 2D grid. The grid is composed of elements called cells, each of which can be either alive or dead at a given moment. The rules of the game are simple (Johnston and Greene 2022):

If a cell is alive, it survives to the next generation if it has 2 or 3 live neighbors; otherwise it dies.
If a cell is dead, it comes to life in the next generation if it has exactly $3$ live neighbors; otherwise it stays dead.

The neighbors of a cell are those eight cells adjacent to it (including diagonals).

Write a program for playing the game of life that meets the following requirements:

It defines a function
```
game_of_life(A: np.ndarray)
```
that accepts a matrix A that encodes the starting state for the game. Use 1 to signify an alive cell and 0 to signify a dead cell. Consider the following details:
1. The game is traditionally played on an infinite grid. However, your program should play the game of life on a torus (doughnut) made from sewing the opposite edges of the starting state A grid. For instance, the neighbors above a cell in the top row are on the bottom row (i.e., neighbors wrap).
2. A visualization is required. A very useful Matplotlib function here is plt.matshow(A), which will display the numerical values of a matrix in a grid. For instance, try the following:
```
plt.matshow([[0,1,1],[1,0,1],[0,0,1]])
```
3. Strongly consider using additional functions to define operations like “evolve one generation,” “kill,” “animate,” and “visualize.”
It calls game_of_life() on matrices corresponding to the following starting states:
1. A $5 \times 5$ grid of cells with the following pattern (blinker): \[ \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0\end{bmatrix} \]
2. A $20\times 20$ grid of cells, all dead ($0$) except a group near the center with the following pattern (glider): \[ \begin{bmatrix}0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 1 & 1\end{bmatrix} \]
3. A $40\times 40$ grid of cells, all dead ($0$) except a group near the center with the pattern that can be loaded as a list from the engcom.data module with the function call
```
engcom.data.game_of_life_starts("gosper_glider")
```

Import packages:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import engcom.data

Functions for Playing the Game

We begin by writing the neighbor_indices() function to get the indices of the neighbors of a cell:

def neighbor_indices(A: np.ndarray, index: tuple) -> list:
    """Returns a list of tuple indices of A that are neighbors of index
    """
    I = np.array(index)  # Index as an array, to treat as a coordinate
    # Define coordinate shifts (x rows, y columns)
    L = np.array([-1, 0])  # Left
    R = np.array([1, 0])  # Right
    U = np.array([0, 1])  # Up
    D = np.array([0, -1])  # Down
    LU = L + U  # Left and up
    LD = L + D  # Left and down
    RU = R + U  # Right and up
    RD = R + D  # Right and down
    shifts = [L, R, U, D, LU, LD, RU, RD]
    # Prewrapped neighbors
    neighbors_pre = []
    for shift in shifts:
        neighbors_pre.append(I + shift)
    # Wrap the neighbors
    x_size, y_size = A.shape
    neighbors = []
    for neighbor in neighbors_pre:
        x, y = neighbor
        neighbors.append((
            x % x_size,  # Wrapped x
            y % y_size,  # Wrapped y
        ))  # Tuples, wrapped
    return neighbors

Next is to write a function to determine if a given cell is alive:

def is_alive(A: np.ndarray, index: tuple) -> bool:
    """Returns True if A[index] == 1; returns False, otherwise"""
    if A[*index] == 1:
        return True
    else:
        return False

Now we can write a function to determine the number of alive neighbors a cell has:

def n_alive_neighbors(A: np.ndarray, index: tuple) -> int:
    """Returns the number of alive neighbors for a tuple index"""
    neighbors = neighbor_indices(A, index)
    n_alive = 0
    for neighbor in neighbors:
        if is_alive(A, neighbor):
            n_alive += 1
    return n_alive

Now define a function to find the next value of a cell, given its current state and those of its neighbors and Life rules:

def evolve_cell(A: np.ndarray, index: tuple):
    """Returns the next value of A[index] based on neighbor rules"""
    n_alive_neighbors_ = n_alive_neighbors(A, index)
    if is_alive(A, index):
        if n_alive_neighbors_ not in (2, 3):
            return 0  # Kill
        else:
            return 1  # Stayin alive
    else:  # Dead
        if n_alive_neighbors_ == 3:  # Otherwise it stays dead
            return 1  # Bring to life
        else:
            return 0  # The dead may never die

Now write a function that finds the next state of the entire world:

def evolve(A: np.array):
    """Returns a new array that is the next generation from A"""
    B = np.zeros(A.shape)  # Initialize
    for i, row in enumerate(A):
        for j, cell in enumerate(row):
            B[i, j] = evolve_cell(A, (i, j))
    return B

Finally we can write the game_of_life() function:

def game_of_life(A: np.ndarray, n_generations=5):
    """Plays a game of life on a torus with starting state A"""
    states = [A]  # Initialize list of states to start with A (gen 0)
    for igen in range(1, n_generations):
        states.append(evolve(states[igen-1]))
    return states

Define Starting States

Now to define the given starting states:

start1 = np.array([
    [0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0],
    [0, 1, 1, 1, 0],
    [0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0],
])
start2 = np.zeros((20, 20))
start2[9:12, 9:12] = np.array([
    [0, 1, 0],
    [0, 0, 1],
    [1, 1, 1],
])
start3 = np.zeros((40, 40))
gg = np.array(engcom.data.game_of_life_starts("gosper_glider"))
start3[15:(15+gg.shape[0]), 2:(2+gg.shape[1])] = gg
starting_states = (start1, start2, start3)

Functions for Visualizing Games

The following function can visualize a single state of the world:

def visualize_state(state: np.ndarray):
    """Visualize a single state of a game"""
    fig, ax = plt.subplots()
    ax.matshow(state)
    plt.show()

Visualize the Games with Animation

The following function will be used to update the animation:

def update(j):
    """Update plot with a new game state"""
    plot.set_data(game_states[j])
    return plot

This function will animate a game (a list of states) and save the animation to a gif:

def animate_game(game_states: list):
    global plot  # Needed global to communicate with update()
    fig = plt.figure(num=0)
    plot = plt.imshow(game_states[0], cmap='gray_r', vmin=0, vmax=1)
    ax = plt.gca()
    ax.set(xticks=[], yticks=[], xticklabels=[], yticklabels=[])
    ax.spines['top'].set_visible(True)
    ax.spines['right'].set_visible(True)
    ax.spines['bottom'].set_visible(True)
    ax.spines['left'].set_visible(True)
    anim = animation.FuncAnimation(
        func=update,
        fig=fig,
        frames=len(game_states), 
        interval = 250,
    )
    # plt.show()
    anim.save(f"animation_{i}.gif")

Now just to loop through the starting states, play the games, and animate:

for i, A in enumerate(starting_states):
    game_states = game_of_life(A, n_generations=81)
    animate_game(game_states)

Visualize the Games with a Grid

This function will arrange visualizations of a game (a list of states) in a grid:

def game_gridder(game_states: list):
    nrows = int(np.floor(np.sqrt(len(game_states))))
    ncols = int(np.ceil(len(game_states)/nrows))
    fig, axs = plt.subplots(num=1, nrows=nrows, ncols=ncols)
    fig.subplots_adjust(
        left=0, right=1, top=1, bottom=0, wspace=0, hspace=0
    )
    for i, state in enumerate(game_states):
        ax = axs.flat[i]
        ax.matshow(game_states[i], cmap='gray_r', vmin=0, vmax=1)
        ax.set(xticks=[], yticks=[], xticklabels=[], yticklabels=[])
        ax.spines['top'].set_visible(True)
        ax.spines['right'].set_visible(True)
        ax.spines['bottom'].set_visible(True)
        ax.spines['left'].set_visible(True)
    n_extra = nrows * ncols - len(game_states)
    for i in range(len(game_states), nrows*ncols):
        axs.flat[i].set_axis_off()
    return fig

Now just to loop through the starting states, play the games, and plotting on a grid:

for i, A in enumerate(starting_states):
    game_states = game_of_life(A, n_generations=81)
    game_gridder(game_states)
    plt.show()

In robotic path planning, it is often important to know if a given point (e.g., a potential location of the robot) is inside of a given polygon (e.g., a shape representing an obstacle). On a plane, a polygon can be defined by a list of $n$ points $(x_i, y_i)$ representing the vertices of the polygon $P$. One algorithm for determining if a given point $R$ is in $P$ is called the winding number algorithm, which computes the winding number $\omega$ as the sum of the angles $\theta_i$ between the vectors from $R$ to consecutive vertices $P_i$ and $P_{i+1}$ of the polygon, denoted $\vec{r}_i$ and $\vec{r}_{i+1}$, as shown in fig. ¿fig:r9-polygon?. In other words, \[ \omega = \frac{1}{2\pi} \sum_{i=0}^{n-1} \theta_i. \qquad{(1)}\]

Figure 3.7: A polygon and vectors from R to two consecutive vertices.

It can be shown that if the winding number is $0$, then $R$ is outside the polygon; otherwise, it is inside. The angle $\phi_i$ of vector $\vec{r}_i = [r_{i_x}, r_{i_y}]$ is \[ \phi_i = \arctan(r_{i_y}/r_{i_x}), \] where we should use np.atan2(riy, rix) for computation. The difference between the angles of two consecutive vectors is \[ \theta_i = \phi_{i+1} - \phi_i \text{ where } |\theta_i| \le \pi. \] The bound $|\theta_i| \le \pi$ must be enforced because the acute angle is used in eq. 1, so if $\phi_{i+1} - \phi_i < -\pi$, we should add $2\pi$ and if $\phi_{i+1} - \phi_i > \pi$, we should subtract $2\pi$.

Write a program that meets the following requirements:

It defines a Polygon class that is constructed with instance attribute vertices, a list of $(x_i, y_i)$ coordinate tuples defining the vertices of the polygon.
The Polygon class has a method plot() that plots the polygon as a closed curve. If a point R is passed to the plot() method, it should appear as a single point on the plot.
The Polygon class has a method is_inside(R) that checks if the point R (a tuple) is inside the polygon using a winding number algorithm. The method should return True if R is inside the polygon and False otherwise. Additional methods can be added to help with the computation of angles and other intermediate quantities.
It tests the Polygon class with the following polygons and points, testing if the points are inside the polygon and plotting the polygon with the points:
1. $P = [(5, 1), (2, 3), (-2, 3.5), (-4, 1), (-2, 1.5), (-2, -2), (-5, -3), (2, -2.5),\allowbreak (5.5, -1)]$ and the points $R_1 = (0, 0)$ and $R_2 = (-4, 0)$
2. $P = [(4, 1), (1, 2), (-1, 1), (-4, 2), (-5, -2), (-3, -2), (-5, -3), (2, -2),\linebreak (5, -2)]$ and the points $R_1 = (0, 0)$ and $R_2 = (-4, 0)$

Restriction: Use only the NumPy and Matplotlib packages.

Load packages:

import numpy as np
import matplotlib.pyplot as plt

Define the Class

Here is a basic Polygon class definition:

class Polygon:
    """A polygon defined by its vertices."""
    def __init__(self, vertices: list):
        self.vertices = vertices  # Vertices attribute

This can be initialized with a list of vertex tuples. Now let’s add a plot() method as follows:

class Polygon:
    """A polygon defined by its vertices."""
    def __init__(self, vertices: list):
        self.vertices = vertices  # Vertices attribute
        self.n = len(vertices)  # Number of vertices
    
    def plot(self, R: tuple | None = None, ax=None):
        """Plot the polygon as a closed curve and a point R"""
        vertices_a = np.array(self.vertices).T
        # Initialize and fill most of the x, y arrays
        x = np.full(self.n + 1, np.nan)  # +1 to repeat the first point
        x[0:self.n] = vertices_a[0]
        y = np.full(self.n + 1, np.nan)  # +1 to repeat the first point
        y[0:self.n] = vertices_a[1]
        # Add the first point to the end
        x[-1] = x[0]
        y[-1] = y[0]
        # Plot
        if ax is None:
            fig, ax = plt.subplots()
        ax.plot(x, y)
        if R is not None:
            ax.scatter(R[0], R[1])
            ax.annotate(
                "$R$", R, 
                xytext=(5, 5), textcoords='offset points'
            )
        return ax
        
    def vectors(self, R:tuple) -> list:
        """Returns the vectors from R to vertices of the polygon"""
        rs = []
        for p in self.vertices:
            rs.append(np.array(p) - np.array(R))  # Difference in x, y
        return rs
    
    def winding_number(self, R: tuple) -> int:
        """Returns the winding number for point R"""
        rs = self.vectors(R)  # Vectors from R to each P vertex
        thetas = np.full(self.n, np.nan)  # Initialize
        for i in range(self.n):
            j = (i+1) % self.n  # Wrap i+1 back to 0
            phi_i = np.arctan2(rs[i][1], rs[i][0])  # Angle of $\vec{r}_i$
            phi_j = np.arctan2(rs[j][1], rs[j][0])  # Angle of $\vec{r}_{i+1}$
            if phi_j - phi_i < -np.pi:  # Handle 2*pi wrapping
                thetas[i] = phi_j - phi_i + 2*np.pi  # $\theta_i$
            elif phi_j - phi_i > np.pi:  # Handle 2*pi wrapping
                thetas[i] = phi_j - phi_i - 2*np.pi  # $\theta_i$
            else:
                thetas[i] = phi_j - phi_i  # $\theta_i$
        omega = np.round(1/(2*np.pi) * np.sum(thetas))
        return omega
    
    def is_inside(self, R: tuple) -> bool:
        """Returns True if R in polygon, False otherwise"""
        omega = self.winding_number(R)
        if omega == 0:
            return False
        else:
            return True

Test the Class

Create the test polygons and points as follows:

polygons = [
    Polygon([
        (5, 1), (2, 3), (-2, 3.5), (-4, 1), (-2, 1.5), 
        (-2, -2), (-5, -3), (2, -2.5), (5.5, -1)
    ]),
    Polygon([
        (4, 1), (1, 2), (-1, 1), (-4, 2), (-5, -2), 
        (-3, -2), (-5, -3), (2, -2), (5, -2)
    ]),
]
points = [(0, 0), (-4, 0)]

Graph the polygons and points as follows:

places = []
fig, axs = plt.subplots(2, 2)
for i, poly in enumerate(polygons):
    for j, R in enumerate(points):
        poly.plot(R, ax=axs[i, j])
        if poly.is_inside(R):
            places.append(f"R = {R} is inside polygon {i}")
        else:
            places.append(f"R = {R} is outside polygon {i}")
plt.show()

Print the test results from is_inside() as follows:

for place in places:
    print(place)

R = (0, 0) is inside polygon 0
R = (-4, 0) is outside polygon 0
R = (0, 0) is inside polygon 1
R = (-4, 0) is inside polygon 1

Create and test the Polygon class described in problem 3.11 with the following augmentation: instead of using the tangent function to compute the angle $\theta_i$ between consecutive vectors, use vector products.

The magnitude of the angle can be found from the dot product (np.dot()) expression \[ \cos\theta_i = \frac{\vec{r}_i \cdot \vec{r}_{i+1}}{|\vec{r}_i|\,|\vec{r}_{i+1}|}, \] where $|\cdot|$ is the vector length and $\cdot$ is the dot product. To determine the sign of the angle, the cross product (np.cross()) can be used because $\vec{r}_i \times \vec{r}_{i+1} = -\vec{r}_{i+1} \times \vec{r}_{i}$. In other words, the direction of rotation between $\vec{r}_i$ and $\vec{r}_{i+1}$ can be found from the sign of their cross product.

In this problem, we will develop a model of a mass’s trajectory on a 2D rectangular region of space with walls off which the mass bounces elastically (i.e., with no loss of kinetic energy). We divide time into discrete moments $t_0, t_1, \cdots, t_{n-1}$ with equal intervals $\Delta t$. The mass’s position at time $t_i$ is given by the vector $\vec{r}_i = \begin{bmatrix} x_i & y_i \end{bmatrix}^\top$. In discrete time, the mass’s velocity at time $t_i$ is given by the vector \[ \vec{v}_i = \frac{\vec{r}_{i} - \vec{r}_{i-1}}{\Delta t}, \] which approximates the derivative of the position vector with respect to time.

We will model the trajectories of several masses. In this problem, there will be no forces acting on the masses, so they will move at constant speed. (See problem 3.14 for an extension with force fields.)

Write a program that meets the following requirements:

Define a class State that represents the state of a mass at a given time. The class should have the following attributes:
1. position, an array representing the position at the current time
2. velocity, an array representing the velocity at the current time
Define a class Region that represents the 2D rectangle in which the mass moves. Use Cartesian coordinates with the origin at the lower-left corner of the rectangle, the $x$-axis increasing to the right, and the $y$-axis increasing upwards. The class should have the following attributes:
1. width, the width of the region
2. height, the height of the region
Define a class Trajectory that represents the trajectory of a mass through time. The class should have the following attributes:
1. initial_state, a State instance representing the mass’s initial state (required to initialize)
2. time, an array representing the corresponding times (required to initialize)
3. region, a Region instance representing the region in which the mass moves (required to initialize)
4. states, a list of State instances representing the mass’s states at each time step
Define an instance method bounce(s: State) -> State for the Trajectory class that updates the mass’s position and velocity when it bounces off a wall. (The angle of incidence equals the angle of reflection, but no trigonometry is needed.)
Define an instance method step_state() for the Trajectory class that advances the mass’s state by one time step. Wall bouncing should be handled here by calling bounce().
Define an instance method compute() for the Trajectory class that advances the mass’s state through time until a specified end time. This method should call step_state() repeatedly and populate the states attribute.
Define class instance methods get_x() and get_y() for the Trajectory class that return arrays of the $x$ and $y$ components of the mass’s positions at each time step.
Define a class Trajectories that represents a collection of Trajectory instances. The class should have the attribute trajectories, a list of Trajectory instances. It should be initialized with a list of initial states, a time array, and a region, from which it constructs the Trajectory instances.
Define an instance method plot() for the Trajectories class that plots the trajectories of all masses in the collection. The method should plot the trajectories as differently colored lines and the initial positions as points.
Test the classes and methods with the following parameters:
1. A time interval of $[0, 10]$ s with $101$ time steps
2. A region with width $10$ and height $5$. (Use units of meters.)
3. A collection of $3$ trajectories with the following initial states:
  1. Mass $0$ at $(1, 1)$ m with velocity $(1, 0)$ m/s
  2. Mass $1$ at $(2, 2)$ m with velocity $(0, 1)$ m/s
  3. Mass $2$ at $(3, 3)$ m with velocity $(1, -1)$ m/s
  4. Mass $3$ at $(4, 4)$ m with velocity $(0.2, 0.6)$ m/s

In this problem, we will extend the model of a mass’s trajectory on a 2D rectangular region of space with walls from problem 3.13 to include force fields acting on the masses. A force field is a vector field $\vec{F}(\vec{r})$ that gives the force acting on the mass at position $\vec{r}$. If $\vec{F}$ is the only force, the acceleration of the mass at position $\vec{r}$ is given by Newton’s second law: \[ \vec{a}(\vec{r}) = \vec{F}(\vec{r})/m, \] where $m$ is the mass of the object. In discrete time, the mass’s velocity at time $t_i$ is given by the vector \[ \vec{v}_{i} = \vec{v}_{i-1} + \vec{a}_i \Delta t, \] where $\vec{a}_i = \vec{F}(\vec{r}_i)/m$.

A force field can be represented by a function that takes a position vector as input and returns a force vector. Consider the following force fields: $$\begin{align} \vec{F}_0(\vec{r}) &= c_0 \frac{\vec{r} - \vec{r}_0}{|\vec{r} - \vec{r}_0|^2}\text{ and} \\ \vec{F}_1(\vec{r}) &= c_1 \frac{\vec{r}_1 - \vec{r}}{|\vec{r}_1 - \vec{r}|^2}, \end{align}$$ where $\vec{r}_0, \vec{r}_1$ are fixed positions vectors, $|\cdot|$ is the vector length, and $c_0, c_1$ are positive constants with units N$\cdot$m that determine the strength of the force fields. $\vec{F}_0$ points radially outward from $\vec{r}_0$ so it is a “repeller,” and $\vec{F}_1$ points radially inward toward $\vec{r}_1$ so it is an “attractor.”

Write a program that meets the following requirements, starting with the classes from problem 3.13:

Extend the Trajectory class to include a mass attribute that represents the mass of the object.
Extend the Trajectory class to include a list of force fields as an attribute. The force fields should be functions that take a position vector as input and return a force vector.
Extend the step_state() method of the Trajectory class to include a damping force given by \[ \vec{F}_\text{damp} = -0.3 |\vec{v}| \vec{v}. \]
Extend the step_state() method of the Trajectory class to include the acceleration due to the sum of the force fields acting on the mass and the damping force.
Extend the Trajectories class constructor to accept a list of mass values and a list of force fields for each trajectory.
Define the force fields $\vec{F}_0$ and $\vec{F}_1$ as functions and test the program with the following parameters:
1. Force field vectors $\vec{r}_0 = (10, 6)$ and $\vec{r}_1 = (10, 4)$ and strength constants $c_0 = c_1 = 0.1$ N$\cdot$m
2. A time interval of $[0, 700]$ s with $701$ time steps
3. A region with width $20$ and height $10$. (Use units of meters.)
4. A collection of trajectories with the following masses and initial states:
  1. Mass $m_0 = 1$ kg at $(3, 4.98)$ m and at rest
  2. Mass $m_1 = 1$ kg at $(3, 4.99)$ m and at rest
  3. Mass $m_2 = 1$ kg at $(3, 5.01)$ m and at rest
  4. Mass $m_3 = 1$ kg at $(3, 5.02)$ m and at rest

Johnston, Nathaniel, and Dave Greene. 2022. Conway’s Game of Life: Mathematics and Construction. Self-published. https://doi.org/10.5281/zenodo.6097284.

Online Resources for Section 3.5

No online resources.