Warning:  Due to planned infrastructure maintenance, the EBRAINS Wiki and EBRAINS Support system will be unavailable for up to three days starting Monday, 14 July. During this period, both services will be inaccessible, and any emails sent to the support address will not be received.

Attention: We are currently experiencing some issues with the EBRAINS Drive. Please bear with us as we fix this issue. We apologise for any inconvenience caused.


Wiki source code of Code description

Last modified by galluzziandrea on 2022/06/20 12:33
Show last authors
1 == Introduction (path and modules): ==
2
3 First of all we check the path and import the necessary modules .
4
5 === Check where I am and place myself in the right folder: ===
6
7 {{code language="python"}}
8 # Import the os module
9 import os
10
11 #Print the current working directory
12 print("Current working directory: {0}".format(os.getcwd()))
13
14 # Change the current working directory
15 os.chdir('/mnt/user/shared/Slow waves in fading anesthesia/Nest3Python3')
16
17 # Print the current working directory
18 print("Current working directory: {0}".format(os.getcwd()))
19 {{/code}}
20
21 === Import the modules necessary for the simulation: ===
22
23 {{code language="python"}}
24 import nest
25 import time
26 from numpy import exp
27 import numpy
28 import math
29 import random
30 import multiprocessing
31 Inizio=time.time()
32 print('tempo di Inizio:',Inizio)
33 {{/code}}
34
35 === ===
36
37 === Define necessary classes to import the Initialization Files: ===
38
39 {{code language="python" title=" "}}
40 class ImportIniLIFCA():
41 #initialize the information to look for in perseo.ini
42 inf=["NeuronType", #still fixed value
43 "DelayDistribType", #still fixed value
44 "SynapticExtractionType", #still fixed value
45 "Life"]
46
47 def __init__(self,files):
48 self.files=files
49
50 def FilesControllo(self):
51 import sys
52 for i in range(0,len(self.files)):
53 if self.FileControllo(self.files[i]):
54 sys.exit(0)
55
56 def FileControllo(self,file1):
57 try:
58 f1=open(file1,"r")
59 f1.close()
60 return 0
61 except ValueError:
62 print("ValueError")
63 return 1
64 except IOError as err:
65 print("OS error: {0}".format(err))
66 return 1
67 except:
68 print("Unexpected error:", sys.exc_info()[0])
69 return 1
70
71 def Estrai_inf(self,stampa=0):
72
73 InfoPerseo=self.EstraiInfoPerseo() #extract info from perseo.ini
74 AppoggioTempM=self.EstraiInfoModuli() #extract info from modules.ini
75 AppoggioTempC=self.EstraiInfoConnectivity() #extract info from connectivity.ini
76 AppoggioTempP=self.EstraiProtocol() #extract info from protocol.ini
77
78 def getKey(item):
79 return item[0]
80 InfoProtocol=AppoggioTempP
81 # I convert the extracted information into a suitable format from tuple to list
82
83 InfoBuildT=[AppoggioTempM[0]]
84 for i in range(0,AppoggioTempM[0]):
85 app1=[int(AppoggioTempM[2][i][0])]
86 app=(app1+list(AppoggioTempM[2][i][3:9])+list(AppoggioTempM[2][i][12])+list(AppoggioTempM[2][i][9:12]))
87 InfoBuildT.append(app)
88 del app
89
90 InfoBuild=[float(InfoBuildT[0])]
91 for i in range(0,int(InfoBuildT[0])):
92 app=[]
93 for j in range(0,11):
94 app.append(float(InfoBuildT[i+1][j]))
95 InfoBuild=InfoBuild+[app]
96 del app
97
98 InfoConnectPop=[AppoggioTempM[0]]
99 for i in range(0,len(AppoggioTempC[1][:])):
100 app=list(AppoggioTempC[1][i])
101 InfoConnectPop.append(app)
102 del app
103
104 InfoConnectNoise=[AppoggioTempM[0]]
105 for i in range(0,AppoggioTempM[0]):
106 app=list(AppoggioTempM[2][i][1:3])
107 InfoConnectNoise.append(app)
108
109
110 if stampa==1: #Print on screen of saved data
111 for i,j in enumerate(InfoPerseo):
112 print(self.inf[i],"=",j)
113 print("\n")
114 print("the network consists of ", AppoggioTempM[0], " neuronal population" )
115 print(AppoggioTempM[1])
116 for i in range(0,AppoggioTempM[0]):
117 print(AppoggioTempM[2][i])
118 print("\n")
119 print(AppoggioTempC[0])
120 for i in range(0,AppoggioTempM[0]**2):
121 print(AppoggioTempC[1][i])
122 print("\n")
123 for i in InfoProtocol:
124 print("SET_PARAM"+str(i))
125
126
127 return InfoPerseo,InfoBuild,InfoConnectPop,InfoConnectNoise,InfoProtocol
128
129 def EstraiProtocol(self):
130 import string
131 f1=open(self.files[3],"r")
132 ProtocolList= []
133 for x in f1.readlines():
134 y=x.split()
135 if len(y):
136 if x[0]!="#" and y[0]=="SET_PARAM":
137 try:
138 ProtocolList.append([float(y[1]),int(y[2]),float(y[3]),float(y[4])])
139 except ValueError:
140 pass
141 f1.close()
142 return ProtocolList
143
144 def EstraiInfoPerseo(self):
145 import string
146 f1=open(self.files[0],"r")
147 InfList= []
148 for x in f1.readlines():
149 y=x.split()
150 if len(y):
151 if x[0]!="#":
152 for findinf in self.inf:
153 try:
154 temp=y.index(findinf)
155 InfList.append(y[temp+2])
156 except ValueError:
157 pass
158 f1.close()
159 return InfList
160
161 def EstraiInfoModuli(self):
162 import string
163 f1=open(self.files[2],"r")
164 NumPop=0
165 for i,x in enumerate(f1.readlines()):
166 y=x.split()
167 if len(y):
168 if x[0]!="#":
169 NumPop=NumPop+1
170 if i==2:
171 ParamList=[]
172 for j in range(1,14):
173 ParamList.append(y[j])
174 f1.close()
175 PopsParamList=[]
176 f1=open(self.files[2],"r")
177 x=f1.readlines()
178 for j in range(0,NumPop):
179 appo=x[4+j];
180 PopsParamList.append(appo.split())
181 f1.close()
182 return NumPop,ParamList,PopsParamList
183
184 def EstraiInfoConnectivity(self):
185 import string
186 f1=open(self.files[1],"r")
187 PopConParamList=[]
188 for i,x in enumerate(f1.readlines()):
189 y=x.split()
190 if len(y):
191 if x[0]!="#":
192 PopConParamList.append(y)
193 if i==1:
194 ParamList=[]
195 for j in range(1,9):
196 ParamList.append(y[j])
197 f1.close()
198 return ParamList,PopConParamList
199 {{/code}}
200
201 === Import the initialization files: ===
202
203 in this section we...
204
205 {{code language="python"}}
206 Salva=1
207 file1="perseo35.ini"
208 file2="c_cortsurf_Pot1.43PotStr148v3.ini"
209 file3="m_cortsurf_Pot1.43.ini"
210 file4="ProtocolExploration36.ini"
211 files=[file1,file2,file3,file4]
212 #define the name of the Output file
213 FileName="dati/Rates_Nest_Run_Milano_Test36_13x13_"+str(nest.Rank())+"_Pot1.43PotStr148v3Long3.dat"
214 #check the existence of the files being read
215 ImpFil=ImportIniLIFCA(files);
216 ImpFil.FilesControllo()
217
218 #extract the information of interest from the files.ini and transfer them to the files:
219 #InfoPerseo,InfoBuild,InfoConnectPop,InfoConnectNoise
220
221 stampa=0; #stampa=1 print output simulation data on screen stampa=0 dont
222 InfoPerseo,InfoBuild,InfoConnectPop,InfoConnectNoise,InfoProtocol=ImpFil.Estrai_inf(stampa)
223
224 # InfoPerseo=["NeuronType","DelayDistribType","SynapticExtractionType","Life" ]
225 # InfoBuild=[numero di popolazioni,
226 # [N,C_ext,\nu_ext,\tau,\tetha,H,\tau_arp,NeuronInitType,\alpha_c,\tau_c,g_c],
227 # [.....],[],...]
228 # InfoConnectPop=[numero di popolazioni,
229 # [post,pre,c,Dmin,Dmax,syn typ,J,DJ],
230 # [.....],[],...]
231 # InfoConnectNoise=[numero di popolazioni,
232 # [J_ext,DJ_ext],
233 # [.....],[],...]
234 # InfoProtocol=[[time,population,param_num,value],
235 # [.....],[],...]
236 {{/code}}
237
238 === [[image:image-20220127173048-1.png||height="508" width="948"]] ===
239
240 === Defining general and nest.kernel parameters ===
241
242 {{code language="python"}}
243 #############################------------------------------------------------------------------------
244 #Clean the Network
245 #############################------------------------------------------------------------------------
246 nest.ResetKernel()
247
248 #############################------------------------------------------------------------------------
249 #insert the introductory parameters of the simulation
250 #############################------------------------------------------------------------------------
251
252
253 dt = 0.1 # the resolution in ms
254 StartMisure=0. # start time of measurements
255 simtime = int(float(InfoPerseo[3])) # Simulation time in ms (200 s)
256 if simtime<=StartMisure: # If the simulation time is less than StartMisure, it is increased by StartMisure
257 simtime=simtime+StartMisure
258 start=0.0 # start time of poissonian processes
259 origin=0.0 # temporal origin
260
261 #############################------------------------------------------------------------------------
262 # Kernel parameters
263 #############################------------------------------------------------------------------------
264 LNT=multiprocessing.cpu_count();
265 nest.SetKernelStatus({"local_num_threads": LNT})
266 nest.SetKernelStatus({"resolution": dt, "print_time": True,
267 "overwrite_files": True})
268
269 #############################------------------------------------------------------------------------
270 #"randomize" the seeds of the random generators
271 #############################------------------------------------------------------------------------
272
273 #msd = int(math.fabs(time.process_time()*1000))
274 #N_vp = nest.GetKernelStatus(['total_num_virtual_procs'])[0]
275 #pyrngs = [numpy.random.RandomState(s) for s in range(msd, msd+N_vp)]
276 #nest.SetKernelStatus({"grng_seed" : msd+N_vp})
277 #nest.SetKernelStatus({"rng_seeds" : list(range(msd+N_vp+1, msd+2*N_vp+1))})
278 {{/code}}
279
280 === Building the network: neuronal populations , Poisson processes and spike detectors ===
281
282 {{code language="python"}}
283 #############################------------------------------------------------------------------------
284 print("Building network")
285 #############################------------------------------------------------------------------------
286
287 startbuild = time.time() #initialize the calculation of the time used to simulate
288
289 NeuronPop=[]
290 NoisePop=[]
291 DetectorPop=[]
292
293 #define and initialize the populations of neurons with the parameters extracted from the.ini files
294 for i in range(1,int(InfoBuild[0])+1):
295 if int(InfoBuild[i][7])==0:
296 app=float(InfoBuild[i][5])
297 else:
298 app=0.
299 app2= nest.Create("aeif_psc_exp", int(InfoBuild[i][0]),params={"C_m": 1.0,
300 "g_L": 1.0/float(InfoBuild[i][3]),
301 "t_ref": float(InfoBuild[i][6]),
302 "E_L": 0.0,
303 "V_reset": float(InfoBuild[i][5]),
304 "V_m": app,
305 "V_th": float(InfoBuild[i][4]),
306 "Delta_T": 0.,
307 "tau_syn_ex": 1.0,
308 "tau_syn_in": 1.0,
309 "a": 0.0,
310 "b": float(InfoBuild[i][10]),
311 "tau_w": float(InfoBuild[i][9]),
312 "V_peak":float(InfoBuild[i][4])+10.0})
313 NeuronPop.append(app2)
314
315 #define and initialize the poisson generators and the spike detectors with the parameters extracted from the.ini files
316
317 for i in range(1,int(InfoBuild[0])+1):
318 app3= nest.Create("poisson_generator",params={"rate": float(InfoBuild[i][1]*InfoBuild[i][2]),
319 'origin':0.,
320 'start':start})
321 NoisePop.append(app3)
322 app4 = nest.Create("spike_recorder",params={ "start":StartMisure})
323 DetectorPop.append(app4)
324
325 endbuild = time.time()
326 {{/code}}
327
328 === [[image:image-20220127165908-2.png||height="659" width="1149"]] ===
329
330 === Connecting the network nodes: neuronal populations, Poisson processes and spike detectors ===
331
332 {{code language="python"}}
333 #############################------------------------------------------------------------------------
334 print("Connecting ")
335 #############################------------------------------------------------------------------------
336
337 startconnect = time.time()
338 Connessioni=[]
339 Medie=[]
340
341 #create and define the connections between the populations of neurons and the poisson generators
342 #and between the populations of neurons and the spike detectors with the parameters extracted from the.ini files
343
344 for i in range(0,int(InfoBuild[0])):
345 nest.Connect(NoisePop[i], NeuronPop[i], syn_spec={'synapse_model': 'static_synapse_hpc',
346 'delay': dt,
347 'weight': nest.math.redraw(nest.random.normal(mean=float(InfoConnectNoise[i+1][0]),
348 std=(float(InfoConnectNoise[i+1][1])*float(InfoConnectNoise[i+1][0]))),
349 min=0., max=float('Inf'))
350 })
351 nest.Connect(NeuronPop[i][:int(InfoBuild[i+1][0])], DetectorPop[i], syn_spec={"weight": 1.0, "delay": dt})
352
353 #create and define the connections between the populations of neurons with the parameters extracted from the.ini files
354
355 for i in range(0,len(InfoConnectPop[1:])):
356
357 conn=nest.Connect(NeuronPop[int(InfoConnectPop[i+1][1])], NeuronPop[int(InfoConnectPop[i+1][0])],
358 {'rule': 'pairwise_bernoulli',
359 'p':float(InfoConnectPop[i+1][2]) },
360 syn_spec={'synapse_model': 'static_synapse_hpc',
361 'delay':nest.math.redraw(nest.random.exponential(beta=float(1./(2.99573227355/(float(InfoConnectPop[i+1][4])-float(InfoConnectPop[i+1][3]))))),
362 min= numpy.max([dt,float(1./float(InfoConnectPop[i+1][4]))]),
363 max= float(1./(float(InfoConnectPop[i+1][3])-dt/2))),
364
365 'weight':nest.random.normal(mean=float(InfoConnectPop[i+1][6]),
366 std=math.fabs(float(InfoConnectPop[i+1][6])*float(InfoConnectPop[i+1][7])))})
367
368
369 endconnect = time.time()
370 {{/code}}
371
372 === ===
373
374 === ===
375
376 === [[image:image-20220127170722-1.png]] ===
377
378 === Simulating: neuronal time evolution. ===
379
380 === ===
381
382 {{code language="python"}}
383 #############################------------------------------------------------------------------------
384 print("Simulating")
385 #############################------------------------------------------------------------------------
386 ###################################################################################################################################################################
387 if Salva:
388 print("I m going to save the data")
389 #x=str(iterazioni)
390 f = open(FileName,"w")
391 if len(InfoProtocol):
392 print("I m going to split the simulation")
393 tempo=0
394 for contatore in range(0,len(InfoProtocol)):
395 appoggio1=int((tempo+InfoProtocol[contatore][0])/1000.)
396 appoggio2=int(tempo/1000.)
397 appoggio3=tempo+InfoProtocol[contatore][0]
398 if (appoggio1-appoggio2)>=1:
399 T1=(1+appoggio2)*1000-tempo
400 nest.Simulate(T1)
401 #Save the Data!!!!
402 ###########################################################
403 Equilibri=[]
404 for i in range(0,int(InfoBuild[0])):
405 Equilibri.append([])
406 a=nest.GetStatus(DetectorPop[i])[0]["events"]["times"]
407 if len(a)>0:
408 Trange=(1000*int(numpy.min(a)/1000.),1000*int(numpy.min(a)/1000.)+1000)
409 hist,Tbin=numpy.histogram(a,200,(Trange[0],Trange[1]))
410 Equilibri[i]=hist*1000./(5.*int(InfoBuild[i+1][0]))
411 else:
412 Trange=(1000*int(tempo/1000.),1000*int(tempo/1000.)+1000)
413 hist=numpy.zeros(200)
414 Tbin=numpy.linspace(Trange[0],Trange[1],num=201)
415 Equilibri[i]=hist
416 nest.SetStatus(DetectorPop[i],{'n_events':0})
417 for j in range(0,len(hist)):
418 f.write(str(Tbin[j])+" ")
419 for i in range(0,int(InfoBuild[0])):
420 f.write(str(Equilibri[i][j])+" ")
421 f.write("\n ")
422 ###########################################################
423 tempo=tempo+T1
424 for contatore2 in range(1,(appoggio1-appoggio2)):
425 nest.Simulate(1000.)
426 #Save the Data!!!!
427 ###########################################################
428 Equilibri=[]
429 for i in range(0,int(InfoBuild[0])):
430 Equilibri.append([])
431 a=nest.GetStatus(DetectorPop[i])[0]["events"]["times"]
432 if len(a)>0:
433 Trange=(1000*int(numpy.min(a)/1000.),1000*int(numpy.min(a)/1000.)+1000)
434 hist,Tbin=numpy.histogram(a,200,(Trange[0],Trange[1]))
435 Equilibri[i]=hist*1000./(5.*int(InfoBuild[i+1][0]))
436 else:
437 Trange=(1000*int(tempo/1000.),1000*int(tempo/1000.)+1000)
438 hist=numpy.zeros(200)
439 Tbin=numpy.linspace(Trange[0],Trange[1],num=201)
440 Equilibri[i]=hist
441 nest.SetStatus(DetectorPop[i],{'n_events':0})
442 for j in range(0,len(hist)):
443 f.write(str(Tbin[j])+" ")
444 for i in range(0,int(InfoBuild[0])):
445 f.write(str(Equilibri[i][j])+" ")
446 f.write("\n ")
447 tempo=tempo+1000.
448 T2=appoggio3-tempo
449 nest.Simulate(T2);
450 tempo=tempo+T2;
451 else:
452 nest.Simulate(InfoProtocol[contatore][0])
453 temp=InfoProtocol[contatore][0]
454 tempo=tempo+temp
455 if InfoProtocol[contatore][2]==4:
456 nest.SetStatus(NoisePop[InfoProtocol[contatore][1]],params={"rate": float(InfoBuild[1+InfoProtocol[contatore][1]][2]*InfoProtocol[contatore][3])})
457 if InfoProtocol[contatore][2]==12:
458 nest.SetStatus(NeuronPop[InfoProtocol[contatore][1]], params={"b": float(InfoProtocol[contatore][3])})
459 else:
460 nest.Simulate(simtime)
461 tempo=simtime
462 if (simtime-tempo)>0.:
463 nest.Simulate(simtime-tempo)
464
465
466 endsimulate = time.time()
467 f.close()
468 else:
469 if len(InfoProtocol):
470 tempo=0
471 for contatore in range(0,len(InfoProtocol)):
472 nest.Simulate(InfoProtocol[contatore][0])
473 temp=InfoProtocol[contatore][0]
474 tempo=tempo+temp
475 if InfoProtocol[contatore][2]==4:
476 nest.SetStatus(NoisePop[InfoProtocol[contatore][1]],params={"rate": float(InfoBuild[1+InfoProtocol[contatore][1]][2]*InfoProtocol[contatore][3])})
477 #print "Population:", InfoProtocol[contatore][1] ,";Parameter:", InfoProtocol[contatore][2] ,"; Value: ",InfoProtocol[contatore][3]
478 if InfoProtocol[contatore][2]==12:
479 nest.SetStatus(NeuronPop[InfoProtocol[contatore][1]], params={"b": float(InfoProtocol[contatore][3])})
480 #print "Population:", InfoProtocol[contatore][1] ,";Parameter:", InfoProtocol[contatore][2] ,"; Value: ",InfoProtocol[contatore][3]
481
482 else:
483 nest.Simulate(simtime)
484 tempo=simtime
485 if (simtime-tempo)>0.:
486 nest.Simulate(simtime-tempo)
487 endsimulate = time.time()
488
489
490 ###################################################################################################################################################################
491
492 #############################------------------------------------------------------------------------
493 #print some information from the simulation
494 #############################------------------------------------------------------------------------
495
496 num_synapses = nest.GetDefaults('static_synapse_hpc')["num_connections"]
497 build_time = endbuild - startbuild
498 connect_time = endconnect - startconnect
499 sim_time = endsimulate - endconnect
500
501 N_neurons=0
502 for i in range(0,int(InfoBuild[0])):
503 N_neurons=N_neurons+int(InfoBuild[i+1][0])
504
505 print(" Network simulation (Python) neuron type:",InfoPerseo[0])
506 print("Number of neurons : {0}".format(N_neurons))
507 print("Number of synapses: {0}".format(num_synapses))
508 print("Building time : %.2f s" % build_time)
509 print("Connecting time : %.2f s" % connect_time)
510 print("Simulation time : %.2f s" % sim_time)
511
512 Fine=time.time()
513 print ("Total Simulation time : %.2f s" % (Fine-Inizio))
514 {{/code}}
515
516 === ===
517
518 === ===
519
520 === [[image:image-20220127171242-1.png]] ===
521
522 === Results: ===
523
524 [[the output of this simulation is...>>https://drive.ebrains.eu/smart-link/215f8213-17e3-468b-b573-e6eaf49d315e/]]
525
526
527
528
529
530
531 ==== ====