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 ==== ====