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