First let me make sure I know what you want to do.
sgratiy wrote:I need to save a vector of CSDs at every time step which results from summign currents for entire population of cells.
Presumably you are calculating current source density by a procedure roughly equivalent to this:
define a set of J volumes vol_j
within the jth volume there are K_j neurites
the kth neurite in volume j generates membrane current i_jk
the current source density in volume j is
In the course of a simulation, each i_jk varies with time so you need to calculate CSD_j after each fadvance(), and so you will want to capture the time course of CSD_j to a vector.
To handle all compartments, you will need J CSD vectors--one for each volume--if your simulation is executed serially. However, if your model is parallelized with MPI, the cells will be distributed over nhost processors. This means that each host will have to capture its own set of J CSD vectors in the course of a simulation. The total number of CSD vectors for the entire model will be J*nhost. After the simulation ends, you will have to "reduce" these to a single set of J CSD vectors (the question of how to do this will be addressed in a subsequent post).
Have I guessed correctly?
For serial implentation I used to modify the following procedure to update my vector of CSDs
I'd do something similar except I'd get rid of localobj tobj and do
Code: Select all
proc advance() {
fadvance()
update_csd()
dlist.append{csdvec.c)
}
or maybe write update_csd as an obfunc that returns the CSD vector and change its name to new_csd(), so proc advance becomes
Code: Select all
proc advance() {
fadvance()
dlist.append{new_csd())
}
but in parallel simulation, as far as I understand, this will not work as simulation control does not have proc advance().
No, proc advance() still exists and can be used to insert calculations that are to be performed after the model's states have been updated. However, on each host you can only add up the currents that are generated by the neurites that exist on that host. This means you must use gids to identify the cells, and the procedure that calculates the CSD must use
if (pc.gid_exists(gid))
statements to make sure that each host adds up only the currents that are produced by neurites that are actually on that host.
For the local variable timestep method . . .
Why bring this up? Are you using local variable time step integration, or fixed time steps (which are inherently global)? If you're using fixed time steps, the caveats about local variable timestep are not relevant. If you are using local variable time step, my first suggestion is stop--use fixed time steps. Local variable time steps will make data reduction inordinately complex, and for what possible benefit?
Two further comments about Vector.record before I give you the chance to reply--
Vector.record can capture the time course of any scalar variable, not just a range variable or some other variable that "belongs" to a section. There is only one small catch: if you try to capture the time course of t itself, e.g. with
objref tvec
tvec = new Vector()
tvec.record(&t)
at least one section must exist before the tvec.record(&t) statement is executed; if not, hoc will produce a "Section access unspecified" error message when it encounters the tvec.record(&t) statement. To make this error message go away, it is sufficient to create a "dummy" section like so
create foo
and then you can execute tvec.record(&t) with impunity. And it will work correctly, too.