<jupyter_start><jupyter_text>Own implementation of the EDMD algorithm (task 2) and application to two examples (tasks 3 and 4) The Extended Dynamic Mode Decomposition (EDMD) is a method to approximate the Koopman operator $\mathbf{K}$ and the correpsonding triples of Koopman eigenvalue, eigenfunction and mode. Where $\{(\mu_k, \varphi_k, v_k)\}^K_{k=1}$ is the set of K triples required to reconstruct the full state. (See Williams et al., 2015)<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import scipy.linalg<jupyter_output><empty_output><jupyter_text>Task 2: Own implementation of the EDMD algorithm Requirements:\Data set with 'labeled data', i.e. $F(x_i) = y_i$\Dictionary of observables, i.e. basisfunctions (best choice depends on data set)<jupyter_code># Calculate Koopman operator K
def calc_K(X, Y, D):
    M = np.shape(X)[1]
    
    # evaluate D with values from X
    phi_x = []
    for i in range(M):
        phi_x.append(D(X[0,i],X[1,i]))
    phi_x = np.array(phi_x)
    
    # evaluate D with values from Y
    phi_y = []
    for i in range(M):
        phi_y.append(D(Y[0,i],Y[1,i]))
    phi_y = np.array(phi_y)
    
    # calulate G and A matrices
    G = 1/M * phi_x.T @ phi_x
    A = 1/M * phi_x.T @ phi_y
    
    # solve for K
    K = np.linalg.pinv(G) @ A
    return K<jupyter_output><empty_output><jupyter_text>Eigenfunction of $\mathbf{K}$: $\phi_j = \Psi \zeta_j$ with $\zeta_j$ being the j-th eigenvector of $\mathbf{K}$<jupyter_code># caluclate eigenvalues, -vectors, -functions
def calc_eig(K, D):
    eigenvalues, left_evecs, right_evecs = scipy.linalg.eig(K, left=True, right=True)
    
    #sort eigenvalues (and associated eigenvectors) from big to small
    idx = eigenvalues.argsort()[::-1]
    eigenvalues = eigenvalues[idx]
    right_evecs = right_evecs[:, idx]
    left_evecs = left_evecs[:,idx]#left eigenvectors are returned for mode calculation
    
    #calculate eigenfunctions with right eigenvectors
    eigenfunctions = lambda x, y : D(x,y) @ right_evecs#the lambda is very useful here
    return eigenvalues, left_evecs, eigenfunctions<jupyter_output><empty_output><jupyter_text>Koopman modes: $\mathbf{V} = (\mathbf{W}^* \mathbf{B})^T$ with $\mathbf{W}^*$ holding the left eigenvectors. As mentioned in the presentation and the report, we are unsure how to obtain a suitable matrix B.<jupyter_code># calculate the Koopman modes
def calc_modes(eigenvectors, B):
    modes = (eigenvectors @ B).T    
    return modes
# perform EDMD algorithm
def EDMD(X, Y, D, B=None):
    K = calc_K(X, Y, D)
    eigenvalues, left_eigenvectors, eigenfunctions = calc_eig(K, D)
    if B is None:#B should be (K,N) where N is dim of data points
        B = np.ones((np.shape(left_eigenvectors)[-1], np.shape(X)[1]))
    modes = calc_modes(left_eigenvectors, B)
    return eigenvalues, eigenfunctions, modes#, K#, left_eigenvectors # if one is interested in them e.g. check shape<jupyter_output><empty_output><jupyter_text>Task 3: Linear example from Williams et al., 2015Analyical solution vs. EDMD approximation with our code\Example: $x(n+1) = Jx(n) = \begin{bmatrix} 0.9 & -0.1 \\ 0.0 & 0.8 \end{bmatrix} x(n)$ Analytical eigenvalues and eigenfunctions:<jupyter_code># analytical solution for eigenfunctions and -values -> formulae from Williams et al., 2015
def sol_func(x, y, i, j):
    return ((x - y)/(2**(1/2)))**i * y**j

def sol_evals(i, j):
    return (0.9)**i * (0.8)**j
# plot first 8 nontrivial eigenfunctions (+ colorbar) like in Williams et al., 2015 - analytical Koopman eigenfunctions
x = np.linspace(-5,5,100)
y = np.linspace(-5,5,100)
X, Y = np.meshgrid(x, y)
fig, axs = plt.subplots(2,4, figsize=(12,5), constrained_layout=True)
cmap = mpl.cm.jet#colormap
ij = [(1,0),(2,0),(0,1),(3,0),(1,1),(4,0),(2,1),(0,2)]#indices for the first 8 nontrivial eigenvalues, -functions
for i in range(1,9):
    values = sol_func(X,Y, *ij[i-1])
    ax = axs.flat[i-1]
    pcm = ax.pcolormesh(np.real(values), cmap=cmap)
    ax.set_xticks([0,25,50,75,100])
    ax.set_xticklabels(np.linspace(-5,5,5))
    ax.set_yticks([0,25,50,75,100])
    ax.set_yticklabels(np.linspace(-5,5,5))
    ax.set_title('$\mu_{%d}=(0.9)^{%d}(0.8)^{%d}=%.3f$' %(i,*ij[i-1],sol_evals(*ij[i-1])))
    ax.set_xlabel('x')
    ax.set_ylabel('y')


norm = mpl.colors.Normalize(vmin=-1, vmax=1)
plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), ax=axs[0,:])
plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), ax=axs[1,:])
plt.show()<jupyter_output><empty_output><jupyter_text>Eigenvalues and eigenfunctions with EDMD:<jupyter_code># Create the data
J = np.array([[0.9, -0.1],[0.0, 0.8]])
# X normal distributed with size (2,100)
X_data = np.random.default_rng().normal(size=(2,100))
# Y as result of J @ X
Y_data = J @ X_data
# no seed for the X data, because it should work for any value distribution and it's more fun
# Create the dictionary like in Williams et al. 2015
# define Hermite polynomials
def H0(x):
    return 1
def H1(x):
    return 2 * x
def H2(x):
    return 4 * x**2 - 2
def H3(x):
    return 8 * x**3 - 12 * x
def H4(x):
    return 16 * x**4 - 48 * x**2 + 12

# define dictionary
D = lambda x,y : [H0(x) * H0(y),H0(x)*H1(y),H0(x) * H2(y),H0(x)*H3(y),H0(x) * H4(y),
                  H1(x) * H0(y),H1(x)*H1(y),H1(x) * H2(y),H1(x)*H3(y),H1(x) * H4(y),
                  H2(x) * H0(y),H2(x)*H1(y),H2(x) * H2(y),H2(x)*H3(y),H2(x) * H4(y),
                  H3(x) * H0(y),H3(x)*H1(y),H3(x) * H2(y),H3(x)*H3(y),H3(x) * H4(y),
                  H4(x) * H0(y),H4(x)*H1(y),H4(x) * H2(y),H4(x)*H3(y),H4(x) * H4(y)]
# in Williams et al, 2015 a concrete matrix B is given, but without thorough explanation
B = np.zeros((25, 2))
B[1,0] = 1
B[5,1] = 1

# apply EDMD to obtain eigenvalues, eigenfunctions, and modes
eigenvalues, eigenfunctions, modes = EDMD(X_data, Y_data, D, B)
# plot first 8 nontrivial eigenfunctions (+ scaling) like in Williams et al. 2015 - EDMD Koopman eigenfunctions
x = np.linspace(-5,5,100)
y = np.linspace(-5,5,100)
X, Y = np.meshgrid(x, y)
fig, axs = plt.subplots(2,4, figsize=(12,5), constrained_layout=True)
cmap = mpl.cm.jet#colormap
for i in range(1,9):
    values = eigenfunctions(X,Y)[i]
    ax = axs.flat[i-1]
    pcm = ax.pcolormesh(np.real(values), cmap=cmap)
    ax.set_xticks([0,25,50,75,100])
    ax.set_xticklabels(np.linspace(-5,5,5))
    ax.set_yticks([0,25,50,75,100])
    ax.set_yticklabels(np.linspace(-5,5,5))
    ax.set_title('$\mu_{%d}$=%.3f' %(i,np.real(eigenvalues[i])))#here im-part is 0
    ax.set_xlabel('x')
    ax.set_ylabel('y')


norm = mpl.colors.Normalize(vmin=-1, vmax=1)
plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), ax=axs[0,:])
plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), ax=axs[1,:])
plt.show()<jupyter_output><empty_output><jupyter_text>Next, compare the full state observable (i.e. the reconstruction) with the original values. Further below is a comparison of the results based on analytically and numerically (EDMD) obtained Koopman eigenvalues, eigenfunctions, and modes.<jupyter_code># Reconstruction of X and Y (and Z) using the Koopman decomposition
t = 1 # one time step -> Y
x_rec = modes @ eigenfunctions(X_data[0,:],X_data[1,:]).T
y_rec = eigenvalues**t * modes @ eigenfunctions(X_data[0,:],X_data[1,:]).T
t = 2 # two time steps -> Z
z_rec = eigenvalues**t * modes @ eigenfunctions(X_data[0,:],X_data[1,:]).T
X_rec = np.array((x_rec[0], x_rec[1]))
Y_rec = np.array((y_rec[0], y_rec[1]))
Z_rec = np.array((z_rec[0], z_rec[1]))
Z_data = J @ Y_data # Z data generated with evolution operator

# (de-)active the data sets of your choice. That's easier than creating thousand subplots with all combinations
plt.scatter(X_rec[0,:], X_rec[1,:], label='X reconstructed')
plt.scatter(Y_rec[0,:], Y_rec[1,:], label='Y reconstructed')
plt.scatter(X_data[0,:], X_data[1,:], label='X original')
plt.scatter(Y_data[0,:], Y_data[1,:], label='Y original')
#plt.scatter(Z_rec[0,:], Z_rec[1,:], label='Z predicted')
#plt.scatter(Z_data[0,:], Z_data[1,:], label='Z original')
plt.legend()
plt.show()
# comparison to analytical reconstruction with eigenvalues, -functions from above and modes from Williams et al.
# according to Williams et al. only the 1st and 3rd nontrivial triples are important
my_evals = np.array((sol_evals(*ij[0]), sol_evals(*ij[2])))
my_modes = np.array([[0, -2**(1/2)],[-1, -1]])
my_funcs = lambda x,y : [sol_func(x,y, *ij[0]), sol_func(x,y, *ij[2])]
Y_anl = my_evals * my_modes @ my_funcs(X_data[0,:], X_data[1,:])

# no plotting of X for better visability
plt.scatter(Y_rec[0,:], Y_rec[1,:], label='Y reconstructed')
plt.scatter(Y_anl[0,:], Y_anl[1,:], label='Y analytic')
plt.scatter(Y_data[0,:], Y_data[1,:], label='Y original')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Task 4: Example from crowd dynamicsUsing the example from exercise 5 task 5: 'MI_timesteps.txt' containing the number of pedestrians in varoius measurement areas at campus Garching. The EDMD algorithm can be applied to a single, multiple, or all measurement areas as demonstrated below. After the data preprocessing, the EDMD algorithm is applied to the single measurement area 1 and to the all areas at once. The resulting eigenvalues and eigenfunctions are visualized as well as the reconstruction of the data via the full state observable.<jupyter_code># load and preprocess data
data = np.loadtxt(fname='MI_timesteps.txt', skiprows=1001)#skipping description + 1000 time steps (burn-in period)
print('Shape of whole data set: '+str(np.shape(data)))
T = data[:,0]#time
X1 = data[:,1]#measurement areas 1-9
#X2 = data[:,2]
#X3 = data[:,3]
#X4 = data[:,4]
#X5 = data[:,5]
#X6 = data[:,6]
#X7 = data[:,7]
#X8 = data[:,8]
#X9 = data[:,9]

# area 1 is used solitair
Y1 = X1[1:]# if X is timestep n, Y is timestep n+1
X1 = X1[:-1]#last element cannot be used
X1 = X1.reshape((np.size(X1),1))
Y1 = Y1.reshape((np.size(Y1),1))
print('Shape of X1: '+str(np.shape(X1)))
print('Shape of Y1: '+str(np.shape(Y1)))

# all areas together
all_X = data[:,1:]
all_Y = all_X[1:,:]
all_X = all_X[:-1,:]
print('Shape of all_X: '+str(np.shape(all_X)))
print('Shape of all_Y: '+str(np.shape(all_Y)))

# plot of the first area 
plt.scatter(T[:-1], X1, label='X original', s=5)
plt.scatter(T[1:], Y1, label='Y original', s=5)
plt.legend()
plt.show()
# apply EDMD to obtain eigenvalues, eigenfunctions, and modes
evals, efuncs, mos = EDMD(X1, Y1, D)
# only X1 - plot first 8 nontrivial eigenfunctions
x = np.linspace(-5,5,100)
y = np.linspace(-5,5,100)
X, Y = np.meshgrid(x, y)
fig, axs = plt.subplots(2,4, figsize=(12,5), constrained_layout=True)
cmap = mpl.cm.jet#colormap
for i in range(1,9):
    values = efuncs(X,Y)[i]
    ax = axs.flat[i-1]
    pcm = ax.pcolormesh(np.real(values), cmap=cmap)
    ax.set_xticks([0,25,50,75,100])
    ax.set_xticklabels(np.linspace(-5,5,5))
    ax.set_yticks([0,25,50,75,100])
    ax.set_yticklabels(np.linspace(-5,5,5))
    ax.set_title('$\mu_{:d}={:.3g}$'.format(i,evals[i]))
    ax.set_xlabel('x')
    ax.set_ylabel('y')


norm = mpl.colors.Normalize(vmin=-1, vmax=1)
plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), ax=axs[0,:])
plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), ax=axs[1,:])
plt.show()
# Reconstruction of X1 and Y1 using the Koopman decomposition
t = 1#one time step -> Y
x_rec = mos @ efuncs(X1,X1).T
y_rec = evals**t * mos @ efuncs(X1,X1).T
X_rec = np.array((x_rec[0]))
Y_rec = np.array((y_rec[0]))
# check for same shapes as above
print('Shape of X_rec: '+str(np.shape(X_rec)))
print('Shape of Y_rec: '+str(np.shape(Y_rec)))

#plot reconstructed X1 and Y1 against the time
plt.scatter(T[:-1], X_rec, label='X reconstructed', s=5)
plt.scatter(T[1:], Y_rec, label='Y reconstructed', s=5)
#plt.scatter(T[:-1], X1, label='X original')# too small compared to the reconstructed data
#plt.scatter(T[1:], Y1, label='Y original')
plt.legend()
plt.show()
# apply EDMD to obtain eigenvalues, eigenfunctions, and modes
all_evals, all_efuncs, all_mos = EDMD(all_X, all_Y, D)
# all areas - plot first 8 nontrivial eigenfunctions
x = np.linspace(-5,5,100)
y = np.linspace(-5,5,100)
X, Y = np.meshgrid(x, y)
fig, axs = plt.subplots(2,4, figsize=(12,5), constrained_layout=True)
cmap = mpl.cm.jet#colormap
for i in range(1,9):
    values = all_efuncs(X,Y)[i]
    ax = axs.flat[i-1]
    pcm = ax.pcolormesh(np.real(values), cmap=cmap)
    ax.set_xticks([0,25,50,75,100])
    ax.set_xticklabels(np.linspace(-5,5,5))
    ax.set_yticks([0,25,50,75,100])
    ax.set_yticklabels(np.linspace(-5,5,5))
    ax.set_title('$\mu_{:d}={:.3g}$'.format(i,all_evals[i]))
    ax.set_xlabel('x')
    ax.set_ylabel('y')


norm = mpl.colors.Normalize(vmin=-1, vmax=1)
plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), ax=axs[0,:])
plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), ax=axs[1,:])
plt.show()
# Reconstruction of X and Y for all measurement areas
fig, axs = plt.subplots(9,1, figsize=(5,27), constrained_layout=True)
# loop over every measurement area
for i in range(0,9):
    ax=axs.flat[i]
    t = 1#one time step -> Y
    x_rec = all_mos @ all_efuncs(all_X[:,i],all_X[:,i]).T
    y_rec = all_evals**t * all_mos @ all_efuncs(all_X[:,i],all_X[:,i]).T
    # Reshaping - little bit ugly, but works
    X_rec = np.array((x_rec[0], x_rec[1], x_rec[2], x_rec[3], x_rec[4], x_rec[5], x_rec[6], x_rec[7], x_rec[8]))
    Y_rec = np.array((y_rec[0], y_rec[1], y_rec[2], y_rec[3], y_rec[4], y_rec[5], y_rec[6], y_rec[7], y_rec[8]))
    # plot
    ax.scatter(T[:-1], X_rec[1,:], label='X reconstructed', s=5)
    ax.scatter(T[1:], Y_rec[0,:], label='Y reconstructed', s=5)
    ax.set_title('Measurement area {:d}'.format(i+1))
    ax.legend()
    
    #saving the images (if desired)
    #extent = ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted())
    #plt.savefig('Area_'+str(i+1)+'.png', bbox_inches=extent.expanded(1.1, 1.2))<jupyter_output><empty_output><jupyter_text>Attempt to use EDMD from datafold code:<jupyter_code>import csv
import pandas as pd

from mpl_toolkits.mplot3d import Axes3D
from datafold.pcfold import TSCDataFrame, GaussianKernel
from datafold.dynfold import DMDFull
from datafold.dynfold.transform import TSCRadialBasis, TSCPolynomialFeatures
from datafold.appfold import EDMD
df = pd.read_csv('MI_timesteps.txt', sep=" ")
df = df.drop(['timeStep'], axis=1)
df = df = df.iloc[1000:]
tsc_data = TSCDataFrame.from_single_timeseries(df)
df['areaDensityCounting-PID5'].plot()
df['areaDensityCounting-PID6'].plot()
df['areaDensityCounting-PID7'].plot()
tsc_data<jupyter_output><empty_output><jupyter_text>EDMD with identity dictionary:<jupyter_code>dmd = DMDFull().fit(X=tsc_data, store_koopman_matrix=True)
dmd_values = dmd.predict(tsc_data.initial_states(), time_values=None)

print("Data snipped with predicted time series data")
dmd_values
dmd_values['areaDensityCounting-PID5'].plot()
dmd_values['areaDensityCounting-PID6'].plot(color = 'green')
dmd_values['areaDensityCounting-PID7'].plot(color = 'orange')<jupyter_output><empty_output><jupyter_text>EDMD with dictionary of polynomials of degree 2:<jupyter_code>dict_step = [
    (
        "polynomial",
        TSCPolynomialFeatures(degree=2),
    )
]

edmd_poly = EDMD(dict_steps=dict_step, include_id_state=True).fit(X=tsc_data)
edmd_poly_values = edmd_poly.predict(tsc_data.initial_states())
# access models in the dictionary
print(edmd_poly.named_steps["polynomial"])

print("")
print("polynomial degrees for data")
print(edmd_poly.named_steps["polynomial"].powers_)

print("")
print("Dictionary space values:")
edmd_p = edmd_poly.reconstruct(tsc_data)#transform(tsc_data)
edmd_p
edmd_p['areaDensityCounting-PID5'].plot()
edmd_p['areaDensityCounting-PID6'].plot(color = 'green')
edmd_p['areaDensityCounting-PID7'].plot(color = 'orange')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16840971033783356, 'fraction_numerical': 0.036001774735374276, 'mean_word_length': 5.022137404580152, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '418345', 'license': 'MIT', 'n_tokens_mistral': 6659, 'n_tokens_neox': 6319, 'n_words': 1443, 'path': 'Final_project/EDMD.ipynb', 'repo_name': 'Schnursenkel/mlcms_group_g', 'size': '15777'}