The Coronavirus Curve – A simple simulator

Here is an implementation for a simple ODE simulation of the COVID-19 curve.

The simulator runs based on a set of first order ordinary differential equations (ODEs). These equations are dependent and measure three values including the change in the number of susceptible cases (S), the change in the number of infected cases (I), and the change in the number of recovered cases (R). It has to be noted that the R set includes both dead or recovered alive patients.

The important factors in this model include two weight values: the transmission rate and the infection rate. In the following program these rates are considered as constants. However, they can change based on many factors such as several geographical differences.

For more details about the simulator equations watch the video (link below).

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

def simulator(y, t, b, c):
S,I,R = y
dydt = [-b*S*I, b*S*I - c*I, c*I]
return dydt

b = 2.2 #3.2 # Transmission rate
c = 0.1 #0.23 # Recovery rate

S_0 = 0.99 # initial percentage of susseptible cases
I_0 = 1- S_0 #0.01 #initial percentage of infected cases
R_0 = 0.0 # initial percentage of recovered

# initial values
y0 = [S_0, I_0, R_0]

tmax = 20 # max time
t = np.linspace(0, tmax, 101)
sol = odeint(simulator, y0, t, args=(b, c))

# plot
plt.plot(t, sol[:, 0], 'b', label='S(t) Susseptible')
plt.plot(t, sol[:, 1], 'r', label='I(t) Infected')
plt.plot(t, sol[:, 2], 'g', label='R(t) Recovered')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.title(['Transmision rate= '+str(b)+' Recovery rate= '+str(c)+' Ratio= '+str(b/c)])
plt.show()

This post is based on this Numberphile video

Star / Asterisks Operator ( * ) in Python – 1

1. Unpacking iterables into a list/tuple

In addition to serving as the symbol for multiplication, the star operator (*) in Python has several other applications. The first application is unpacking iterables into a list/tuple

Example: Adding (appending) elements to an array

>>> n = [2, ,3, 4, 5, 6, 1]
>>> n
[2, 3, 4, 5, 6, 1]
>>> m = [n, 2]
>>> m
[[2, 3, 4, 5, 6, 1], 2]

instead,

>>> n = [2,3,4,5,6,1]
>>> n
[2, 3, 4, 5, 6, 1]
>>> m = [*n, 2]
>>> m
[2, 3, 4, 5, 6, 1, 2]

Another example:

>>> chars = ['a', 'b', 'c', 'd']
>>> print(chars)
['a', 'b', 'c', 'd']
>>> print(chars[0], chars[1], chars[2], chars[3])
a b c d
>>>print(*chars)
a b c d

2. Packing function arguments

The star operator can be used for packing arguments sent to a function. In other words, we can develop functions that can accept any number of arguments. For instance,

>>> def collect_all(*n):
return sum(n)

>>> collect_all(1)
1
>>> collect_all(1,2,3)
6
>>>

The star or asterisk operator has other applications that we will discuss in another post.

3D Rotations using Rodrigues Rotation Formula

The problem we are addressing here is the rotation of a general 3D vector about a given axis of rotation denoted by by radians.

The Rodrigues Rotation Formula is as follows:

Example: for a sanity check we can consider that will result in

Example:

Rotate about by .

import numpy as np
v = np.array([2.,0.,0.])
n = np.array([0.,0.,1.])
theta = np.pi/2
vp = (1-np.cos(theta))*(np.dot(v,n))*n + np.cos(theta)*v + np.sin(theta)*(np.cross(n,v))
>>> vp
array([1.2246468e-16, 2.0000000e+00, 0.0000000e+00])

 

 

 

Extension of 2D Complex Exponential Formula to 3D Rotation

In our previous post, we used , where and are 2D vectors and the formula rotates vector by .

This idea can be extended to 3D by converting the vectors to quaternions. Assume and . Also, we consider to be a unit vector normal to the plane where and   live.

The formula can be written as

where

Example:

Rotate about axis by .

Answer:

First write the vectors as quaternions.

Then

which can be easily verified.

Python code:

from pyquaternion import Quaternion as Q
v = Q((0,0,0,1))
n = Q((0,0,1,0))
vp = Q.exp(np.pi/2 * n) * v
>>> Quaternion(0.0, 1.0, 0.0, 6.123233995736766e-17)

2D Vector Rotation using Complex Exponential

Complex exponential formula represents the relationship between a vector and the rotated version of that vector by as

where .

Example:

Rotate vector by .

Answer:

which is a vector towards Y-axis (Im axis).

In python, we can compute this as follows:

v = 1. + 0.j
vp = (np.exp((np.pi/2)*1.j))*v
print(vp)
(6.123233995736766e-17+1j)

Note: in python, we use `j` to represent the unit imaginary number.

What is Kronecker Product?

Kronecker Product is a generalization of the outer product of two arbitrary size matrices and results in a block matrix.

and then

Example:

Python example using numpy

import numpy as np
A = np.array([[1.,2.],[3.,4.]])
B = np.array([[0.,-1.],[-1.,2.]])
C = np.kron(A,B)

[[ 0. -1. 0. -2.]
[-1. 2. -2. 4.]
[ 0. -3. 0. -4.]
[-3. 6. -4. 8.]]

 

 

Quaternions – Use them for Rotation and Transformation

This is the second tutorial on Quaternions. In the previous post, I explained how to use basics of quaternions from the pyquaternion python library. Check the previous post here.

In this post, I will show you how to perform spatial rotation using quaternions.

First, import the modules we need for this tutorial:

import numpy as np
import matplotlib.pyplot as plt
from pyquaternion import Quaternion
# This import registers the 3D projection, but is otherwise unused.
from mpl_toolkits.mplot3d import Axes3D

Let’s start by creating a circle in 3D coordinate:

s = np.linspace(0,2*np.pi,100)
x = np.cos(s)
y = np.sin(s)
z = np.zeros_like(x)
ax.plot(x,y,zs=0,zdir='z',label='transformation with Quaternions')
plt.show()

This code will result in a circle in X-Y plane plotted in XYZ space.

Now let’s define a quaternion for rotation. The best way is to use the ‘axis’ – ‘angle’ method as follows:

q = Quaternion(axis=[0, 1, 0], angle=np.pi / 2)

This quaternion will rotate a point by 90 degrees about Y-axis

Now let’s apply it to the circle. To do this, first concatenate all the points into a 3d-array and then use a for-loop to rotate each point as follows:

P = np.stack((x,y,z))
# make an empty matrix
Pr = np.zeros((len(s),3))
# rotate points one by one
for i in range(len(s)):
Pr[i][:] = q2.rotate([x[i], y[i], z[i]])

# transpose for shape consistency
Pr = Pr.transpose()

Now, we only need to plot the rotated matrix as follows:

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(x,y,zs=0,zdir='z',label='transformation with Quaternions')
ax.plot(Pr[:][0], Pr[:][1], Pr[:][2],'r')
plt.show()

An the result is the red circle in the following figure which is the blue circle rotated 90 degrees about the Y-axis

Quaternions – How to generate and operate them in Python?

Python Package

There are several python libraries (modules) that you can install. Here, I am going to use pyquaternion

Installation

This package can be installed as follows

pip install pyquaternion

Usage

After installing the package, use it as follows:

from pyquaternion import Quaternion

 Forming Quaternions

This package offers several methods for creating quaternions, some of which are very useful. I list the most useful ones here:

1. From `axis`   and `angle`:

q1 = Quaternion(axis=[0., 0., 1.], angle=3.14)

2. From list/tuple/array

q2 = Quaternion([1., 0., 0., 0.])
q3 = Quaternion((1., 0., 0., 0.))
q4 = Quaternion(np.array([1.0, 0., 0., 0.]))

3. From elements Q=[w,x,y,z]

q5 = Quaternion(w=1.0, x=0., y=0., z=0.)

 

4. using `scalar` – `vector` method

q6 = Quaternion(scalar=1.0, vector=(0., 0., 0.))
q7 = Quaternion(scalar=1.0, vector=[0., 0., 0.])
q8 = Quaternion(scalar=1.0, vector=np.array([0., 0., 0.]))

 

Conversion

1. We can convert a rotation matrix into a quaternion as follows:

R = np.eye(3)
q9 = Quaternion(matrix=R)

 

2. We can also convert a transformation matrix into a quaternion as follows:

T = np.eye(4)
q10 = Quaternion(matrix=T)

Calculate Norm/Magnitude

q11 = Quaternion(np.array([1., 0., 0., 0.]))
print(q11.norm)
1.0
print(q11.magnitude)
1.0

Calculate Inverse

print(q11.inverse)
1.000 -0.000i -0.000j -0.000k

Calcualte Conjugate

print(q11.inverse)
1.000 -0.000i -0.000j -0.000k

Calculate Normalized

q12 = Quaternion(np.array([2.0, 1., 1., 1.]))
print(q12.normalised)
0.756 +0.378i +0.378j +0.378k

In the next post, we will discuss more advanced topics in playing with Quaternions.

In the next post, I discuss how to use quaternions for spatial rotation and transformation.

Run local jupyter notebook on a Pytorch container

Introduction

When, I was learning Pytorch I had already two separate Python version (2.x and 3.x) together with a lot of packages installed on my machine and I did not want to risk messing up my laptop. So, I decided instead to use a docker image that provides me with the latest installation of Pytorch on a linux-based system and isolates it from everything else on my machine.

I didn’t know much about docker and wanted to use a bunch of interactive python notebooks on my machine to use pytorch installed on the docker image.

My plan was to use an existing Pytorch image that comes with a bunch of other python packages (as an image) and run my local python files on that image without installing Pytorch or any other dependencies on my machine.

This can be done using Docker. The existing Pytorch docker image can be accessed for free here. And then I just needed to figure out a way to connect my local files to the docker image and tell my files to use python packages from that image.

Here is how:

0- you need to have docker installed (tutorials here)

 

1- Open a terminal (powershell in windows is a good choice)

In this tutorial, I use the docker for windows and hence the terminal I use is the windows powershell. Similar instructions can be followed for linux and MacOS.

 

2- pull the docker image for pytorch

docker pull pytorch/pytorch

 

3- run the docker image

docker run -it --name pytorch1 -v current_dir_including_ipynb_files:/workspace -p 5000:8888 -p 5001:6006 pytorch/pytorch

• –name gives a name to the container rather than a random one
• -v maps the folder in the local machine to a folder inside the container (mounts the first path to the second – divided by the : symbol )
• -p defines the ports on the local machine and the container

 

4- if the docker runs correctly, you should be inside the workspace folder of the container

root@ca09f087f880:/workspace#

 

5- if the mounting is done correctly you should see all the files in the local folder inside the workspace folder

root@ca09f087f880:/workspace# ls

file1.ipynb    file2.ipynb

 

7- check the installed packages in the container

pip freeze

 

6- the pytorch does not come with the jupyter package and you need to install it

pip instal jupyter

 

7- Now you want to start the jupyter notebook inside the container but it needs to talk to the local machine

jupyter notebook --ip 0.0.0.0 --port 8888 --allow-root &

 

8- This gives you a link that you can put into your browser to access the notebook

root@ca09f087f880:/workspace# [I 19:03:50.503 NotebookApp] Writing notebook server cookie secret to /root/.local/share/jupyter/runtime/notebook_cookie_secret
[I 19:03:50.813 NotebookApp] Serving notebooks from local directory: /workspace
[I 19:03:50.814 NotebookApp] The Jupyter Notebook is running at:
[I 19:03:50.814 NotebookApp] http://(ca09f087f880 or 127.0.0.1):8888/?token=643e3ac56e3de36398e7e9c27bc4e61b9b8158aab0d4c1f9
[I 19:03:50.814 NotebookApp] Use Control-C to stop this server and shut down all kernels (twice to skip confirmation).
[W 19:03:50.819 NotebookApp] No web browser found: could not locate runable browser.
[C 19:03:50.819 NotebookApp]

To access the notebook, open this file in a browser:
file:///root/.local/share/jupyter/runtime/nbserver-43-open.html
Or copy and paste one of these URLs:
http://(ca09f087f880 or 127.0.0.1):8888/?token=643e3ac56e3de36398e7e9c27bc4e61b9b8158aab0d4c1f9

 

9- you need to replace the port before using the url

change 127.0.0.1:8888 into 127.0.0.1:5000 and everything should work correctly.

10- Copy and paste the modified url into your browser and you should see something like this:

 

Stop the docker

1- exit from the container by pressing ctr+p, ctrl+q in the powershell

 

2- get the container ID using

docker ps

or

docker ps -l

 

CONTAINER ID IMAGE COMMAND CREATED STATUS PORTS NAMES
162ccf9ac9dc pytorch/pytorch "/bin/bash" 31 minutes ago Up 31 minutes 0.0.0.0:5001->6006/tcp, 0.0.0.0:5000->8888/tcp pytorch1

 

3- stop the container using

docker stop 162ccf9ac9dc

 

4- remove the docker using

docker rm 162ccf9ac9dc

 

 

Important Note:

Every time you run the docker, you need to install the additional dependencies (in this case, Jupyter) because the container is actually an image and cannot save changes.

 

Linear Transformation of Random Variables

Linear Transformation of Random Variables is a very basic yet important operation in Statistics results of which are extensively used in different methods and algorithms in Machine Learning.

Assume is a linear function,

,

where, is a random variable, and are a matrix and a vector respectively that form the linear transformation. If we calculate the mean and covariance for the transformation, we get:

1- for Mean

where

2- for Covariance

Both of these results are used in Machine Learning extensively.