Movie Run for a user defined array

[sgratiy]

I need to be able to plot an array whose values change as a function of time during the simulation using Movie Run. In the code below an array in question 'csdvec' is a linear transformation of i_membrane(x). Currently, I run the simulation for a certain duration of time and then reexecute the code listed below, but this is far worse than having a nice time dynamics offered by Movie Run. How can I see the time course of my array on the graph?

Code: Select all

double csdvec[Nzpos+1] 			// vector of CSDs

	i = 0	
	forall {
	    	for (l,0) {
			i = i + 1 				// global segment index
			for fi = 0, Nzpos { csdvec[fi] = csdvec[fi] + wmat[fi][i]*i_membrane(l)*area(l) }		
			print i, i_membrane(l)
		}	
	}

objref g	
g = new Graph()		
g.size(-120, 850, 0, 10)
for fi=0,Nzpos {	
	x = fi*dy-120
	y = csdvec[fi] 	
	g.line(x, y)	
}
g.flush()

[ted]

How can I see the time course of my array on the graph

Instead of doubles and Graph.line, I'd use Vectors and the Vector class's plot method. And I'd let the standard run system automate the graph updates.

Code: Select all

objref g, csdvec, xvec
csdvec = new Vector(Nzpos+1)
xvec = csdvec.c
xvec.indgen()
xvec.mul(dy).sub(120)
csdvec.plot(g, xvec)

proc advance() { // custom proc advance()--see p. 161 in NEURON Book
  fadvance()
  update_csdvec()
}

proc update_csdvec() {
   . . . // statements to update values in csdvec
  g.flush() // redraws plot of csdvec in g
}

I'm probably misunderstanding something about how you're calculating the contents of csdvec

Code: Select all

	i = 0	
	forall {
	    	for (l,0) {
			i = i + 1 				// global segment index
			for fi = 0, Nzpos { csdvec[fi] = csdvec[fi] + wmat[fi][i]*i_membrane(l)*area(l) }		
			print i, i_membrane(l)
		}	
	}

but if wmat is two dimensional, shouldn't csdvec also be two dimensional? In which case what is the meaning of a two dimensional plot of csdvec, unless you mean to do a two dimensional plot of each "row" or "column" of csdvec.

[sgratiy]

Ted,

thanks it worked. However, I have some issues with saving this window into a session file for reuse. The first time I execute my code with a custom advance() procedure I created the graph using 'g = new Graph()' statement. Then, I saved my session which included a plot of 'csdvec' into a .ses file and then added that session file at the end of my init.hoc file. The next time I executed init.hoc, I got an extra window. The reason is that the .ses file saved the 'csdvec' plot window with 'save_window_' object reference rather than 'g' object reference given in my .hoc file. Therefore, upon execution of the init.hoc, I got two references: 'g' and 'save_window' pointing to window objects. To fix that I had to include a statement ' g=save_window_' and remove the statement 'g =new Graph()' to prevent interpreter from creating an extra window. This does not look like a good way of programming to me, and I wonder if there is a more elegant solution. The problem is that all object references to windows in .ses file are called 'save_window_' therefore it is easy to point graph to a wrong window. On top of all that, I was not able to save 'Color/Brush' properties for this window. I would change the line color, but upon reload of the session, the line color turns gray.

Given these issues, is possible plot a user defined vector via GUI without writing a custom advance() procedure, and experiencing the above difficulties?

The code of the file which does window flush is given below:

Code: Select all

objref g, csdvec, xvec
csdvec = new Vector(Nzpos+1)
xvec = csdvec.c
xvec.indgen()
xvec.mul(dy)
//g = new Graph()		
g = save_window_


csdvec.plot(g, xvec)
double csddouble[Nzpos+1]

proc update_csdvec() { local j, fi

double csddouble[Nzpos+1] 			// double array of CSDs

	j = 0	
	forall {
	    	for (l,0) {
			j = j + 1 				// global segment index
			for fi = 0, Nzpos { csddouble[fi] = csddouble[fi] + wmat[fi][j]*i_membrane(l)*area(l) }		
		}	
	}

	for fi = 0, Nzpos { csdvec.x[fi] = csddouble[fi]} // pack csd array of doubles into a vector		

  	g.flush() // redraws plot of csdvec in g

}

proc advance() { // custom proc advance()--see p. 161 in NEURON Book
  fadvance()
  update_csdvec()
}

Thanks

[ted]

I have some issues with saving this window into a session file for reuse.

So don't reuse the session file. Use hoc statements to create the Graph. This gives you complete control over the name of the objref that refers to the Graph. It also gives you complete control over the location, size, axis scaling, and line colors used in the Graph.

I would be tempted to put the code that creates the Graph, plus the code for the custom proc advance(), into a hoc file with a special name something like csd.hoc. Then I would put a load_file("csd.hoc") statement in the "instrumentation" section of my main program file. If I didn't like the location on the screen where the Graph is drawn, I'd move it to where I want it, then use the Print & File Window Manager to save the Graph to a session file all by itself, and finally I'd read the ses file to discover the correct screen coordinates to use in the .view() statement in csd.hoc. Likewise, if I needed to change the size, axis scaling, or line color, I'd run the program, make the changes to the Graph, save to a ses file, and examine that file to discover the Graph method(s) to call, and the parameter values to use, in csd.hoc.

got two references: 'g' and 'save_window' pointing to window objects. To fix that I had to include a statement ' g=save_window_' and remove the statement 'g =new Graph()' to prevent interpreter from creating an extra window. This does not look like a good way of programming to me

so don't do it. Each tool has its most appropriate use. ses files are perfect for saving and recreating about 90% of the graphs that I use, but I find that the remaining 10% require customizations. For that 10% I may start with a ses file, but then I steal reusable code from the ses file to make my customizations.

The problem is that all object references to windows in .ses file are called 'save_window_'

because it's computer-generated code. The first two things I do to the code I steal from a ses file is
(1) ignore the ses file's "boilerplate" (stuff that NEURON puts there to help with the automatic recreation of objects saved from the GUI), and
(2) identify those statements that I need to reuse, and replace the computer-generated "generic" names in those statements (like save_window_) with the variable or objref names that I want to use. For example, see
How to use hoc to plot a variable vs. time
http://www.neuron.yale.edu/phpBB/viewto ... f=30&t=552

On top of all that, I was not able to save 'Color/Brush' properties for this window. I would change the line color, but upon reload of the session, the line color turns gray.

So after you change the color, save the Graph to a ses file, examine the statements in that file, and discover* the statement that is responsible for making the line red. Then go back to your hoc code and change it accordingly.

*--of course, the statement is documented somewhere in the Programmer's Reference, but often it is much easier to change a Graph and see what happens to the ses file. That tells me the name of the method, and which argument in the method call, and makes it easier for me to read the documentation.

is possible plot a user defined vector via GUI without writing a custom advance() procedure

"All things are possible through programming" but some strategies and implementations are more idiomatic than others. Let me know if you have questions about anything I have written here.

Code: Select all

double csddouble[Nzpos+1] 			// double array of CSDs
 . . .
			for fi = 0, Nzpos { csddouble[fi] = csddouble[fi] + wmat[fi][j]*i_membrane(l)*area(l) 
 . . .
	for fi = 0, Nzpos { csdvec.x[fi] = csddouble[fi]} // pack csd array of doubles into a vector	

I'm still not sure what the statement
for fi = 0, Nzpos { csddouble[fi] = csddouble[fi] + wmat[fi][j]*i_membrane(l)*area(l)
is doing, and what makes this multiply nested loop an improvement over the use of transfer resistances (are you generating a 2 dimensional map of extracellular potentials? doesn't seem likely, because csddouble is one dimensional), but presumably you know what you're doing.

That aside, I'd expect this code to execute more quickly if you did everything with Vectors (and maybe a Matrix, or even just an array of Vectors, for wmat). That would eliminate the "packing" step, and the Vector multiplication would occur at compiled code speed.

[sgratiy]

Ted,

I followed your suggestion and implemented graphs via hoc statements. Also, thanks for pointing out the Matrix class for me. I have now implemented my code using strictly Vector and Matrix classes. However, I have a minor problem: The description of Graph class states that when using .plot function it is only necessary to call g.flush() to display further changes to the vector. My code refuses to do that and simply does not display the updated vector, but it would update the vector when I call 'plot' function. According to the programmers reference, calling g.flush() is more efficient and it would be nice to have this working. Please see my code below. When I remove comment before the statement 'csdvec.plot(g,xvec,2,2)' the graph is updated as needed.

Code: Select all

objref g, csdvec, xvec, i_segment
csdvec = new Vector(Nzpos+1)
i_segment = new Vector(Nseg_all+1)
xvec = csdvec.c
xvec.indgen()
xvec.mul(dy)

g = new Graph(0)
g.size(0,1100,-0.63,0.06)
g.view(0, -0.63, 1100, 0.69, 1022, 371, 525.12, 200.32)
g.label(0.658147, 0.891374, "CSD", 2, 1, 0, 0, 3)

csdvec.plot(g, xvec,3,2)
proc update_csd() { local j, fi

	j = 0
	forall {
	    	for (l,0) {
			j = j + 1 				// global segment index
			i_segment.x[j] = i_membrane(l)*area(l)
		}	
	}
	csdvec = wmat.mulv(i_segment) // get CSD
	csdvec.div(dy)		  // remove dependence on descritization		

//	csdvec.plot(g, xvec,2,2)
}

proc advance() { // custom proc advance()--see p. 161 in NEURON Book
  fadvance()
  update_csd()
  g.flush()

}  

Now about calculation of csdvec: it is assumed that cortical column is homogeneous in the radial direction (which is true to some extent) and therefore CSD is only a function of the cortical depth 'z'. Therefore, CSD(z) can be found from i_membrane(z) passing though all branches of the cell at the depth 'z'. This relationship is described by my 'wmat' matrix which maps currents from each segment to the CSD at a given cortical depth. In the previous implementation of this mapping the first two external loops were implemented to loop over all segments of the cell in order to find their contribution to CSD at a certain depth. But now I create the Vector 'i_segment' (total membrane current through each segment) and multiply it by 'wmat'. I hope this explains it.
I do not know if this approach is better or worse than using transfer resistances. In my way, I can calculate CSD/LFPs with an arbitrarily fine vertical scale needed for my calculations. I do not know if it is easy and computationally efficient to calculate transfer resistances. I would like to hear your opinion on that.

Thanks.

[ted]

I have now implemented my code using strictly Vector and Matrix classes.

That should help to speed things up.

However, I have a minor problem: The description of Graph class states that when using .plot function it is only necessary to call g.flush() to display further changes to the vector. My code refuses to do that and simply does not display the updated vector, but it would update the vector when I call 'plot' function. According to the programmers reference, calling g.flush() is more efficient and it would be nice to have this working.

Here's the problem: csdvec is not a Vector--it is merely an objref, that is, a pointer to a Vector. This statement
wmat.mulv(i_segment)
very efficiently calculates the product of wmat and i_segment, and returns a pointer to a new Vector. Consequently, the statement
csdvec = wmat.mulv(i_segment)
breaks the link between csdvec and whatever Vector it previously referenced--but that previously referenced Vector was the one whose contents were plotted in the Graph. That's why a new csdvec.plot(g, xvec . . . ) statement is necessary.

Thanks for explaining the calculation of csdvec--would you call it a "mean field" approach? At some level this is mathematically equivalent to the use of transfer resistances, but for fields generated by a large number of cells that are more or less synchronous (or which can be lumped into synchronous subsets) it leads to a computation that is certainly more efficient.

[sgratiy]

Yes, this is a mean field calculation, because it assumes that we calculate the CSD from an entire population of identical cells distributed along a cortical layer with some density, rather than from a single cell.

Thanks,
Sergey