Tutorial A11 – Solutions¶
1 Write a function that takes three vectors \(\mathbf{a}\), \(\mathbf{b}\), \(\mathbf{c}\) in 2D (NumPy arrays) and checks whether they form a triangle, i.e. if \(a + b = \pm c\).
[1]:
import numpy as np
[2]:
def is_triangle(a, b, c):
"""Check if a, b, c form a triangle
Args:
a, b, c (numpy.ndarray): Input vectors
Returns:
True or False
"""
if np.allclose(a + b, np.absolute(c)):
return True
return False
[3]:
is_triangle(
np.array([0, 1]), # a
np.array([1, 0]), # b
np.array([1, 1]) # c
)
[3]:
True
2 Consider the vector \(\mathbf{v} = \begin{pmatrix}0& 0& 1\end{pmatrix}^\mathrm{T}\) and the rotation matrix
\begin{equation*} \mathbf{R} = \begin{pmatrix} 1 & 0 & 0\\ 0 & 0 & -1\\ 0 & 1 & 0 \end{pmatrix} \end{equation*}
Rotate \(\mathbf{v}\) by \(\mathbf{R}\) to obtain \(\mathbf{v}^\prime\)
[4]:
v = np.array([0, 0, 1])
R = np.array([[1, 0, 0],
[0, 0, -1],
[0, 1, 0]])
v_ = np.dot(R, v)
print(v_)
[ 0 -1 0]
Compute the scalar product of \(\mathbf{v}\) and \(\mathbf{v}^\prime\)
[5]:
print(np.dot(v, v_))
0
By which angle did we rotate along which axis?
90° around x-axis
\begin{equation*} \mathbf{R} = \begin{pmatrix} 1 & 0 & 0\\ 0 & \cos \theta & -\sin \theta\\ 0 & \sin \theta & \cos \theta \end{pmatrix} \end{equation*}
3 Consider the second Laplace spherical harmonic
\[Y_1^{-1}(\theta, \phi) = \frac{1}{2}\sqrt{\frac{3}{2\pi}}\sin\theta\mathrm{e}^{-i\phi}\]
Compute \(Y\) on a grid of 10 points per dimension (\(0 \leq \theta \leq \pi\), \(0 \leq \phi < 2\pi\)).
[6]:
theta = np.linspace(0, np.pi, 10)
phi = np.linspace(0, 2*np.pi, 10)
print(r"$\theta$", theta)
print(r"$\phi$", phi)
$\theta$ [0. 0.34906585 0.6981317 1.04719755 1.3962634 1.74532925
2.0943951 2.44346095 2.7925268 3.14159265]
$\phi$ [0. 0.6981317 1.3962634 2.0943951 2.7925268 3.4906585
4.1887902 4.88692191 5.58505361 6.28318531]
[7]:
def Y(theta, phi):
return 1/2 * np.sqrt(3 / (2*np.pi)) * np.sin(theta) * np.exp(-1j * phi)
T, P = np.meshgrid(theta, phi)
y = Y(T, P)
print(y)
print(y.shape)
[[ 0.00000000e+00+0.00000000e+00j 1.18165959e-01+0.00000000e+00j
2.22079358e-01+0.00000000e+00j 2.99206710e-01+0.00000000e+00j
3.40245317e-01+0.00000000e+00j 3.40245317e-01+0.00000000e+00j
2.99206710e-01+0.00000000e+00j 2.22079358e-01+0.00000000e+00j
1.18165959e-01+0.00000000e+00j 4.23108304e-17+0.00000000e+00j]
[ 0.00000000e+00+0.00000000e+00j 9.05203759e-02-7.59556140e-02j
1.70122659e-01-1.42749860e-01j 2.29205638e-01-1.92326366e-01j
2.60643034e-01-2.18705474e-01j 2.60643034e-01-2.18705474e-01j
2.29205638e-01-1.92326366e-01j 1.70122659e-01-1.42749860e-01j
9.05203759e-02-7.59556140e-02j 3.24119765e-17-2.71968776e-17j]
[ 0.00000000e+00+0.00000000e+00j 2.05193034e-02-1.16370752e-01j
3.85636759e-02-2.18705474e-01j 5.19567000e-02-2.94661088e-01j
5.90829793e-02-3.35076226e-01j 5.90829793e-02-3.35076226e-01j
5.19567000e-02-2.94661088e-01j 3.85636759e-02-2.18705474e-01j
2.05193034e-02-1.16370752e-01j 7.34719860e-18-4.16680338e-17j]
[ 0.00000000e+00-0.00000000e+00j -5.90829793e-02-1.02334722e-01j
-1.11039679e-01-1.92326366e-01j -1.49603355e-01-2.59120612e-01j
-1.70122659e-01-2.94661088e-01j -1.70122659e-01-2.94661088e-01j
-1.49603355e-01-2.59120612e-01j -1.11039679e-01-1.92326366e-01j
-5.90829793e-02-1.02334722e-01j -2.11554152e-17-3.66422540e-17j]
[ 0.00000000e+00-0.00000000e+00j -1.11039679e-01-4.04151381e-02j
-2.08686334e-01-7.59556140e-02j -2.81162338e-01-1.02334722e-01j
-3.19726014e-01-1.16370752e-01j -3.19726014e-01-1.16370752e-01j
-2.81162338e-01-1.02334722e-01j -2.08686334e-01-7.59556140e-02j
-1.11039679e-01-4.04151381e-02j -3.97591751e-17-1.44711563e-17j]
[-0.00000000e+00+0.00000000e+00j -1.11039679e-01+4.04151381e-02j
-2.08686334e-01+7.59556140e-02j -2.81162338e-01+1.02334722e-01j
-3.19726014e-01+1.16370752e-01j -3.19726014e-01+1.16370752e-01j
-2.81162338e-01+1.02334722e-01j -2.08686334e-01+7.59556140e-02j
-1.11039679e-01+4.04151381e-02j -3.97591751e-17+1.44711563e-17j]
[-0.00000000e+00+0.00000000e+00j -5.90829793e-02+1.02334722e-01j
-1.11039679e-01+1.92326366e-01j -1.49603355e-01+2.59120612e-01j
-1.70122659e-01+2.94661088e-01j -1.70122659e-01+2.94661088e-01j
-1.49603355e-01+2.59120612e-01j -1.11039679e-01+1.92326366e-01j
-5.90829793e-02+1.02334722e-01j -2.11554152e-17+3.66422540e-17j]
[ 0.00000000e+00+0.00000000e+00j 2.05193034e-02+1.16370752e-01j
3.85636759e-02+2.18705474e-01j 5.19567000e-02+2.94661088e-01j
5.90829793e-02+3.35076226e-01j 5.90829793e-02+3.35076226e-01j
5.19567000e-02+2.94661088e-01j 3.85636759e-02+2.18705474e-01j
2.05193034e-02+1.16370752e-01j 7.34719860e-18+4.16680338e-17j]
[ 0.00000000e+00+0.00000000e+00j 9.05203759e-02+7.59556140e-02j
1.70122659e-01+1.42749860e-01j 2.29205638e-01+1.92326366e-01j
2.60643034e-01+2.18705474e-01j 2.60643034e-01+2.18705474e-01j
2.29205638e-01+1.92326366e-01j 1.70122659e-01+1.42749860e-01j
9.05203759e-02+7.59556140e-02j 3.24119765e-17+2.71968776e-17j]
[ 0.00000000e+00+0.00000000e+00j 1.18165959e-01+2.89423126e-17j
2.22079358e-01+5.43937551e-17j 2.99206710e-01+7.32845080e-17j
3.40245317e-01+8.33360677e-17j 3.40245317e-01+8.33360677e-17j
2.99206710e-01+7.32845080e-17j 2.22079358e-01+5.43937551e-17j
1.18165959e-01+2.89423126e-17j 4.23108304e-17+1.03631646e-32j]]
(10, 10)
4 Given 1000 randomly places points in 2D, compute all the pairwise distances according to the formula below.
\(d_{ij} = \sqrt{(x_i-x_j)^2+(y_i-y_j)^2}\) with \(i, j \in (1,2,...,1000)\)
Use three different Methods for this and compare the times:
a) A naive double for loop
How many distances are larger than 0.5?
[8]:
import time
import numpy as np
samples = 1000 # number of points
p = (
np.random.random(samples), # x
np.random.random(samples), # y
)
x = p[0]
y = p[1]
[9]:
t0 = time.time()
D = np.zeros((samples, samples))
# Array container to store the distances
for i in range(samples):
for j in range(samples):
if i >= j:
# We make use of the symmetry
# and only calculate whats necessary
continue
D[i, j] = D[j, i] = np.sqrt((x[i] - x[j])**2 + (y[i] - y[j])**2)
print(f"Double for loop needs {time.time() - t0:.6f} seconds")
print("Number of elements larger than 0.5: ", (D > 0.5).sum())
Double for loop needs 1.968029 seconds
Number of elements larger than 0.5: 523098
Advanced b) A vectorized version using meshgrids
[10]:
t0 = time.time()
X, Y = np.meshgrid(x, y)
# Get all combinations of x and y
D = np.sqrt((X.T - X)**2 + (Y.T - Y)**2)
print(f"Meshgrid needs {time.time() - t0:.6f} seconds")
print("Number of elements larger than 0.5: ", (D > 0.5).sum())
Meshgrid needs 0.025411 seconds
Number of elements larger than 0.5: 523098
Advanced c) Outer subtraction
[11]:
t0 = time.time()
D = np.sqrt(
np.subtract.outer(x, x)**2 + np.subtract.outer(y, y)**2
)
print(f"Substract.outer needs {time.time() - t0:.6f} seconds")
print("Number of elements larger than 0.5: ", (D > 0.5).sum())
Substract.outer needs 0.020332 seconds
Number of elements larger than 0.5: 523098
Advanced d): The cdist function from the scipy module (scipy.spatial.distance.cdist)
[12]:
from scipy.spatial.distance import cdist
P = np.array(p).T
# We need our points in a different format, though
[13]:
t0 = time.time()
D = cdist(P, P)
print(f"Scipy's cdist needs {time.time() - t0:.6f} seconds")
print("Number of elements larger than 0.5: ", (D > 0.5).sum())
Substract.outer needs 0.009225 seconds
Number of elements larger than 0.5: 523098