Attention: The EBRAINS IDM/IAM will be down tomorrow, Wednesday 17nd December, from 17:00 CET for up to 30 minutes for maintenance. Please be aware that this will affect all services that require login or authentication.


Wiki source code of Code description

Version 15.1 by galluzziandrea on 2022/01/27 17:01
Hide last authors
galluzziandrea 5.1 1 == Introduction (path and modules): ==
galluzziandrea 1.1 2
galluzziandrea 3.1 3 First of all we check the path and import the necessary modules .
galluzziandrea 1.1 4
galluzziandrea 5.1 5 === Check where I am and place myself in the right folder: ===
galluzziandrea 3.1 6
galluzziandrea 2.2 7 {{code language="python"}}
galluzziandrea 3.1 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()))
galluzziandrea 2.2 19 {{/code}}
galluzziandrea 1.1 20
galluzziandrea 5.1 21 === Import the modules necessary for the simulation: ===
galluzziandrea 1.1 22
23 {{code language="python"}}
galluzziandrea 3.1 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)
galluzziandrea 1.1 33 {{/code}}
34
galluzziandrea 4.1 35 === ===
galluzziandrea 2.2 36
galluzziandrea 5.1 37 === Define necessary classes to import the Initialization Files: ===
galluzziandrea 2.2 38
galluzziandrea 7.1 39 {{code language="python" title=" "}}
galluzziandrea 5.1 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
galluzziandrea 8.1 238 === Defining general and nest.kernel parameters ===
galluzziandrea 5.1 239
galluzziandrea 8.1 240 {{code language="python"}}
241 #############################------------------------------------------------------------------------
242 #Clean the Network
243 #############################------------------------------------------------------------------------
244 nest.ResetKernel()
galluzziandrea 5.1 245
galluzziandrea 8.1 246 #############################------------------------------------------------------------------------
247 #insert the introductory parameters of the simulation
248 #############################------------------------------------------------------------------------
galluzziandrea 5.1 249
galluzziandrea 8.1 250
251 dt = 0.1 # the resolution in ms
252 StartMisure=0. # start time of measurements
253 simtime = int(float(InfoPerseo[3])) # Simulation time in ms (200 s)
254 if simtime<=StartMisure: # If the simulation time is less than StartMisure, it is increased by StartMisure
255 simtime=simtime+StartMisure
256 start=0.0 # start time of poissonian processes
257 origin=0.0 # temporal origin
258
259 #############################------------------------------------------------------------------------
260 # Kernel parameters
261 #############################------------------------------------------------------------------------
262 LNT=multiprocessing.cpu_count();
263 nest.SetKernelStatus({"local_num_threads": LNT})
264 nest.SetKernelStatus({"resolution": dt, "print_time": True,
265 "overwrite_files": True})
266
267 #############################------------------------------------------------------------------------
268 #"randomize" the seeds of the random generators
269 #############################------------------------------------------------------------------------
270
271 #msd = int(math.fabs(time.process_time()*1000))
272 #N_vp = nest.GetKernelStatus(['total_num_virtual_procs'])[0]
273 #pyrngs = [numpy.random.RandomState(s) for s in range(msd, msd+N_vp)]
274 #nest.SetKernelStatus({"grng_seed" : msd+N_vp})
275 #nest.SetKernelStatus({"rng_seeds" : list(range(msd+N_vp+1, msd+2*N_vp+1))})
276 {{/code}}
277
278 === Building the network: neuronal populations , Poisson processes and spike detectors ===
279
280 {{code language="python"}}
281 #############################------------------------------------------------------------------------
282 print("Building network")
283 #############################------------------------------------------------------------------------
284
285 startbuild = time.time() #initialize the calculation of the time used to simulate
286
287 NeuronPop=[]
288 NoisePop=[]
289 DetectorPop=[]
290
291 #define and initialize the populations of neurons with the parameters extracted from the.ini files
292 for i in range(1,int(InfoBuild[0])+1):
293 if int(InfoBuild[i][7])==0:
294 app=float(InfoBuild[i][5])
295 else:
296 app=0.
297 app2= nest.Create("aeif_psc_exp", int(InfoBuild[i][0]),params={"C_m": 1.0,
298 "g_L": 1.0/float(InfoBuild[i][3]),
299 "t_ref": float(InfoBuild[i][6]),
300 "E_L": 0.0,
301 "V_reset": float(InfoBuild[i][5]),
302 "V_m": app,
303 "V_th": float(InfoBuild[i][4]),
304 "Delta_T": 0.,
305 "tau_syn_ex": 1.0,
306 "tau_syn_in": 1.0,
307 "a": 0.0,
308 "b": float(InfoBuild[i][10]),
309 "tau_w": float(InfoBuild[i][9]),
310 "V_peak":float(InfoBuild[i][4])+10.0})
311 NeuronPop.append(app2)
312
313 #define and initialize the poisson generators and the spike detectors with the parameters extracted from the.ini files
314
315 for i in range(1,int(InfoBuild[0])+1):
316 app3= nest.Create("poisson_generator",params={"rate": float(InfoBuild[i][1]*InfoBuild[i][2]),
317 'origin':0.,
318 'start':start})
319 NoisePop.append(app3)
320 app4 = nest.Create("spike_recorder",params={ "start":StartMisure})
321 DetectorPop.append(app4)
322
323 endbuild = time.time()
324 {{/code}}
325
galluzziandrea 13.1 326 (% class="wikigeneratedid" %)
327 === [[image:image-20220127165908-2.png||height="659" width="1149"]] ===
328
galluzziandrea 9.1 329 === Connecting the network nodes: neuronal populations, Poisson processes and spike detectors ===
330
331 {{code language="python"}}
332 #############################------------------------------------------------------------------------
333 print("Connecting ")
334 #############################------------------------------------------------------------------------
335
336 startconnect = time.time()
337 Connessioni=[]
338 Medie=[]
339
340 #create and define the connections between the populations of neurons and the poisson generators
341 #and between the populations of neurons and the spike detectors with the parameters extracted from the.ini files
342
343 for i in range(0,int(InfoBuild[0])):
344 nest.Connect(NoisePop[i], NeuronPop[i], syn_spec={'synapse_model': 'static_synapse_hpc',
345 'delay': dt,
346 'weight': nest.math.redraw(nest.random.normal(mean=float(InfoConnectNoise[i+1][0]),
347 std=(float(InfoConnectNoise[i+1][1])*float(InfoConnectNoise[i+1][0]))),
348 min=0., max=float('Inf'))
349 })
350 nest.Connect(NeuronPop[i][:int(InfoBuild[i+1][0])], DetectorPop[i], syn_spec={"weight": 1.0, "delay": dt})
351
352 #create and define the connections between the populations of neurons with the parameters extracted from the.ini files
353
354 for i in range(0,len(InfoConnectPop[1:])):
355
356 conn=nest.Connect(NeuronPop[int(InfoConnectPop[i+1][1])], NeuronPop[int(InfoConnectPop[i+1][0])],
357 {'rule': 'pairwise_bernoulli',
358 'p':float(InfoConnectPop[i+1][2]) },
359 syn_spec={'synapse_model': 'static_synapse_hpc',
360 'delay':nest.math.redraw(nest.random.exponential(beta=float(1./(2.99573227355/(float(InfoConnectPop[i+1][4])-float(InfoConnectPop[i+1][3]))))),
361 min= numpy.max([dt,float(1./float(InfoConnectPop[i+1][4]))]),
362 max= float(1./(float(InfoConnectPop[i+1][3])-dt/2))),
363
364 'weight':nest.random.normal(mean=float(InfoConnectPop[i+1][6]),
365 std=math.fabs(float(InfoConnectPop[i+1][6])*float(InfoConnectPop[i+1][7])))})
366
367
368 endconnect = time.time()
369 {{/code}}
370
371 === ===
372
373 === ===
374
375 === Simulating: neuronal time evolution. ===
376
377 === ===
378
379 {{code language="python"}}
380 #############################------------------------------------------------------------------------
381 print("Simulating")
382 #############################------------------------------------------------------------------------
383 ###################################################################################################################################################################
384 if Salva:
385 print("I m going to save the data")
386 #x=str(iterazioni)
387 f = open(FileName,"w")
388 if len(InfoProtocol):
389 print("I m going to split the simulation")
390 tempo=0
391 for contatore in range(0,len(InfoProtocol)):
392 appoggio1=int((tempo+InfoProtocol[contatore][0])/1000.)
393 appoggio2=int(tempo/1000.)
394 appoggio3=tempo+InfoProtocol[contatore][0]
395 if (appoggio1-appoggio2)>=1:
396 T1=(1+appoggio2)*1000-tempo
397 nest.Simulate(T1)
398 #Save the Data!!!!
399 ###########################################################
400 Equilibri=[]
401 for i in range(0,int(InfoBuild[0])):
402 Equilibri.append([])
403 a=nest.GetStatus(DetectorPop[i])[0]["events"]["times"]
404 if len(a)>0:
405 Trange=(1000*int(numpy.min(a)/1000.),1000*int(numpy.min(a)/1000.)+1000)
406 hist,Tbin=numpy.histogram(a,200,(Trange[0],Trange[1]))
407 Equilibri[i]=hist*1000./(5.*int(InfoBuild[i+1][0]))
408 else:
409 Trange=(1000*int(tempo/1000.),1000*int(tempo/1000.)+1000)
410 hist=numpy.zeros(200)
411 Tbin=numpy.linspace(Trange[0],Trange[1],num=201)
412 Equilibri[i]=hist
413 nest.SetStatus(DetectorPop[i],{'n_events':0})
414 for j in range(0,len(hist)):
415 f.write(str(Tbin[j])+" ")
416 for i in range(0,int(InfoBuild[0])):
417 f.write(str(Equilibri[i][j])+" ")
418 f.write("\n ")
419 ###########################################################
420 tempo=tempo+T1
421 for contatore2 in range(1,(appoggio1-appoggio2)):
422 nest.Simulate(1000.)
423 #Save the Data!!!!
424 ###########################################################
425 Equilibri=[]
426 for i in range(0,int(InfoBuild[0])):
427 Equilibri.append([])
428 a=nest.GetStatus(DetectorPop[i])[0]["events"]["times"]
429 if len(a)>0:
430 Trange=(1000*int(numpy.min(a)/1000.),1000*int(numpy.min(a)/1000.)+1000)
431 hist,Tbin=numpy.histogram(a,200,(Trange[0],Trange[1]))
432 Equilibri[i]=hist*1000./(5.*int(InfoBuild[i+1][0]))
433 else:
434 Trange=(1000*int(tempo/1000.),1000*int(tempo/1000.)+1000)
435 hist=numpy.zeros(200)
436 Tbin=numpy.linspace(Trange[0],Trange[1],num=201)
437 Equilibri[i]=hist
438 nest.SetStatus(DetectorPop[i],{'n_events':0})
439 for j in range(0,len(hist)):
440 f.write(str(Tbin[j])+" ")
441 for i in range(0,int(InfoBuild[0])):
442 f.write(str(Equilibri[i][j])+" ")
443 f.write("\n ")
444 tempo=tempo+1000.
445 T2=appoggio3-tempo
446 nest.Simulate(T2);
447 tempo=tempo+T2;
448 else:
449 nest.Simulate(InfoProtocol[contatore][0])
450 temp=InfoProtocol[contatore][0]
451 tempo=tempo+temp
452 if InfoProtocol[contatore][2]==4:
453 nest.SetStatus(NoisePop[InfoProtocol[contatore][1]],params={"rate": float(InfoBuild[1+InfoProtocol[contatore][1]][2]*InfoProtocol[contatore][3])})
454 if InfoProtocol[contatore][2]==12:
455 nest.SetStatus(NeuronPop[InfoProtocol[contatore][1]], params={"b": float(InfoProtocol[contatore][3])})
456 else:
457 nest.Simulate(simtime)
458 tempo=simtime
459 if (simtime-tempo)>0.:
460 nest.Simulate(simtime-tempo)
461
462
463 endsimulate = time.time()
464 f.close()
465 else:
466 if len(InfoProtocol):
467 tempo=0
468 for contatore in range(0,len(InfoProtocol)):
469 nest.Simulate(InfoProtocol[contatore][0])
470 temp=InfoProtocol[contatore][0]
471 tempo=tempo+temp
472 if InfoProtocol[contatore][2]==4:
473 nest.SetStatus(NoisePop[InfoProtocol[contatore][1]],params={"rate": float(InfoBuild[1+InfoProtocol[contatore][1]][2]*InfoProtocol[contatore][3])})
474 #print "Population:", InfoProtocol[contatore][1] ,";Parameter:", InfoProtocol[contatore][2] ,"; Value: ",InfoProtocol[contatore][3]
475 if InfoProtocol[contatore][2]==12:
476 nest.SetStatus(NeuronPop[InfoProtocol[contatore][1]], params={"b": float(InfoProtocol[contatore][3])})
477 #print "Population:", InfoProtocol[contatore][1] ,";Parameter:", InfoProtocol[contatore][2] ,"; Value: ",InfoProtocol[contatore][3]
478
479 else:
480 nest.Simulate(simtime)
481 tempo=simtime
482 if (simtime-tempo)>0.:
483 nest.Simulate(simtime-tempo)
484 endsimulate = time.time()
485
486
487 ###################################################################################################################################################################
488
489 #############################------------------------------------------------------------------------
490 #print some information from the simulation
491 #############################------------------------------------------------------------------------
492
493 num_synapses = nest.GetDefaults('static_synapse_hpc')["num_connections"]
494 build_time = endbuild - startbuild
495 connect_time = endconnect - startconnect
496 sim_time = endsimulate - endconnect
497
498 N_neurons=0
499 for i in range(0,int(InfoBuild[0])):
500 N_neurons=N_neurons+int(InfoBuild[i+1][0])
501
502 print(" Network simulation (Python) neuron type:",InfoPerseo[0])
503 print("Number of neurons : {0}".format(N_neurons))
504 print("Number of synapses: {0}".format(num_synapses))
505 print("Building time : %.2f s" % build_time)
506 print("Connecting time : %.2f s" % connect_time)
507 print("Simulation time : %.2f s" % sim_time)
508
509 Fine=time.time()
510 print ("Total Simulation time : %.2f s" % (Fine-Inizio))
511 {{/code}}
512
513 === ===
514
galluzziandrea 9.2 515 === Results: ===
galluzziandrea 5.1 516
galluzziandrea 9.2 517 the output of this simulationo is...
galluzziandrea 9.1 518
519
520
521
522
galluzziandrea 9.2 523
galluzziandrea 4.1 524 ==== ====