-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgenerate_data.py
418 lines (355 loc) · 12.1 KB
/
generate_data.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
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
## Written by Madeline Galbraith
## Last edited: July 2022
import numpy as np
import pandas as pd
## Allow the calculation to be parallelized
from mpi4py import MPI
myrank = MPI.COMM_WORLD.Get_rank()
nprocs = MPI.COMM_WORLD.Get_size()
stat = MPI.Status()
###################
###################
###################
## Parameters for notch system
N0 = 5.0e+2 ##molecules/hour
D0 = 1.0e+3 ##molecules/hour
I0 = 2.0e+2 ##molecules
kc = 5.0e-4 #1/hour/molecule
kt = 5.0e-5 #1/hour/molecule
nI = 2.0e+0 #hill coeff
lambdaN = 2.0e+0 # no dimension
lambdaD = 0.0e+0 # no dimension
gamma = 1.0e-1 #1/hour
gamma_i = 5.0e-1 #1/hour
###################
###################
###################
def generate_lattice(nrow,ncol,nlayer,special,maxICS):
### Create the initial conditions for
### -'1D'
### 2D square lattice
### nlayer >1 means a 2d system
### nrow and ncol is the length and width of system
nlayer=1
a= np.zeros((nlayer,nrow,ncol))
if not special:
a= np.random.uniform(0,maxICS,[nlayer,nrow,ncol])
elif special=='SP':
a[nlayer/2][nrow/2][ncol/2]=3000
return a
###################
###################
###################
def repackage(total):
## rebuild dictionary so its easier to parse for
test = total[0]
keys = test.keys()
update={}
for k in keys:
update[k]=[]
for i in range(len(total)):
for k in keys:
update[k]+=[total[i][k]]
return update
###################
###################
###################
def shiftedHill(X,X0,nX,lambdY):
return lambdY+(1.0 - lambdY)/(1.0+(X/X0)**nX)
###################
###################
###################
def ts_update(dt,X,tempX,noiseType,noiseX):
if noiseType=='shot' and noiseX!=0:
sigma=np.sqrt(dt*X)*np.random.normal(0,noiseX,X.shape)
elif noiseType=='white' and noiseX!=0:
sigma=np.random.normal(0,noiseX,X.shape)*dt
else:
sigma=0.
X = X+dt*tempX+sigma
decX = np.sign(X)+1 ## decX is now 0 if X neg, 1 if X=0, and 2 if X pos
X = X*decX/2. ## zeros out neg X vals, doesn't modif if X==0 or X is positive
return X
###################
###################
###################
def getExt(X,lattice):
if '1D' in lattice:
## the value of Next and Dext is equal to the value of N or D in the neighbor
if X.shape[1]>1:
## if the array has two columns
return np.roll(X,1,axis=1)
elif X.shape[2]>1:
## if the array has two rows
return np.roll(X,1,axis=2)
elif lattice=='2Dsquare':
## the value of Next and Dext is equal to the 1/4 of N or D in the neighbor as each cell has 4 neighbors
return 1./4.*(np.roll(X,1,axis=1)+np.roll(X,-1,axis=1)+np.roll(X,1,axis=2)+
np.roll(X,-1,axis=2))
###################
###################
###################
def run_simulation(tmax,delta_t,lattice,latticesize,noiseType,noiseAmp,special=None,ICS=None,step4save=None):
if not step4save:
step4save=1
if not ICS:
## get the initial conditions
## if special is not set than conditions are randomized
Ni=generate_lattice(latticesize,latticesize,1,special,10000)
Di=generate_lattice(latticesize,latticesize,1,special,10000)
Ii=generate_lattice(latticesize,latticesize,1,special,2000)
varbls = {'N':Ni,'D':Di,'I':Ii}
else:
varbls=ICS
total=[]
for i in range(int(tmax/delta_t)):
## Nc,Dc, and Ic are the Notch, Delta, and NICD values at time=t
Nc,Dc,Ic = varbls['N'],varbls['D'],varbls['I']
## calculate the value of Next and Dext, depends on whether the system is 2 cell or multicell
Next = getExt(Nc,lattice)
Dext = getExt(Dc,lattice)
## temp_X =dX/dt, must multiply by dt to get the dX which will be added to value at time=t
temp_N = N0*shiftedHill(Ic,I0,nI,lambdaN) -Nc*(kc*Dc+gamma+kt*Dext)
temp_D = D0*shiftedHill(Ic,I0,nI,lambdaD) -Dc*(kc*Nc+gamma+kt*Next)
temp_I = kt*Nc*Dext-gamma_i*Ic
## get value of N,D, and I at time=t+dt
## also include the noise calculation
Nc = ts_update(delta_t,Nc,temp_N,noiseType,noiseAmp)
Dc = ts_update(delta_t,Dc,temp_D,noiseType,noiseAmp)
Ic = ts_update(delta_t,Ic,temp_I,noiseType,0)## use noiseI=0 because we don't add noise here
varbls={'N':Nc,'D':Dc,'I':Ic}
if i%step4save==0:
total.append(varbls)
total = repackage(total)
return total
###################
###################
###################
#-----------------------------------
def add(arr,val):
for el in arr:
if el[0]==val[0] and el[1]==val[1]:
return False
return True
def flip(ics,rv,cv):
rv = rv[0]
cv = cv[0]
maxICS=10000
for key in ['I','N','D']:
ics[key][0][rv][cv]=np.random.uniform(0,maxICS,1)[0]
return ics
#-----------------------------------
def getLatticeICS(pattern,recv,senv,shapeV):
ics={}
ics['N']=np.zeros(shapeV)
ics['I']=np.zeros(shapeV)
ics['D']=np.zeros(shapeV)
if pattern==0:
## only one corner is receivers
for k in ics:
ics[k][0,:8,:8]=recv[k]
ics[k][0,:8,8:]=senv[k]
ics[k][0,8:]=senv[k]
elif pattern==1:
## only one corner is senders
for k in ics:
ics[k][0,:8,:8]=senv[k]
ics[k][0,:8,8:]=recv[k]
ics[k][0,8:]=recv[k]
elif pattern==2:
## all senders
for k in ics:
ics[k][0,:,:]=senv[k]
elif pattern==3:
## all receivers
for k in ics:
ics[k][0,:,:]=recv[k]
elif pattern==4:
## one line receivers
for k in ics:
ics[k][0,:,:]=senv[k]
ics[k][0,2,:]=recv[k]
elif pattern==5:
## one line senders
for k in ics:
ics[k][0,:,:]=recv[k]
ics[k][0,2,:]=senv[k]
elif pattern==6:
## half senders
for k in ics:
ics[k][0,:8,:]=senv[k]
ics[k][0,8:,:]=recv[k]
elif pattern==7:
## nucleating site
randx,randy = np.random.randint(0,16,2)
for k in ics:
ics[k][0,:,:]=recv[k]
ics[k][0,randx,randy]=senv[k]
return ics
#-----------------------------------
##########################################
############ Start simulation ###########
##########################################
## Code will work with or without mpi
## setup to run the different types of noise on different processes
nAmpList={'white':{},'shot':{}}
nTypeL={}
if nprocs==1:
nTypeL[0]=['white','shot']
nAmpList['white'][0] = list(np.arange(0,2100,100))
nAmpList['shot'][0] = list(np.arange(0,21,1))
elif nprocs<=3:
nTypeL[0]=['white']
nTypeL[1]=['shot']
nAmpList['white'][0] = list(np.arange(0,2100,100))
nAmpList['shot'][1] = list(np.arange(0,21,1))
elif nprocs>=4:
nTypeL[0]=['white']
nTypeL[1]=['white']
nTypeL[2]=['shot']
nTypeL[3]=['shot']
nAmpList['white'][0] = list(np.arange(0,1100,100))
nAmpList['white'][1] = list(np.arange(1100,2100,100))
nAmpList['shot'][2] = list(np.arange(0,11,1))
nAmpList['shot'][3] = list(np.arange(11,21,1))
if myrank not in nTypeL.keys():
## if the process won't do anything than exit to save computation time
exit()
## get seeds so data is consistent across runs
seedList=pd.read_csv(filepath_or_buffer='seedsForSims.txt',header=None).values[:,0]
## parameters for simulation
mistakeList=[13,64,230]
sdevList=[50,1500,5000]
pattList=[0,7]#0-7 for patterns listed above]
latsizeList=[16]
tmax=10000
dt =0.1
numberOfSims=20
lattice='2Dsquare'
## Set these flags to determine which simulation(s) you want to run
genMistakes=False
genDeviations=False
genPatterns=False
genRandom=True
for latticesize in latsizeList:
for nType in nTypeL[myrank]:
for nAmp in nAmpList[nType][myrank]:
if genMistakes:
for mistakes in mistakeList:
for seedInd in range(numberOfSims):
ss = seedList[seedInd]
np.random.seed(ss)
nrow=latticesize
ncol=latticesize
## Initial condition: Salt and pepper with discrete perturbations (mistakes in the lattice)
## First generate the salt and pepper pattern
datsp=run_simulation(tmax,dt,lattice,latticesize,nType,0,special='SP')
##
ics={}
for k in datsp:
ics[k] = datsp[k][-1]
v=[]
while len(v)<mistakes:
newEl = np.array([np.random.randint(0,nrow,1),np.random.randint(0,nrow,1)])
t=add(v,newEl)
if t:
v+=[newEl]
v = np.array(v)
for k in range(len(v)):
ics = flip(ics,v[k][0],v[k][1])
## solve the model
dat1=run_simulation(tmax,dt,lattice,latticesize,nType,nAmp)
titleT = "data/mistakes/traj_"+str(nrow)+"x"+str(ncol)+"x1_"+str(nType)+"_n"+str(nAmp)+"_m"+str(mistakes)+"_s"+str(seedInd)+".dat"
##Put the data into a form that is easier to output
for k in dat1:
for i in range(len(dat1[k])):
dat1[k][i] = list(np.around(dat1[k][i].flatten(),decimals=3))
df_Final = pd.DataFrame(dat1)
df_Final['run'] = seedInd
df_Final.to_csv(path_or_buf=titleT,index=False)
elif genDeviations:
for sdev in sdevList:
for seedInd in range(numberOfSims):
ss = seedList[seedInd]
np.random.seed(ss)
nrow=latticesize
ncol=latticesize
## Initial condition: Salt and pepper with continuous perturbations (Gaussian added to lattice)
## First generate the salt and pepper pattern
datsp=run_simulation(tmax,dt,lattice,latticesize,nType,0,special='SP')
ics={}
for k in datsp:
ics[k] = datsp[k][-1]
ics[k]=np.random.normal(ics[k],sdev,ics[k].shape)
ics[k] = (ics[k]>=0)*ics[k]
## solve the model
dat1=run_simulation(tmax,dt,lattice,latticesize,nType,nAmp)
titleT = "data/dev/traj_"+str(nrow)+"x"+str(ncol)+"x1_"+str(nType)+"_n"+str(nAmp)+"_d"+str(sdev)+"_s"+str(seedInd)+".dat"
##Put the data into a form that is easier to output
for k in dat1:
for i in range(len(dat1[k])):
dat1[k][i] = list(np.around(dat1[k][i].flatten(),decimals=3))
df_Final = pd.DataFrame(dat1)
df_Final['run'] = seedInd
df_Final.to_csv(path_or_buf=titleT,index=False)
elif genPatterns:
datsp=run_simulation(tmax,dt,lattice,latticesize,nType,0,special='SP')
recv ={'N':np.max(datsp['N'][-1]),'D':np.min(datsp['D'][-1]),'I':np.max(datsp['I'][-1])}
senv ={'N':np.min(datsp['N'][-1]),'D':np.max(datsp['D'][-1]),'I':np.min(datsp['I'][-1])}
shapeV = datsp['N'][-1].shape
for pattN in pattList: ## >=0 for pattern
for seedInd in range(numberOfSims):
ss = seedList[seedInd]
np.random.seed(ss)
nrow=latticesize
ncol=latticesize
## Generate the pattern:
## 0: only one corner is receivers, 1: only one corner is senders
## 2: all senders, 3: all receivers
## 4: one line receivers, 5: one line senders
## 6: half senders
## 7: nucleating condition
## First generate the salt and pepper pattern
ics=getLatticeICS(pattN,recv,senv,shapeV)
## solve the model
dat1=run_simulation(tmax,dt,lattice,latticesize,nType,nAmp)
titleT = "data/pattern/traj_"+str(nrow)+"x"+str(ncol)+"x1_"+str(nType)+"_n"+str(nAmp)+"_p"+str(pattN)+"_s"+str(seedInd)+".dat"
##Put the data into a form that is easier to output
for k in dat1:
for i in range(len(dat1[k])):
dat1[k][i] = list(np.around(dat1[k][i].flatten(),decimals=3))
df_Final = pd.DataFrame(dat1)
df_Final['run'] = seedInd
df_Final.to_csv(path_or_buf=titleT,index=False)
if latticesize !=16:
for seedInd in range(numberOfSims):
ss = seedList[seedInd]
np.random.seed(ss)
nrow=latticesize
ncol=latticesize
## solve model for randomized initial conditions
dat1=run_simulation(tmax,dt,lattice,latticesize,nType,nAmp)
titleT = "data/lattices/traj_"+str(nrow)+"x"+str(ncol)+"x1_"+nType+"_n"+str(nAmp)+"_s"+str(seedInd)+".dat"
##Put the data into a form that is easier to output
for k in dat1:
for i in range(len(dat1[k])):
dat1[k][i] = list(np.around(dat1[k][i].flatten(),decimals=3))
df_Final = pd.DataFrame(dat1)
df_Final['run'] = seedInd
df_Final.to_csv(path_or_buf=titleT,index=False)
elif genRandom:
for seedInd in range(numberOfSims):
ss = seedList[seedInd]
np.random.seed(ss)
nrow=latticesize
ncol=latticesize
## solve model for randomized initial conditions
dat1=run_simulation(tmax,dt,lattice,latticesize,nType,nAmp)
titleT = "data/random/traj_"+str(nrow)+"x"+str(ncol)+"x1_"+nType+"_n"+str(nAmp)+"_s"+str(seedInd)+".dat"
##Put the data into a form that is easier to output
for k in dat1:
for i in range(len(dat1[k])):
dat1[k][i] = list(np.around(dat1[k][i].flatten(),decimals=3))
df_Final = pd.DataFrame(dat1)
df_Final['run'] = seedInd
df_Final.to_csv(path_or_buf=titleT,index=False)