Data analysis

[erhervert]

Hello again,

I'd like to know if there is a way to customize the output ASCII file obtained from a simulation for further analysis. For example, if I execute my model n times and then I used the print & file window manager to select the V(t) plot, what I get is a file with two columns, one x and one y containing all the results. But, what I'd like to do is to get a file with a single x (time column) and several y's (y0, y1, y2,..., yn), one y from each run. Do you know if it's possible to automate that? Or if there is a way to average all the vectors?

Thanks in advance

[ted]

Set up Vector recording of the variable you are interested in.
Create another Vector that will eventually hold the average over several runs.
For the sake of discussion, let's say the two Vectors are called tmp and avg.

Now automate the process of generating several runs and calculating the average.

Code: Select all

proc batchavg() { local i
  avg = new Vector() // discard whatever is already in avg
  for  i = 1,$1 {
    run()
    if (i==1) {
      avg = tmp.c // copy tmp to avg
    } else {
      avg.add(tmp) // add tmp to avg, element by element
    }
    avg.div($1) // divide each element by number of runs to get the average
  }
}

Example of usage:
batchavg(5)
runs 5 simulations and fills a Vector called avg with the average time course of the recorded variable.

[erhervert]

Ok, I incorporate that into my program and this is what I got:

Code: Select all

///////////////////////////
//* simulation control *//
///////////////////////////
dt = 0.025
tstop = 50
v_init = -65

objref r
r = new Random()
r.normal(.5, .05)

objref avg

proc batchavg() { local i, index
  for i=0, $1 {
    index = r.repick()
    syn.onset = (((10-5)*index)+5)
    run()
  }
  avg = new Vector()
  if (i==1) {
    tmp = avg.c()
    } else {
      avg.add(tmp)
    }
  avg.div($1)
}

Apparently it is working well, but I still having the same two colums (x,y). Do I have to create a second graph to plot the vector avg vs time to see the resulting data? How can I do it? I have been trying to plot the vector as follows, based on the programmers reference example (http://www.neuron.yale.edu/neuron/stati ... .html#plot) but it is not working:

Code: Select all

objref avg, g
g = new Graph()
g.size(0,50,-80,200)
avg.plot(g, 0.025)
graphList[0].append(g)

I'm not sure if it is correct.

Thanks

[ted]

First a comment: if you don't prevent negative values from being assigned to your synapse's start time, you're going to have some runs that produce nonuseful results.

If you want to see the averaged result in a graph, you have to plot it _after_ all of the runs have finished--that is, after batchavg is executed. The graphLists are irrelevant to this. Use the Vector class's plot method. Also, you will probably want the x values to be the times that correspond to the data points, so it is necessary to record t to its own Vector (do that with code next to the statements that set up recording of v or whatever other variable you are recording).

The same considerations apply if you want to see a printout of time and the averaged values.

[erhervert]

First a comment: if you don't prevent negative values from being assigned to your synapse's start time, you're going to have some runs that produce nonuseful results.

Sorry, I did it already using the tsynon() function as you recommend but I posted the uncorrected code (by mistake), this is the good one:

Code: Select all

objref r
r = new Random()
{r.normal(.5, .05)}

func tsynon() {
  ton = -1
  while (ton<0) {
    ton = (((10-5)*r.repick())+5)
  }
  return ton
}

I also did the plot of average vs time using the vector class's plot method. In addition, I incorporate some code to write the results in a txt file at the end of the simulation. But I'm not sure if it is calculating the average correctly, the resulting signal is smaller in size than a single AP and I was expecting something similar to a compound action potential (the sum of the action potentials).

This is the code, could you tell me what do you think?:

Code: Select all

objref avg, time, savdata, tempmatrix, g

time = new Vector()
{time.record(&t)}
avg = new Vector()
{avg.record(&soma.v(0.5))}

proc batchavg() { local i
  for i=0, $1 {
    syn.onset = tsynon()
    run()
  }
  if (i==1) {
    tmp = avg.c()
  } else {
      avg.add(tmp)
    }
  g = new Graph()
  g.size(0,50,-10,5)
  avg.plot(g, time)
  avg.div($1)
  savdata = new File()
  savdata.wopen("results.dat")
  savdata.printf("t avg\n")
  tempmatrix = new Matrix()
  tempmatrix.resize(avg.size(),2)
  tempmatrix.setcol(0, time)
  tempmatrix.setcol(1, avg)
  tempmatrix.fprint(savdata, " %g")
  savdata.close()
}

I'd like to know also if there is a way to make measures on the average trace. I'm interested in things like the peak amplitude, dv/dt and the area under the curve.

Thanks for your help Ted.

[ted]

Looks good. My only suggestion would be to split the "write data to file" stuff out of batchavg and into its own proc. Then you'd pull them together with

Code: Select all

strdef fname
fname = "name_for_your_file.dat"
proc doall() {
  batchavg(...)
  savetofile(fname)
}

Rationale: modular code is easier to write and debug. Also makes it easier to reuse for other purposes.

The Vector class has lots of procs and funcs that are very useful for analyzing data. Check it out in the Programmer's Reference (the "Quick Index" is particularly convenient).

[erhervert]

Done

Code: Select all

objref savdata, tempmatrix

proc savetofile() {
  savdata = new File()
  savdata.wopen("results.dat")
  savdata.printf("t avg\n")
  tempmatrix = new Matrix()
  tempmatrix.resize(avg.size(),2)
  tempmatrix.setcol(0, time)
  tempmatrix.setcol(1, avg)
  tempmatrix.fprint(savdata, " %g")
  savdata.close()
}

proc doall() {
  batchavg(10)
  savetofile()
}

Thanks again Ted.

[ted]

Hey, I've got it easy. You're the one who has to actually write code and fix it when it doesn't work.

[erhervert]

I found it addictive, feels great when the things start working good.