/*
A primitive error function class for multiple run fitter must have a
func efun(y_vec, x_vec)
which is the simulation result of the run. Also, it should register
itself with a call to parmfitness_efun_append("classname").
Finally it should supply rfile and wfile procedures for saving its
state.
*/
/*********************************************************************************
TrajDens_Fitness
Calculate the error in trajectory density between model and target. This
is an objective function based on the phase portrait of the trace, proposed
in
LeMasson, G and Maex, R. Introduction to equation solving and parameter
fitting. In: Computational Neuroscience: Realistic Modeling for
Experimentalists. E. de Schutter, ed. CRC Press, 2001.
Briefly, we partition the phase plane (dV/dt vs. V) into a rectangular grid.
In each gridblock, we count the number of times the voltage trace occupies
that block. This gives a large matrix containing the trajectory density of
the voltage trace.
The trajectory density error is the sum, over all grid blocks, of the mean
squared difference between the model and target densities. The error is
normalized by the number of points contained in the trajectory, as well as
the size of the grid blocks.
Christina Weaver, Dec. 2005
christina.weaver@mssm.edu
*********************************************************************************/
install_vector_fitness()
parmfitness_efun_append("TrajDens_Fitness")
begintemplate TrajDens_Fitness
// PUBLIC declarations
//
public efun, set_data, xdat, ydat, have_data, use_x
public dV_tgt, dV_mdl, tag
public unmap, map, vbox, rfile, wfile
public clone, build, mintstop, before_run
public save_context, restore_context
public get_outname, set_r
public boundary, outname
public scale, Vbds, DVbds, DMtx_mdl, DMtx_tgt
public EFUN_DBG, VERBOSE, GBAR_SMRY
public Vbds, DVbds, delV, delDV
// VARIABLE declarations
//
// general variables
objref xdat, ydat // authoritative target data
objref dV_mdl, dV_tgt // first derivative of voltage data.
objref g, vbox, tobj, this
strdef tstr, mserrlabel, modelabel, scalelabel, tag, tmpstr
objref boundary
// variables related to Trajectory Density
objref Vbds, DVbds, DMtx_mdl, DMtx_tgt, DMtx2_mdl, DMtx2_tgt
objref NewGridInd, DMtx_diff
objref xmesh, ymesh, left_ptr, right_ptr, mid_ptr
objref GridTgt_x, GridTgt_y, GridMdl_x, GridMdl_y
objref GridTgt2_x, GridTgt2_y, GridMdl2_x, GridMdl2_y
// output variables
strdef outname, outMname, tab, tmpname
objref xtmp, ytmp, tmp_idx
objref dbgfile, dbgfile2
EFUN_DBG = 0
tab = " "
if( EFUN_DBG ) { print "Printing efun debug messages" }
mLidx = tMidx = mRidx = 0
tLidx = tMidx = tRidx = 0
nV = nDV = 0
external clipboard_set, clipboard_get
i=0
proc before_run() {}
/*******************************************************************************
efun() The main error function called by the Multiple Run Fitter.
This is an analog of calc_trajectory_density_error from the matlab code.
INPUT: $o1 y-values of model output (e.g soma.v(.5) )
$o2 t-values of model output (time, ms)
The model data is always changing when this function is called. Target data,
on the other hand, is unchanged unless new target data is loaded (or if
parameters defining the Density Grid are changed). Therefore, the Trajectory
Density Matrix of the target data does not need to be evaluated at each step,
but the Trajectory Density Matrix of the model data does. Then, the mean
squared difference between the elements of the two matrices are evaluated.
********************************************************************************/
func efun() { local tot_err, one_err, i
if( EFUN_DBG ) { print_Mfile_header($o1,$o2) }
// calculate trajectory density for model and (if needed) target
//
left_ptr = new Pointer(&mLidx)
mid_ptr = new Pointer(&mMidx)
right_ptr = new Pointer(&mRidx)
get_dVdt($o2,$o1,left_ptr,right_ptr,dV_mdl,mid_ptr)
calc_trajectory_density($o2,$o1,DMtx_mdl,GridMdl_x,GridMdl_y,left_ptr,right_ptr,dV_mdl,mid_ptr,DMtx2_mdl,GridMdl2_x,GridMdl2_y)
if( GridTgt_x.size() == 0 ) {
left_ptr = new Pointer(&tLidx)
mid_ptr = new Pointer(&tMidx)
right_ptr = new Pointer(&tRidx)
get_dVdt(ydat,xdat,left_ptr,right_ptr,dV_tgt,mid_ptr)
calc_trajectory_density(xdat,ydat,DMtx_tgt,GridTgt_x,GridTgt_y,left_ptr,right_ptr,dV_tgt,mid_ptr,DMtx2_tgt,GridTgt2_x,GridTgt2_y)
}
// merge all components of GridTgt into GridMdl
merge_GridLists(GridMdl_x,GridMdl_y,GridTgt_x,GridTgt_y)
merge_GridLists(GridMdl2_x,GridMdl2_y,GridTgt2_x,GridTgt2_y)
if( VERBOSE ) {
DMtx_diff = new Matrix(DMtx_mdl.nrow, DMtx_mdl.ncol)
}
tot_err = 0
for i = 0, GridMdl_x.size()-1 {
one_err = DMtx_tgt.x[GridMdl_x.x[i]][GridMdl_y.x[i]] - DMtx_mdl.x[GridMdl_x.x[i]][GridMdl_y.x[i]]
if( VERBOSE ) {
DMtx_diff.x[GridMdl_x.x[i]][GridMdl_y.x[i]] = one_err
}
tot_err += one_err*one_err
}
// to normalize by the number of grid blocks
// tot_err = tot_err / (nV*nDV)
printf("Trajectory Density Error = %g * %g = %g\n",tot_err,scale,tot_err*scale)
if( EFUN_DBG ) {
print_Mfile_tail($o1,$o2,tot_err)
}
if (use_gui) {
sprint(mserrlabel, "%g", tot_err*scale)
redraw($o2, $o1)
}
return tot_err * scale
}
/***************************************************************************
calc_trajectory_density
Calculate the trajectory density of the input trace. This is an
objective function based on the phase portrait of the trace, proposed
in
LeMasson, G and Maex, R. Introduction to equation solving and
parameter fitting. In: Computational Neuroscience: Realistic
Modeling for Experimentalists, E. de Schutter, ed. CRC Press,
2001.
input $o1 vector of time values (ms)
$o2 vector of V values (mV)
$o3 trajectory density matrix
$o4 X-coordinates of nonzero matrix entries
$o5 Y-coordinates of nonzero matrix entries
$o6 pointer to first time point within specified boundary
$o7 pointer to last time point within specified boundary
$o8 dV/dt will be stored in this vector
$o9 pointer to middle time point between specified boundaries
***************************************************************************/
proc calc_trajectory_density() { local nx, ny, val, ix, iy, Lidx, i
// july 18, 2011: started modifying this to allow for two regions for phase plane fit, within a single dV boundary. For the pre-CNS work, this is not worth the effort. Just add two
define_DMtx($o3,$o4,$o5)
Lidx = $o6.val
//
// loop through the trace, adding counts as necessary
for i = 0, $o8.size()-1 {
// first determine the indices; recall dVdt vector is 1 shorter
// than the V vector; start the V vector one index ahead.
ix = int( ($o2.x[Lidx+1+i]- Vbds.x[0]) / delV)
iy = int( ($o8.x[i] -DVbds.x[0]) / delDV)
// NOTE: if indices are outside of [0, nV-1] or [0, nDV-1], truncate.
if( ix < 0 ) ix = 0
if( iy < 0 ) iy = 0
if( ix >= nV ) ix = nV-1
if( iy >= nDV ) iy = nDV-1
/****
if( VERBOSE ) {
printf("[%g, %g] -> [%g, %g] -> (%d, %d)\n",$o2.x[Lidx+1+i],$o8.x[i],($o2.x[Lidx+1+i]- Vbds.x[0]) / delV,($o8.x[i] -DVbds.x[0]) / delDV,ix,iy)
dbgfile.printf("%% [%g, %g] -> [%g, %g] -> (%d, %d)\n",$o2.x[Lidx+1+i],$o8.x[i],($o2.x[Lidx+1+i]- Vbds.x[0]) / delV,($o8.x[i] -DVbds.x[0]) / delDV,ix,iy)
}
****/
// now increment the matrix
if( $o3.x[ix][iy] == 0 ) {
$o4.append(ix)
$o5.append(iy)
}
$o3.x[ix][iy] += 1
}
// normalize by the number of data points
$o3.muls(1/$o8.size())
}
/***************************************************************************
get_dVdt
Calculate the first derivative of the voltage data. Also identify
starting and ending points of the voltage vector within the specified
boundary.
input $o1 vector of time values (ms)
$o2 vector of V values (mV)
$o3 pointer to first time point within specified boundary
$o4 pointer to last time point within specified boundary
$o5 dV/dt will be stored in this vector
$o6 pointer to midpoint between specified boundaries
***************************************************************************/
proc get_dVdt() { local lpt, rpt, mpt
lpt = $o1.indwhere(">=",boundary.x[0])
rpt = $o1.indwhere(">=",boundary.x[2])
mpt = $o1.indwhere(">=",boundary.x[1])
if( rpt == -1 ) {
rpt = $o1.size()-1
} else { rpt -= 1 }
$o5 = new Vector()
$o5.deriv($o2.c(lpt,rpt),$o1.x[1]-$o1.x[0],2)
// $o5.deriv($o2.c(lpt,rpt),$o1.x[1]-$o1.x[0],1)
if( VERBOSE ) {
dbgfile.printf("dVmanual = [\n")
for i = lpt, rpt {
dbgfile.printf("%12.9f %12.9f %12.9f \n",$o2.x[i]-$o2.x[i-1],$o1.x[1]-$o1.x[0],($o2.x[i]-$o2.x[i-1])/($o1.x[1]-$o1.x[0]))
}
dbgfile.printf("];\n")
}
$o3.val = lpt
$o4.val = rpt
$o5.val = mpt
}
/*************************************************************************
merge_GridLists
input $o1 Master list of x-indices (GridMdl_x)
$o2 Master list of y-indices (GridMdl_y)
$o3 x-indices to merge (GridTgt_x)
$o4 y-indices to merge (GridTgt_y)
*************************************************************************/
proc merge_GridLists() { local ix, jx, mtch
NewGridInd = new Vector()
// first, look for repeat occurrences
for i = 0, $o3.size()-1 {
mtch = 0
for j = 0, $o1.size()-1 {
if( $o1.x[j] == $o3.x[i] && $o2.x[j] == $o4.x[i] ) {
mtch = 1
}
if( mtch ) break
}
if( mtch == 0 ) NewGridInd.append(i)
}
// now, append new elements the Master vectors
for i = 0, NewGridInd.size()-1 {
$o1.append($o3.x[NewGridInd.x[i]])
$o2.append($o4.x[NewGridInd.x[i]])
}
}
proc save_context() { local i
$o1.pack(tag, scale, ydat, \
ydat.label, xdat, boundary)
$o1.pack(delV, delDV)
$o1.pack(Vbds)
$o1.pack(DVbds)
}
proc restore_context() { local i
$o1.unpack(tag, &scale, ydat)
$o1.unpack(tstr, xdat)
ydat.label(tstr)
set_data(xdat, ydat)
$o1.unpack(boundary)
$o1.unpack(&delV, &delDV)
$o1.unpack(Vbds)
$o1.unpack(DVbds)
}
proc init() {local i
print "initializing TrajDens_Fitness, last modified 6 Dec 05"
EFUN_DBG = 0
VERBOSE = 0
GBAR_SMRY = 0
sprint(tag, "%s", this)
sscanf(tag, "%[^[]", tag)
use_x = 0
use_gui = 0
have_data = 0
xdat = new Vector(0)
ydat = new Vector(0)
boundary = new Vector(3)
outname = "TrajDens_data"
scale = 5e4
delV = 5
delDV = 25
//
// Define grid boundaries for density measure
Vbds = new Vector(2)
Vbds.x[0] = -80
Vbds.x[1] = 50
DVbds = new Vector(2)
DVbds.x[0] = -600
DVbds.x[1] = 600
// These two pairs of vectors store indices of non-zero DMtx* elements
GridMdl_x = new Vector()
GridMdl_y = new Vector()
GridTgt_x = new Vector()
GridTgt_y = new Vector()
GridMdl2_x = new Vector()
GridMdl2_y = new Vector()
GridTgt2_x = new Vector()
GridTgt2_y = new Vector()
define_DMtx(DMtx_mdl,GridMdl_x,GridMdl_y)
define_DMtx(DMtx_tgt,GridTgt_x,GridTgt_y)
define_DMtx(DMtx2_mdl,GridMdl2_x,GridMdl_y)
define_DMtx(DMtx2_tgt,GridTgt2_x,GridTgt_y)
sprint(scalelabel, "scale=%g del=(%g, %g)",\
scale,delV,delDV)
if (have_data) {
boundary.x[0] = xdat.x[0]
boundary.x[1] = xdat.x[int(0.5*(xdat.size()-1))]
boundary.x[2] = xdat.x[xdat.size()-1]
} else {
boundary.x[0] = boundary.x[1] = boundary.x[2] = 0
}
//print "\tdone initializing TrajDens_Fitness"
}
proc clone() {
$o1 = new TrajDens_Fitness()
$o1.have_data = have_data
if (have_data) {
$o1.set_data(xdat, ydat)
}
$o1.boundary = boundary.c
$o1.Vbds = Vbds.c
$o1.DVbds = DVbds.c
$o1.scale = scale
$o1.DMtx_mdl = DMtx_mdl.c
$o1.DMtx_tgt = DMtx_tgt.c
$o1.GridMdl_x = GridMdl_x.c
$o1.GridMdl_y = GridMdl_y.c
$o1.GridTgt_x = GridTgt_x.c
$o1.GridTgt_y = GridTgt_y.c
$o1.delV = delV
$o1.delDV = delDV
sprint(scalelabel, "scale=%g del=(%g, %g)",\
scale,delV,delDV)
}
proc redraw() { local i, xpsn
if (use_gui) {
g.erase()
if (have_data) {
g.label(.8, .95)
ydat.plot(g, xdat, 2, 1)
for i=0, boundary.size() - 1 {
g.beginline(3, 1)
g.line(boundary.x[i], g.size(3))
g.line(boundary.x[i], g.size(4))
}
}
if (numarg() == 2) {
$o2.line(g, $o1)
}
g.flush()
}
}
func mintstop() {
return boundary.x[2]
}
proc set_data() {local i
have_data = 0
i = $o1.indwhere(">=", 0)
if (i < 0 || $o1.size < 1) return
// copy $o1 into xdat without label;
// copy $o2 into ydat with its label string
xdat = $o1.c(i)
ydat = $o2.cl(i)
boundary.x[0] = xdat.x[0]
boundary.x[1] = xdat.x[int(0.5*(xdat.size-1))]
boundary.x[2] = xdat.x[xdat.size-1]
have_data = 1
if (use_gui) {
g.size(xdat.min(), xdat.max(), ydat.min(), ydat.max())
redraw()
}
}
proc clipboard_data() {
sprint(tstr, "%s.set_data(hoc_obj_[1], hoc_obj_[0])", this)
if(execute1(tstr) == 0) {
continue_dialog("No data in the Vector clipboard. Select a Graph line first")
}
}
proc build() {
if (use_gui) return
use_gui = 1
vbox = new VBox(3)
vbox.ref(this)
sprint(tstr, "execute(\"%s.unmap()\")", this)
vbox.dismiss_action(tstr)
vbox.save("save()")
vbox.intercept(1)
g = new Graph(0)
xpanel("", 1)
g.menu_tool("Adjust", "adjust_region")
xmenu("Select")
xbutton("Data from Clipboard", "clipboard_data()")
xbutton("Set weight", "set_weight()")
xbutton("Set grid spacing","set_grid_spacing()")
xbutton("Set grid boundaries","set_grid_bounds()")
xbutton("Region panel", "region_panel()")
xbutton("Set output info", "output_panel()")
xmenu()
xvarlabel(modelabel)
sprint(scalelabel, "scale= %g ", scale)
xvarlabel(scalelabel)
mserrlabel="MeanSqErr xxxxxxxxxxx"
xvarlabel(mserrlabel)
xpanel()
g.view(0, -80, 5, 40, 0,300,0,200)
if (have_data) {
g.size(xdat.min(), xdat.max(), ydat.min(), ydat.max())
redraw()
}
vbox.intercept(0)
}
proc adjust_region() {local x
//print $1, $2, $3
if ($1 == 2) { // press
adjust = pick_region($2)
set_r()
}
if (adjust == -1) {
return
}
if ($1 == 1) { // drag
boundary.x[adjust] = $2
if (adjust < boundary.size()-1) if ($2 > boundary.x[adjust+1]){
boundary.x[adjust] = boundary.x[adjust+1]
}
set_r()
}
if ($1 == 3) { // release
x = g.view_info(g.view_info, 11, $2)
if (boundary.size() > 2) {
if (x < 0) {
boundary.remove(0)
}else if (x > 1) {
boundary.remove(boundary.size()-1)
}else{
x = boundary.x[adjust]
tobj = boundary.c.indvwhere("==", x)
if (tobj.size() > 1) {
boundary.remove[adjust]
}
}
}
boundary.sort()
set_r()
adjust = -1
}
}
func pick_region() {local vi, d, i, j, x, m
vi = g.view_info
x = g.view_info(vi, 13, $1)
for i=0, boundary.size() - 1 {
d = x - g.view_info(vi, 13, boundary.x[i])
if (d < 12) { // before or on i
break
}
}
/***
if (i == boundary.size()) {
boundary.append($1)
return i
}
***/
if (d > -12) { // actual selection of line
return i
}
/***
boundary.insrt(i, $1)
if (i == 0) {
weight.insrt(1, weight.x[1])
}else{
weight.insrt(i, weight.x[i])
}
***/
return i
}
proc set_weight() {local wt
wt = scale
sprint(tmpstr, "%g", wt)
while (1) {
if (string_dialog("Enter function weight",tmpstr)){
if (sscanf(tmpstr, "%g", &scale) == 1) {
sprint(scalelabel, "scale=%g",scale)
return
}
}else{
break
}
}
scale = wt
}
proc set_grid_spacing() {local vsp, dvsp
vsp = delV
dvsp = delDV
sprint(tmpstr, "%g %g", vsp,dvsp)
while (1) {
if (string_dialog("Enter space separated V- and dV/dt-direction density grid spacing",tmpstr)){
if (sscanf(tmpstr, "%g %g", &delV, &delDV) == 2) {
sprint(scalelabel, "delV=%g delDV=%g",delV,delDV)
return
}
}else{
break
}
}
delV = vsp
delDV = dvsp
}
proc set_grid_bounds() {local vL, vR, dvL, dvR
vL = Vbds.x[0]
vR = Vbds.x[1]
dvL = DVbds.x[0]
dvR = DVbds.x[1]
dvsp = delDV
sprint(tmpstr, "%g %g %g %g", vL,vR,dvL,dvR)
while (1) {
if (string_dialog("Enter space separated V- and dV/dt-direction min/max density grid bounds",tmpstr)){
if (sscanf(tmpstr, "%g %g %g %g",&Vbds.x[0],&Vbds.x[1],&DVbds.x[0],&DVbds.x[1]) == 4) {
sprint(scalelabel, "Vbds=[%g %g]; DVbds=[%g %g]",Vbds.x[0],Vbds.x[1],DVbds.x[0],DVbds.x[1])
return
}
}else{
break
}
}
Vbds.x[0] = vL
Vbds.x[1] = vR
DVbds.x[0] = dvL
DVbds.x[1] = dvR
}
proc map() {
if (!use_gui) build()
if (numarg() > 1) {
vbox.map($s1, $2, $3, $4, $5)
}else if (numarg() == 1){
vbox.map($s1)
}else{
vbox.map()
}
redraw()
}
proc unmap() {
}
proc save() {local i
vbox.save("load_file(\"e_trajdens.hoc\", \"TrajDens_Fitness\")}\n{")
vbox.save("ocbox_=new TrajDens_Fitness(1)")
vbox.save("}\n{object_push(ocbox_)}\n{build()")
if (object_id(xdat)) {
sprint(tstr, "xdat = new Vector(%d)", xdat.size)
vbox.save(tstr)
sprint(tstr, "ydat = new Vector(%d)", ydat.size)
vbox.save(tstr)
sprint(tstr, "ydat.label(\"%s\")", ydat.label)
vbox.save(tstr)
sprint(tstr, "for i=0,%d { xdat.x[i]=fscan() ydat.x[i]=fscan()}}",\
xdat.size - 1)
vbox.save(tstr)
for i=0,xdat.size-1 {
sprint(tstr, "%g %g", xdat.x[i], ydat.x[i])
vbox.save(tstr)
}
vbox.save("{set_data(xdat, ydat)}")
}else{
vbox.save("}")
}
vbox.save("{object_pop()}\n{")
g.save_name("ocbox_.g", 1)
}
func pick_x() {local i
// NOTE 'indwhere' (or even manual search above)
// does NOT always find the correct index!
// often returns the first index greater than the sought value,
// even if there is an element that is exactly equal to the sought value.
i = xdat.indwhere(">=", $1)
if (i == -1) i = 0
return i
}
proc wfile() { local chooseAP
// Whenever adding more data, be sure to update the line count here!
// add 2 lines (to report vector size, plus one more), + vector size
// extra 6 lines are for boundary, ap window, ap subwindow
// extra lines for which APs will be used.
$o1.printf("TrajDens_Fitness xdat ydat boundary Vbds DVbds (lines=%d) ",\
4 + 2*xdat.size() + 2 + 2 + 2 )
$o1.printf(" %g %g %g\n",\
scale,delV,delDV)
$o1.printf("|%s|\n", ydat.label)
$o1.printf("%d\n", xdat.size())
xdat.printf($o1)
ydat.printf($o1)
$o1.printf("2\n%g\n%g%g\n\n", boundary.x[0],boundary.x[1],boundary.x[2])
$o1.printf("2\n%g\n%g\n", Vbds.x[0],Vbds.x[1])
$o1.printf("2\n%g\n%g\n", DVbds.x[0],DVbds.x[1])
}
proc rfile() {local i, n
scale = $o1.scanvar
delV = $o1.scanvar
delDV = $o1.scanvar
sprint(scalelabel, "scale=%g",scale)
$o1.gets(tstr)
if (sscanf(tstr, "|%[^|]", tstr) == 1) {
ydat.label(tstr)
}
n = $o1.scanvar
if (n > 0) {
xdat.resize(n) ydat.resize(n)
xdat.scanf($o1, n)
ydat.scanf($o1, n)
set_data(xdat, ydat)
}
n = $o1.scanvar
boundary.resize(n)
boundary.scanf($o1, n)
n = $o1.scanvar
Vbds.resize(n)
Vbds.scanf($o1, n)
n = $o1.scanvar
DVbds.resize(n)
DVbds.scanf($o1, n)
}
proc region_panel() {
xpanel("Region boundaries")
xpvalue("first interval startpoint", &boundary.x[0], 1, "set_r()")
xpvalue("mid-interval separation", &boundary.x[1], 1, "set_r()")
xpvalue("second interval endpoint", &boundary.x[2], 1, "set_r()")
xpanel()
}
proc output_panel() {
xpanel("Output info")
xcheckbox("Write M-file", &EFUN_DBG)
xcheckbox("Write verbose M-file", &VERBOSE)
xcheckbox("Write conductance summary (Av-Ron model only)", &GBAR_SMRY)
//xcheckbox("Write M-file", &EFUN_DBG, "EFUN_DBG = !EFUN_DBG")
//xcheckbox("Write verbose M-file", &VERBOSE, "{ VERBOSE = !VERBOSE EFUN_DBG = 1}")
xbutton("Set output filename","get_outname(outname)")
if( VERBOSE ) { EFUN_DBG = 1 }
if( !ismembrane("fn") ) { GBAR_SMRY = 0 }
xpanel()
}
proc set_r() {local i, t, tmin, tmax, n
if (have_data){
// make sure regions are within data boundaries
tmin = xdat.x[0]
tmax = xdat.x[xdat.size() - 1]
n = boundary.size()
for i=0, n-1 {
t = boundary.x[i]
if (t < tmin) {
boundary.x[i] = tmin
}
if (t > tmax) {
boundary.x[i] = tmax
}
}
}
redraw()
}
proc get_outname() {
string_dialog("Enter output filename",$s1)
}
/***************************************************************************
define_DMtx()
initialize the trajectory density matrix and the GridList vectors.
input $o1 trajectory density matrix
$o2 X-coordinates of nonzero matrix entries
$o3 Y-coordinates of nonzero matrix entries
***************************************************************************/
proc define_DMtx() {
//
// set up the grid
val = (Vbds.x[1]-Vbds.x[0]) / delV
if( int(val) == val ) {
nV = val
} else {
nV = int(val) + 1
}
val = (DVbds.x[1]-DVbds.x[0]) / delDV
if( int(val) == val ) {
nDV = val
} else {
nDV = int(val) + 1
}
if( VERBOSE ) {
printf("\tV bounds %g - %g will be divided into %d bins of size %g\n",\
Vbds.x[0],Vbds.x[1],nV,delV)
printf("\tDV bounds %g - %g will be divided into %d bins of size %g\n",\
DVbds.x[0],DVbds.x[1],nDV,delDV)
}
$o1 = new Matrix(nV,nDV)
$o2 = new Vector()
$o3 = new Vector()
}
proc print_Mfile_header() {
print "\nPrinting error function data to M-file ",outname
print "Scale factor = ", scale
sprint(outMname,"%s.m",outname)
dbgfile = new File()
dbgfile.wopen(outMname)
if( GBAR_SMRY ) {
dbgfile.printf("gbar_na=%g; gbar_k=%g; gbar_kca=%g;\n",gnabar_fn,gkbar_fn,gbar_kca)
dbgfile.printf("gbar_ka=%g; gbar_ca=%g; gbar_nap=%g;\n",gbar_ka,gbar_cahi,gbar_nap)
dbgfile.printf("Kp_cad=%g; Rca_cad=%g; ca0=%g;\n",Kp_cad,Rca_cad,cainf_cad)
}
nlbl = 1
dbgfile.printf("lbl(%d) = {'scale = %g, delV = %g, delDV = %g'};\n",nlbl,scale,delV, delDV)
nlbl += 1
dbgfile.printf("scale = %g; delV = %g; delDV = %g;\n",scale,delV, delDV)
dbgfile.printf("lbl(%d) = {'Vbds [%g, %g]'};\n",nlbl,Vbds.x[0],Vbds.x[1])
nlbl += 1
dbgfile.printf("lbl(%d) = {'DVbds [%g, %g]'};\n",nlbl,DVbds.x[0],DVbds.x[1])
nlbl += 1
dbgfile.printf("Vbds = [%g %g];\n",Vbds.x[0],Vbds.x[1])
dbgfile.printf("DVbds = [%g %g];\n",DVbds.x[0],DVbds.x[1])
dbgfile.printf("lbl(%d) = {' '};\n",nlbl)
nlbl += 1
if( GBAR_SMRY ) {
dbgfile.printf("lbl(%d) = {'g_{na} = %g g_{k-dr} = %g'};\n",nlbl,gnabar_fn,gkbar_fn)
nlbl += 1
dbgfile.printf("lbl(%d) = {'g_{k-ca} = %g g_{k-a} = %g'};\n",nlbl,gbar_kca,gbar_ka)
nlbl += 1
dbgfile.printf("lbl(%d) = {'g_{ca} = %g g_{nap} = %g'};\n",nlbl,gbar_cahi,gbar_nap)
nlbl += 1
dbgfile.printf("lbl(%d) = {'Kp = %g Rca = %g'};\n",nlbl,Kp_cad,Rca_cad)
nlbl += 1
dbgfile.printf("lbl(%d) = {'[Ca]_0 = %g'};\n",nlbl,cainf_cad)
nlbl += 1
}
dbgfile.printf("boundary = [%g %g %g];\n",boundary.x[0],boundary.x[1],boundary.x[2])
if( VERBOSE ) {
dbgfile.printf("xmodel = [")
$o2.printf(dbgfile,"%12.9f\n ")
dbgfile.printf("];\n")
dbgfile.printf("ymodel = [")
$o1.printf(dbgfile,"%12.9f\n ")
dbgfile.printf("];\n")
dbgfile.printf("xexpt = [")
xdat.printf(dbgfile,"%12.9f\n ")
dbgfile.printf("];\n")
dbgfile.printf("yexpt = [")
ydat.printf(dbgfile,"%12.9f\n ")
dbgfile.printf("];\n")
}
nstr = 1
}
proc print_Mfile_tail() { local e
e = $3
dbgfile.printf("e_final = %g;\n",e)
dbgfile.printf("str(%d) = {'Total Error %g'};\n",nstr,e)
nstr += 1
dbgfile.printf("figure(1);\n")
dbgfile.printf("h = axes('Position',[0 0 1 1],'Visible','off');\n")
dbgfile.printf("ttl = '%s: Total Error %g';\n",outname,e)
if( VERBOSE ) {
dbgfile.printf("axes('Position',[.1 .5 .8 .4]);\n")
dbgfile.printf("x1=[boundary(1) boundary(1)];\n")
dbgfile.printf("x2=[boundary(2) boundary(2)];\n")
dbgfile.printf("x3=[boundary(3) boundary(3)];\n")
dbgfile.printf("y1=[min(yexpt) max(yexpt)];\n")
dbgfile.printf("pl = plot(xexpt,yexpt,'--',xmodel,ymodel,'-',x1,y1,'r',x2,y1,'r',x3,y1,'r');\n")
dbgfile.printf("set(pl(1),'LineWidth',2); set(pl(2),'LineWidth',2); \n")
dbgfile.printf("legend('target','model');\n")
dbgfile.printf("title(ttl);\n")
}
dbgfile.printf("set(gcf,'CurrentAxes',h);\n")
dbgfile.printf("text(.05,.25,lbl,'FontSize',12);\n")
dbgfile.printf("text(.5,.25,str,'FontSize',12);\n")
if( VERBOSE ) {
dbgfile.printf("orient(gcf,'tall');saveas(gcf,'%s_V.fig');\n\n",outname)
/** set up the trajectory density plots, plus the phase portrait. **/
xmesh = new Vector()
ymesh = new Vector()
xmesh.indgen(Vbds.x[0],Vbds.x[1]-delV,delV)
ymesh.indgen(DVbds.x[0],DVbds.x[1]-delDV,delDV)
// now print out the xmesh data
dbgfile.printf("xmesh = [")
xmesh.printf(dbgfile)
dbgfile.printf("];\n")
dbgfile.printf("ymesh = [")
ymesh.printf(dbgfile)
dbgfile.printf("];\n")
// plot phase portraits of target and model, together.
dbgfile.printf("dV_mdl = [")
dV_mdl.printf(dbgfile,"%12.9f\n ")
dbgfile.printf("];\n")
dbgfile.printf("dV_tgt = [")
dV_tgt.printf(dbgfile, "%12.9f\n ")
dbgfile.printf("];\n")
dbgfile.printf("mLidx = %d; mRidx = %d; tLidx = %d; tRidx = %d;\n",\
mLidx+1,mRidx+1,tLidx+1,tRidx+1)
dbgfile.printf("\nfigure(2); \n")
dbgfile.printf("plot(yexpt(tLidx:tRidx),dV_tgt,'b.',ymodel(mLidx:mRidx),dV_mdl,'r.');\n")
//dbgfile.printf("plot(yexpt(tLidx+1:tRidx),dV_tgt,'b.',ymodel(mLidx+1:mRidx),dV_mdl,'r.');\n")
// show gridlines for density measure
dbgfile.printf("hold on;\n")
dbgfile.printf("for i = 1 : length(xmesh)\n")
dbgfile.printf("\tplot([xmesh(i) xmesh(i)], DVbds, 'k-');\n")
dbgfile.printf("end;\n\n")
dbgfile.printf("for i = 1 : length(ymesh)\n")
dbgfile.printf("\tplot(Vbds, [ymesh(i) ymesh(i)], 'k-');\n")
dbgfile.printf("end;\n\n")
dbgfile.printf("xlabel('membrane potential (mV)'); ylabel('dV/dt (mV/ms)');\n")
dbgfile.printf("legend('Target','Model');\n")
dbgfile.printf("phase_ttl = '%s : Voltage trajectory phase plane'; title(phase_ttl);\n",outname)
dbgfile.printf("figname = 'phaseplane_%sVsTgt.fig'; saveas(gcf,figname);\n",outname)
dbgfile.printf("hold off;\n")
// also print out GridList and DMtx_diff, the density difference between the two traces.
// (recall, indices in NEURON start with 0, matlab start with 1)
//
dbgfile.printf("DMtx_diff = [\n%% ")
DMtx_diff.fprint(dbgfile)
dbgfile.printf("];\n\n")
dbgfile.printf("GridList = [\n")
for i = 0, GridMdl_x.size()-1 {
dbgfile.printf("%d %d\n",GridMdl_x.x[i]+1,GridMdl_y.x[i]+1)
}
dbgfile.printf("];\n\n")
dbgfile.printf("figure(3);\n")
dbgfile.printf("print_trajectory_density(DMtx_diff,GridList,'model',xmesh,ymesh,0);\n")
dbgfile.printf("%% Trajectory Density Error = %g\n",e)
dbgfile.printf("tot_err = %g;\n",e)
dbgfile.close()
}
}
endtemplate TrajDens_Fitness