Exporting vectors to .dat file afte each time step

[sgratiy]

I need to be able to export the spatial CSD vector at each time step to the external file for plotting with Matlab. I implemented the solution, but wanted to ask for the feedback on whether or not this is the optimal way of doing this. The CSD vector is calculated at each time step with the update_csd() procedure which is called by a custom advance() procedure. Then, I store CSD vector at each time step into the matrix column. This matrix is then saved to a file within a custom run() procedure. I have placed the definition of the matrix within the run() to be able to rerun the simulation using the Run Control panel. Obviously, this code is not going to work with the variable time step, because I would not know how many columns to allocate in the matrix. I suppose for the variable time step the only solution is to output each vector directly into the file. However, I was not able to accomplish the latter because I could not figure out how to save vector generated at each time step in the row format. Could you please advice on how to store vectors as rows?

Code: Select all

objref fobj2    	

proc advance() {	 		// custom proc advance() 
	fadvance()
	update_csd() 	 		// calculate CSD vector after each time step
	it=it+1		 		// column index for csdmat	
	csdmat.setcol(it,csdvec) 	// populate it-th column of the output matrix with csdvec(dt*it)
}


proc run() { 				   // custom proc run() 
	stdinit()			   // standard run system procedure
		
	fobj2 = new File("./Data/csd.dat") // redefine the reference to the file for each run
	fobj2.wopen()

	it=0				   // start at the beginning on the matrix for each run	
	Nt = (tstop/dt+1) 		   // number of time steps (= number of columns in csdmat)
	csdmat = new Matrix(Nzpos+1,Nt+1)  // create an empty matrix with enough columns to store the data

	continuerun(tstop)		   // standard run system procedure	

	csdmat.fprint(0,fobj2)		   // output the csdmat to ASCII file	
	fobj2.close()
	it=0 				   // reinitialize, otherwise the Movie Run will not work

}

[ted]

There are a few tricks that would help you solve your problem with a bit more elegance and flexibility.

Instead of changing run(), I'd leave it alone and declare a new proc
proc myrun() {
. . . do whatever needs doing before a run . . .
run()
. . . whatever needs doing after the run . . .
}
That avoids having to call stdinit() or continuerun(). That's one less line of code, which means I have fewer opportunities to make a mistake.
Also, since no data is written to the output file until _after_ completion of the run, I can put all of the file-related statements in a single block. That makes the code simpler, and simpler code is easier to debug. In fact, I think I'd split the file-related stuff into its own proc. Short procs are ease to read and understand. Also, I'd use a localobj to manage the file. That eliminates the need to declare a top-level variable, keeps code a bit cleaner and simpler.
Instead of using a Matrix to hold all of the output Vectors, I'd simply append them to a List. This eliminates the need to know in advance how many time points will be saved.
Here's a draft that implements these ideas:

Code: Select all

objref dlist
dlist = new List()
proc advance() { localobj tobj
  fadvance()
  update_csd()
  tobj=csdvec.c // prevent future changes to csdvec
    // from affecting contents of Vectors that are already in dlist
  dlist.append(tobj)
}
proc save_results() { local i  localobj fobj
  // statements that open output file associated with fobj
  for i=0,dlist.count()-1 { // for one row at a time
    // statements that print the contents of dlist.o(i) in whatever format you like
  }
  // statement that closes output file
}
proc myrun() {
  dlist = new List() // ensure that dlist is empty at start of new run
  run()
  save_results()
}

After I got all this working, I might go on to make one more change if I was thinking about using variable time step integration: instead of using a custom proc advance(), use events to control when the csdvec is generated and saved. This would ensure sampling at regular intervals.

Code: Select all

DT = 3.14159... // or whatever you want for the sampling interval
proc sample_csdvec() { localobj tobj
  update_csd()
  tobj=csdvec.c // prevent future changes to csdvec
    // from affecting contents of Vectors that are already in dlist
  dlist.append(tobj)
  cvode.event("sample_csdvec()", t + DT) // execute sample_csdvec DT ms from now
}
// now we need to call sample_csdvec() at the start of the run
objref fih
fih = new FInitializeHandler("sample_csdvec(DT)") // execute sample_csdvec() at t = 0,
  // right after the mechanism INITIAL blocks have been executed
objref dlist
dlist = new List()
// comment out the custom proc advance
// proc advance() { localobj tobj
//  fadvance()
//  update_csd()
//  tobj=csdvec.c // prevent future changes to csdvec
//    // from affecting contents of Vectors that are already in dlist
//  dlist.append(tobj)
// }
proc save_results() { local i  localobj fobj
  // statements that open output file associated with fobj
  for i=0,dlist.count()-1 { // for one row at a time
    // statements that print the contents of dlist.o(i) in whatever format you like
  }
  // statement that closes output file
}
proc myrun() {
  dlist = new List() // ensure that dlist is empty at start of new run
  run()
  save_results()
}

I think that should work--but typos and stupid or even clever mistakes are always possible.

[sgratiy]

Thanks Ted, I followed your suggestions to use List and the custom myrun() procedure. Very elegant indeed.

[ted]

Elegance is good, but most important is that it works and is easy to implement and maintain.

[sgratiy]

It looks like I was too quick to declare a victory. I am having issues with specifying the proper output format such that it would not stick adjacent elements in the vector together in the ASCII file. I have written the following 'save_results' procedure which prints each raw vector from the list at once.

Code: Select all

proc save_results() { local it localobj fobj2

	fobj2 = new File("./Data/csd.dat") // redefine the reference to the file for each run
	fobj2.wopen()
	for it = 0,dlist.count()-1 {	   // for one row at a time
		dlist.o(it).printf(fobj2,"%11.4e")	
	}	
	fobj2.close()
}

It woks fine for positive numbers, but for negative it does not put space between the adjacent vector elements. I tried to specify a longer width in the format, but this does not fix the problem. So, to circumvent the issue , I written the list of vectors into the matrix and printed the latter into the file.

Code: Select all

proc save_results2() { local it localobj fobj2, csdmat

	fobj2 = new File("./Data/csd.dat") // redefine the reference to the file for each run

	csdmat = new Matrix(csdvec.size(),dlist.count())  // create an empty matrix with enough columns to 
	for it = 0,dlist.count()-1 {	   // for one row at a time
		csdmat.setcol(it,dlist.o(it)) 
	}	
	fobj2.wopen()
	csdmat.fprint(0,fobj2)		   // output the csdmat to ASCII file	
	fobj2.close()
}

I know that this is not the best way in terms of memory handing, but I do not know how to print the list of vectors into the file otherwise. Please advise.

[ted]

Just insert a space or two in the format string. Instead of "%11.4e" you could simply specify "%11.4e "