-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathphasegrids_reku.py
219 lines (142 loc) · 6.48 KB
/
phasegrids_reku.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
# Code to generate the grids of phasemaps at different times for the 2D ReKU model.
# The code can generate 5x5 grids. Make changes to generate alternate sized grids.
# We use multiprocessing module in Python to parallelize the computation.
# This code is particularly fast and efficient when run on a multiple core CPU server.
#
#
# Author: Kaushik Roy
#
# Current affiliation: Dept. of Biosciences, Rice University, TX 77005, USA
# For correspondence, please email: [email protected]
#
#
# Import libraries
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
from scipy.stats import truncnorm
from multiprocessing import Pool
import time
# =============== Initialization Functions ===============
# Legacy seed for reproducibility
np.random.seed(12345)
# Define the omega distributions
# Uniform distribution
def omega_uniform(w_min, w_max, N): # N is the number of oscillators along each axis
omega = np.random.uniform(w_min,w_max, size=(N,N))
return omega
# Truncated normal distribution
def omega_trunc_normal(w_min, w_max, mean, scale, N):
# Define the parameters of the truncated normal distribution
a = (w_min - mean) / scale # Lower bound
b = (w_max - mean) / scale # Upper bound
# Generate random samples from the truncated Gaussian distribution
omega = truncnorm.rvs(a, b , loc = mean , scale = scale , size = (N,N))
return omega
# Homogenous initial phases
def phase_homogenous(N,angle): # N is the number of oscillators along each axis
phase = np.ones((N,N))*angle # the angle variable can be set to any value between -pi and pi
return phase
# Initial phase distributions of width less than 2*pi
# We use the following nomenclature for the distributions
# narrowp1: np.random.uniform(0, np.pi/2, (N,N))
# narrowp2: np.random.uniform(0, np.pi, (N,N))
# narrowp3: np.random.uniform(-np.pi/2, np.pi/2, (N,N))
# narrowp4: np.random.uniform(-np.pi/4, np.pi/4, (N,N))
# widerp1: np.random.uniform(-np.pi/8, np.pi, (N,N))
# widerp2: np.random.uniform(-np.pi/4, np.pi, (N,N))
# widerp3: np.random.uniform(-3*np.pi/8, np.pi, (N,N))
# widerp4: np.random.uniform(-np.pi/2, np.pi, (N,N))
def phase_narrow(N):
phase = np.random.uniform(-np.pi/2 , np.pi/2 , (N,N)) # narrowp3
return phase
# completely random initial phases
def phase_random(N):
phase = np.random.uniform(-np.pi, np.pi, (N, N)) # randomp
return phase
# =============== Model Functions ===============
# define the maximum function
def max_sin(x):
return np.maximum(np.sin(x), 0)
# 2D ReKU model
def ReKU_model_2d(phi,N, K, omega):
phi_up = np.roll(phi, -1, axis=0) # getting nearest neighbors
phi_up [-1, :]= phi[-1, :] # imposing open boundary conditions
phi_down = np.roll(phi, 1, axis=0)
phi_down[0, :] = phi[0, :]
phi_left = np.roll(phi, -1, axis=1)
phi_left[:, -1] = phi[:, -1]
phi_right = np.roll(phi, 1, axis=1)
phi_right[:, 0] = phi[:, 0]
dphi_dt = np.zeros((N, N))
dphi_dt += max_sin(phi_up-phi) + max_sin(phi_down-phi)
dphi_dt += max_sin(phi_left-phi) + max_sin(phi_right-phi)
dphi_dt = omega + K * dphi_dt # computing the phase velocity
return dphi_dt
# =============== Simulation Functions ===============
# function that computes the phases at the end of the time integration
def phase_map_ReKU(N, T, dt, K, phi_initial, omega_initial):
num_steps = int(T / dt) + 1
X = np.zeros((num_steps, N, N)) # Storing the phases
omega = omega_initial.copy() # ensuring that the same initial
phi = phi_initial.copy() # conditions are used throughout
for t in range(num_steps):
k1 = ReKU_model_2d(phi, N, K, omega)
k2 = ReKU_model_2d(phi + 0.5 * dt * k1, N, K, omega)
k3 = ReKU_model_2d(phi + 0.5 * dt * k2, N, K, omega)
k4 = ReKU_model_2d(phi + dt * k3, N, K, omega)
phi += (dt / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
phi = np.angle(np.exp(1j * phi)) # Ensuring that phi is between -pi and pi
X[t, :, :] = phi
return X[-1, :, :]
# =============== Plotting and Utility Functions ===============
# Wrap the phase map function to allow passing a single tuple of parameters
def generate_plot_data_ReKU(args):
N, T, dt, K, phi_initial, omega_initial = args
return phase_map_ReKU(N, T, dt, K, phi_initial, omega_initial)
# Define the model parameters
w_min = 2*np.pi/180
w_max = 2*np.pi/150
# for use with truncated normal distributions
#mean = 2*np.pi/165
#scale = 2*np.pi*(1/165 - 1/160)
# for use with homogenous initial phases
#in_phase = np.pi/4
N = 50 # number of oscillators: N^2
dt = 0.01 # time step
# adjust the range as necessary
a_values = np.linspace(0.5, 3.0, 25) # range of a in K = a \Delta_{\omega}
# define the initial frequency and phase distributions for the computation
omega_initial = omega_uniform(w_min, w_max, N) # uniform distribution
# uncomment below if using truncated normal distribution
#omega_initial = omega_trunc_normal(w_min, w_max, mean, scale, N)
phi_initial = phase_narrow(N) # change the limits in original function for different results
# uncomment below for widest possible initial phase distribution
#phi_initial = phase_random(N)
start_time = time.time()
# Define a list of desired T values
T_values = [1000, 1200, 1500]
# Create a PDF file for the plots at the current T value
# customize output filename as necessary
pdf_filename = "phasemapsreKU_narrowp3_uniformf_5x5.pdf"
pdf_pages = PdfPages(pdf_filename)
for T in T_values:
# Prepare a list of parameters for each iteration
params = [(N, T, dt, a * (w_max - w_min) , phi_initial, omega_initial) for a in a_values]
# Use multiprocessing to generate plot data
with Pool() as p:
plot_data = p.map(generate_plot_data_ReKU, params)
# Prepare the plot for the current T value
fig, axes = plt.subplots(5, 5, figsize=(20, 20))
for idx, data in enumerate(plot_data):
ax = axes[idx // 5, idx % 5] # Get the current axis
ax.imshow(data, cmap="plasma", origin="lower", extent=[0, N, 0, N])
ax.axis("off")
ax.set_title(f'a={a_values[idx]:.2f}', fontsize=16)
plt.tight_layout()
pdf_pages.savefig(fig, bbox_inches='tight')
plt.close()
pdf_pages.close()
print(f"Phase maps saved to {pdf_filename}")
end_time = time.time()
print("CPU computation took", end_time - start_time, "seconds")