Using For Loops to Call Multiple Simulation Runs

Anything that doesn't fit elsewhere.
Post Reply
Yaeger
Posts: 33
Joined: Mon Aug 19, 2013 4:36 pm

Using For Loops to Call Multiple Simulation Runs

Post by Yaeger »

Hi, I am trying to write a hoc file that will iterate through a user defined number of runs but I keep getting errors that Neuron says are at the end of the hoc file but which do not make much sense. The reason for iterating through the number of runs is that I am changing a list of event times on each run. For instance:

/Applications/NEURON-7.3/nrn/x86_64/bin/nrniv: syntax error
in mosinit2.hoc near line 65
file = new File()

If I move that command to earlier in the file, I get a syntax error pointing to the last line of the hoc file again.

My protocol for populating the vector is contained in a hoc file that is loaded by the mosinit hoc file that I use to run the simulation. This protocol is called interval:

Code: Select all

// create vector for spike times
objref svec
svec = new Vector()
proc interval () { local i, x
	low = $1
	reps = $3
	x = low
	dx = $2
	for i = 1,reps {
		svec.append(x)
		x = x + dx
		}
}


// create a vector stim that generates events at times in svec
objref vecstim
vecstim = new VecStim()
vecstim.play(svec)
Then in the mosinit file I load the file containing the above code and the rest of the simulation parameters and then:

Code: Select all

dt = 0.025
tstop = 10000
v_init = -60

proc init() {
	finitialize(-60)
	frecord_init()
}

proc integrate() {
	while (t<tstop) {
		fadvance()
	}
}

proc go() {
	init()
	integrate()
}

proc newvec() {
	objref grcspk, nc, nil
	grcspk = new Vector()
	grc nc = new NetCon(&v(0.5), nil)
	nc.threshold = 0
	nc.record(grcspk)
}
	
proc counter() { local firstspk, lastspk, totspk, spkfreq
	totspk = grcspk.size
	firstspk = grcspk.x[0]
	lastspk = grcspk.x[totspk-1]
	if (totspk > 1) {
		spkfreq = (totspk - 1)/(lastspk - firstspk)
	} else{
		spkfreq = 0
}
	spkcnts.append(spkfreq)
}

proc gocount() { local ii
	for ii = 1,2 {
	newvec()
	interval(3000, ii, 100)
	go()
	counter()
}

gocount()
newvec is to create a new spike-counting vector each run of the simulation.
counter finds the average spike frequency and appends this to a vector.
gocount launches the simulation.

Neuron does not seem to be happy with this code but it is not giving me error messages that make much sense. How can I run the simulation with iteratively calls so that I can increment my call to the interval procedure?
Yaeger
Posts: 33
Joined: Mon Aug 19, 2013 4:36 pm

Re: Using For Loops to Call Multiple Simulation Runs

Post by Yaeger »

I found part of the mistake. I was missing a second } in the gocount procedure:

Code: Select all

proc gocount() { local ii
	for ii = 1,2 {
	newvec()
	interval(3000, ii, 100)
	go()
	counter()
	objref svec, vecstim
	}
}
Yaeger
Posts: 33
Joined: Mon Aug 19, 2013 4:36 pm

Re: Using For Loops to Call Multiple Simulation Runs

Post by Yaeger »

Now I get the error message:

Code: Select all

net_send td-t = -396 SelfEvent target=VecStim[0] 3000 flag=1
/Users/danielyaeger/Desktop/Dan_Feedback_Network/x86_64/special: line 13:  6385 Abort trap: 6           "${NRNIV}" -dll "/Users/danielyaeger/Desktop/Dan_Feedback_Network/x86_64/.libs/libnrnmech.so" "$@"
nrngui exit status was 134
So for some reason the VecStim is saying that the first event occurred in negative time (the first event occurs at 3000)?
Yaeger
Posts: 33
Joined: Mon Aug 19, 2013 4:36 pm

Re: Using For Loops to Call Multiple Simulation Runs

Post by Yaeger »

Alternatively, does anyone have a better strategy for how to use a hoc file to go through mutiple simulation runs in which a parameter is changed using a loop?
Yaeger
Posts: 33
Joined: Mon Aug 19, 2013 4:36 pm

Re: Using For Loops to Call Multiple Simulation Runs

Post by Yaeger »

Ok, now I almost have something that works, but I need to close a hoc file in my code. How do I close selected hoc files? There is a load_file("file.hoc") command but I cannot find the equivalent close file command for hoc files.
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Using For Loops to Call Multiple Simulation Runs

Post by ted »

Yaeger wrote:does anyone have a better strategy for how to use a hoc file to go through mutiple simulation runs in which a parameter is changed using a loop?
Sorry this reply is so delayed. Much of the following may now be moot, but perhaps there are some items that will be useful in the future.
The original error message

Code: Select all

syntax error
in mosinit2.hoc near line 65
file = new File()
suggests that "file" was not declared to be an objref.
I am trying to write a hoc file that will iterate through a user defined number of runs
. . .
changing a list of event times on each run
The code you present looks workable, with one exception: proc interval() will fail to do what you want, because each time it is called it merely tacks a new sequence of values onto the end of svec. This is probably what caused the "net_send td-t = -396" error message.

Also, the code misses a couple of opportunities to be shorter and simpler--shorter code offers fewer opportunities for programmer-caused errors, and is easier to maintain. For one thing, the Vector class has methods that make short work of filling a vector with numerical values. For another, it is rarely necessary to write one's own proc init() or "proc integrate()" or a proc go() that calls init or "integrate"--just use NEURON's standard run system, which you get by

Code: Select all

load_file("nrngui.hoc")
or, if you don't want to see the NEURON Main Menu toolbar,

Code: Select all

load_file("stdgui.hoc")
at the start of your main hoc file. Don't worry about the "gui" in these file names--your program will execute just fine, even if NEURON was compiled without Interviews.

The revised procedures (these eliminate several of the procs in your examples):

Code: Select all

// $1 start time, $2 DT, $3 number of intervals
proc before_run() { // instead of "interval"
  svec.resize($3+1)
  svec.indgen($2) // 0, DT, 2*DT...
  svec.add($1)
}

proc batrun() { local i
  for i=1,$1 {
    before_run()
    run()
    counter()
  }
}
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Using For Loops to Call Multiple Simulation Runs

Post by ted »

Yaeger wrote:I need to close a hoc file
Actually, no. More likely you were writing to a text file, and you want to close that file. You'll find what you need in the Programmer's Reference--here
http://www.neuron.yale.edu/neuron/stati ... /file.html
if you're using the File class, and here
http://www.neuron.yale.edu/neuron/stati ... html#wopen
if you're using wopen.
Yaeger
Posts: 33
Joined: Mon Aug 19, 2013 4:36 pm

Re: Using For Loops to Call Multiple Simulation Runs

Post by Yaeger »

Thank you Ted. That code ran beautifully, and, as reccomended, I am using the standard run procedure. I have a couple of questions on combining the standard run system with hoc code though. If I write dt = 0.001 in my hoc file that launches the simulation, does that instruct the run system to use this dt as the time step? Also, can you specify the integration mode (CVODE, backward euler, etc) in hoc code when using the standard run system?
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Using For Loops to Call Multiple Simulation Runs

Post by ted »

Yaeger wrote:If I write dt = 0.001 in my hoc file that launches the simulation, does that instruct the run system to use this dt as the time step?
Here's how to answer this for yourself:
make sure that the statement
dt = 0.001
is executed before run() is called, and print the value of dt after run() has been executed. This code will do it:

Code: Select all

proc batrun() { local i, dtold
  dtold = dt
  dt = 0.001
  for i=1,$1 {
    before_run()
    run()
    counter()
  }
  print "executed ", $1, "runs with dt ", dt
  dt = dtold // restore whatever it was before batrun() is called
}
can you specify the integration mode (CVODE, backward euler, etc) in hoc code when using the standard run system?
For fixed step integration the options are backward Euler and two variants of Crank-Nicholson. secondorder specifies which of these will be used.

For adaptive integration, it's essential that each state variable's error tolerance is properly specified. In all but the most trivial of models, the values to specify can only be determined empirically. It is most convenient to use the VariableTimeStep tool.

Code: Select all

1. NEURONMainMenu / VariableStepControl
2. In the VariableTimeStep window click on "Use variable dt" and "Atol Scale Tool"
3. In the Absolute Tolerance Scale Factors panel click on "Analysis Run", then click on "Rescale"
4, Run a simulation and examine the results.
   Inadequate error control will manifest as strange wiggles in plots of v vs. time.
   If you see anything like that
REPEAT
  reduce the global Absolute Tolerance (in the VariableTimeStep window)--try a factor of 10
  do another Analysis Run (Abs Tol Scale Factors panel)
  click on Rescale
  run a simulation
UNTIL glitches no longer are seen in the v vs. time graph
5. Close the Absolute Tolerance Scale Factors panel.
6. Save the VariableTimeStep window to a session file.
   Might as well call the file something very mnemonic, like vts.ses
   If you don't know how to save a selected NEURON GUI tool to a ses file,
   see http://www.neuron.yale.edu/neuron/static/docs/saveses/saveses.html
7. In your hoc code insert this statement
load_file("vts.ses")
   before the statement that calls batrun().
"But I don't want to see any GUI windows when I actually run a batch of simulations."
Don't worry. If you launch NEURON by executing
nrniv whatever.hoc
where whatever.hoc (your hoc file) contains the load_file("stdgui.hoc") statement, you won't see any GUI windows.

"Well, how can I see GUI windows when I am tuning up the VariableTimeStep tool?"
Launch your program by executing
nrngui whatever.hoc
Then bring up a RunControl panel, and a Voltage axis graph, and you can launch a simulation with the RunControl panel's Init & Run button.

"OK, but I'm using dt = 0.001 and tstop is 1e4, and the first stimulus pulse doesn't occur until t = 1000, so I'm sitting here forever waiting for something to happen."
Do you have to wait a long time for your model to reach a steady state? Or is the model spontaneously active, and takes a long time to converge on a limit cycle? You could eliminate a lot of wasted time with a custom initialization. Or you could just take the brute force approach of using the RunControl panel to change tstop to something reasonable, like 100 ms, and bring up a PointProcessManager configured as an IClamp that delivers a reasonable stimulus at t = 10 ms or so. This wouldn't change the dynamical properties of the model's system equations, but it would allow you to execute Analysis Runs in a reasonable length of time.
Post Reply