/*
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
*******************************************************************/
/* error is independent of position of Action potential peak */
begintemplate TrajDens_fit2
// PUBLIC declarations
//
public read_ExpDat
public efun, g, set_data, set_modelx, xdat, ydat, idx_, peak, have_data
public rfile, wfile, use_x, peakscale, xscale, yscale, ntag, tag
public clone, mintstop
public find_APtimes, xyAPtimes, modelAPtimes, findFR, get_outname
public set_r
public efun_1AP, efun_withFRCV, efun_old, efun_yfit, efun_isiW
public datFR, modFR // set of all instantaneous FR's
public datFRmean, datCVmean, modFRmean, modCVmean // mean values
public datFRsd, datCVsd, modFRsd, modCVsd // SD values
public boundary, apwin, outname
public t_hat, pow1, pow2, frscale, cvscale, ey, ey2
public frm_fac, frs_fac, cvm_fac, cvs_fac, frcv_int
public xdat_rel_list, ydat_rel_list, idx_rel_list
public time_scale, miss_scale, shape_scale, brst_thold, brst_pnlty
public set_apwin_noGUI, efun_mismatch_avgAP, efun_bower
public EFUN_DBG, VERBOSE, GBAR_SMRY, TIME_DBG
public PASTE_DBG
public track_APtimes, allList
public bound_2exp
// public variables for time-modulated FR error
public tot_frscale, nf_frscale, slope_scale, intcp_scale, intcp2_scale, slope2_scale
public exp_intcpt, exp_slope, expsa, expsb, expchi
public mod_intcpt, mod_slope, modsa, modsb, modchi
public datFit, modFit, efr
public e_dly, delay_pnlty, ebrst, e_slope
//public variables for AHP measurement
public neg_p, trgh_sz, trgh_ps, thr_ac
public apA, apS, apD // AP count in Axon, Soma, distal Dendrite
public spktA, spktS, spktD // spike times in Axon, Soma, distal Dendrite
public f_exp, exp_vec, local_vec, local_sf
// general variables
//
objref xdat, ydat // authoritative data
objref g, vbox, tobj, this
strdef tstr, mserrlabel, modelabel, scalelabel, ntag, tmpstr, tag
objref boundary
// variables related to shape calculations
objref xdat_rel_list, ydat_rel_list, idx_rel_list // subset of data points to fit
objref xyAPtimes, modelAPtimes // AP times of (xdat, ydat) and of the model result
objref tgtAPtrough, modelAPtrough // troughs after APs of (xdat, ydat) and of the model result
objref tgtAPtrough_V, modelAPtrough_V // troughs after APs of (xdat, ydat) and of the model result
objref apwin
objref left_ptr, right_ptr, n_ptr, mid_ptr
objref Png, Pts, Pps, Pac, nw_ptr
objref dvdt, allList
// variables for FR & CV calculations
objref datFR, modFR
objref frtmp, frptr, cvptr, frptrsd, cvptrsd // mean, stdev of FR & CV
objref frindx, frtmpsd, cvtmp
//variables for time-modulated FR calculations
objref ptra, ptrb, ptrsiga, ptrsigb, ptrchi
objref datFit, modFit
// variables to measure AHP
objref dV, iAHP, iAHPn, iAHPa, AHPchg
objref APth, APst, APnd, npts
objref trough_sz,trough_posn
objref AHPout
strdef AHPname
// variables to check AP firing order
objref apA, apS, apD // AP count in Axon, Soma, distal Dendrite
objref spktA, spktS, spktD // spike times in Axon, Soma, distal Dendrite
strdef dist_name
// output variables
strdef outname, outMname, tab, tmpname
objref xtmp, ytmp, tmp_idx
objref dbgfile, dbgfile2
objref f_exp, exp_vec[2], local_sf
datFRmean = datCVmean = modFRmean = modCVmean = 0
datFRsd = datCVsd = modFRsd = modCVsd = 0
exp_APleft = exp_APright = mod_APleft = mod_APright = 0
exp_APmid = mod_APmid = 0
n_expAP = n_modAP = 0
n_expWin = n_modWin = 0
EFUN_DBG = 0
PASTE_DBG = 0
SAVE_DBG = 0
tab = " "
if( SAVE_DBG ) { print "Printing save debug messages" }
if( PASTE_DBG ) {
dbgfile = new File()
dbgfile.wopen("paste_curve.m")
dbgfile.printf("%% trying to understand paste mechanism\n")
dbgfile.close()
}
if( EFUN_DBG ) { print "Printing efun debug messages" }
i=0
local_sf = new StringFunctions()
/*****************************************************
for now: assume data are stored as NRN .dat files.
Code modified from clipboard_retrieve(), found in ~/neuron/share/nrn/lib/hoc/stdlib.hoc
input $s1 .dat name where exp data is found
//$o2 read time vector
//$o3 read voltage vector
*****************************************************/
proc read_ExpDat() { local b, ii, jj, size
// read experimental data here, into (xdat,ydat). Same as 'retrieve from file' in NEURON menu.
f_exp = new File()
bsize = 200
b = f_exp.ropen($s1)
if (b) {
for ii = 0, 1 {
exp_vec[ii] = new Vector()
}
// first line may be label:str or single number
// otherwise start the data immediately
f_exp.gets(tstr)
if (sscanf(tstr, "label:%s\n", tstr) == 1) {
f_exp.gets(tstr)
}
ii = sscanf(tstr, "%g", &size)
jj = 1
if (ii == 1) {
for ii=0, 1 {
exp_vec[ii].buffer_size(size)
exp_vec[ii].resize(0)
}
jj = 0
}else{
size = 1e20
}
while(!f_exp.eof()) {
if (jj >= bsize) {
bsize *= 2
for ii=0, 1 {exp_vec[ii].buffer_size(bsize)}
}
for ii=0,1 {
exp_vec[ii].append(f_exp.scanvar())
}
jj += 1
if (jj > size) { break }
}
f_exp.close()
} else {
printf("Reading experimental data from %s failed",$s1)
}
//$o2 = exp_vec[0].c
//$o3 = exp_vec[1].c
xdat = exp_vec[0].c
ydat = exp_vec[1].c
have_data = 1
}
proc before_run() {}
/*******************************************************
func efun() Aug 04
Let Nt = number of APs in experiment
Nm = number of APs in model
YFit_i = result from efun_yfit(model spike i)
Then, the error function calculated here is the E(M) described in our paper (in
preparation). Namely, there are three components to this function:
1. for each spike which occurs in both exp & model (i.e. number is
min(Nt,Nm)), calculate the absolute difference in spike timing, and divide
by scale factor Ai. Tally the sum of these errors.
2. for each missing/extra spike, add P/Nt, where P=tstop
3. for each model spike, add Yfit_i, making sure that it is scaled
appropriately as a fraction of the window size (e.g. multiply scale factor
from 0 to 1)
Notes on YFit_i: we already truncate the window appropriately if it occurs at
the beginning or end of the simulation. But also, we do not want to include
the previous or next AP in the window. In such cases, set the cutoff to be the
recorded "spiketime", namely the time at which the voltage crosses -20 mV.
This will avoid any error due to the interfering AP.
If there are no model "APs", YFit_i is determined as the
YFit from the first model peak, even though it is not
actually an AP (then Nm in denominator is 1, but Nm in
numerator is 0).
See efun_yfit below for a description of its evaluation.
This describes the shape error. The error function also has the follwoing terms:
* Bursting penalty for ISIs smaller than a specified threshold
* Mean and standard deviation of instantaneous firing rate
* Mean and standard deviation of coefficient of variation
* Firing delay penalty
*******************************************************/
func efun() { local e, etm, i, A, mnind, mxind, exind
//printf("EFUN (%g,%g,%g,%g,%g,%g,%g,%g)\n",gnabar_fn,gkbar_fn,gbar_kca,Kp_cad,Rca_cad,gbar_ka,gbar_cahi,gbar_nap)
/**
EFUN_DBG=1
GBAR_SMRY=1
VERBOSE=1
sprint(outname,"data_Na%.18f",gnabar_fn)
**/
// INPUT:
// $o1 y-values of model output (e.g soma.v(.5) )
// $o2 t-values of model output (time, ms)
if( EFUN_DBG ) {
print "\nPrinting error function data to M-file ",outname
print "Scale factor = ", shape_scale
sprint(outMname,"%s.m",outname)
dbgfile = new File()
dbgfile.wopen(outMname)
dbgfile.printf("%% finding FR for experiment, then model\n")
}
// NOTE AP times gives an index into the respective time vector,
// representing the time when the AP passes threshold, NOT the
// peak. Use nextpeak.
if( xyAPtimes.size == 0 ) {
left_ptr = new Pointer(&exp_APleft)
mid_ptr = new Pointer(&exp_APmid)
right_ptr = new Pointer(&exp_APright)
n_ptr = new Pointer(&n_expAP)
nw_ptr = new Pointer(&n_expWin)
if( EFUN_DBG ) printf("----------Finding experimental AP times:\n")
find_APtimes(xdat,ydat,xyAPtimes,"Exp AP",1,left_ptr,right_ptr,n_ptr,nw_ptr,mid_ptr)
} else {
if( xdat.x[xyAPtimes.x[exp_APleft]] < boundary.x[0] || xdat.x[xyAPtimes.x[exp_APright]] > boundary.x[1]) {
printf("Time boundary changed; repeat target AP search.\n")
left_ptr = new Pointer(&exp_APleft)
mid_ptr = new Pointer(&exp_APmid)
right_ptr = new Pointer(&exp_APright)
n_ptr = new Pointer(&n_expAP)
nw_ptr = new Pointer(&n_expWin)
if( EFUN_DBG ) printf("----------New time bds: Finding experimental AP times:\n")
find_APtimes(xdat,ydat,xyAPtimes,"Exp AP",1,left_ptr,right_ptr,n_ptr,nw_ptr,mid_ptr)
}
}
tgtAPtrough = new Vector()
tgtAPtrough_V = new Vector()
if( VERBOSE ) printf("calling APtrough for TARGET\n")
if( xyAPtimes.size() > 0 ) {
find_APtrough(xdat,ydat,xyAPtimes.c(exp_APleft,exp_APright),boundary,tgtAPtrough,tgtAPtrough_V)
}
//printf("Boundary is %g - %g and expAPtimes go from %g to %g\n",boundary.x[0],boundary.x[1],xdat.x[xyAPtimes.x[exp_APleft]],xdat.x[xyAPtimes.x[exp_APright]])
n_AP = 0
for i = 0, xyAPtimes.size-1 {
if( xdat.x[xyAPtimes.x[i]] < boundary.x[0] ) continue
if( xdat.x[xyAPtimes.x[i]] > boundary.x[1] ) { break }
if( xdat.x[xyAPtimes.x[i]] > $o2.x[$o2.size-1] ) { break }
n_AP += 1
}
//print "Found n_AP = ",n_AP
frptr = new Pointer(&datFRmean)
cvptr = new Pointer(&datCVmean)
frptrsd = new Pointer(&datFRsd)
cvptrsd = new Pointer(&datCVsd)
findFR(xdat.ind(xyAPtimes),datFR,frptr,frptrsd,cvptr,cvptrsd)
left_ptr = new Pointer(&mod_APleft)
mid_ptr = new Pointer(&mod_APmid)
right_ptr = new Pointer(&mod_APright)
n_ptr = new Pointer(&n_modAP)
nw_ptr = new Pointer(&n_modWin)
if( EFUN_DBG ) printf("----------Finding model AP times:\n")
find_APtimes($o2,$o1,modelAPtimes,"Model AP",1,left_ptr,right_ptr,n_ptr,nw_ptr,mid_ptr)
if( EFUN_DBG ) printf("Found n_modWin %d\n",n_modWin)
modelAPtrough = new Vector()
modelAPtrough_V = new Vector()
if( VERBOSE ) printf("calling APtrough for MODEL from %d to %d\n",mod_APright,mod_APleft)
if( mod_APright>=0 && mod_APleft>=0 && mod_APright-mod_APleft >= 0 ) {
find_APtrough($o2,$o1,modelAPtimes.c(mod_APleft,mod_APright),boundary,modelAPtrough,modelAPtrough_V)
} else {
if( VERBOSE ) printf("No model APs, so no troughs identified.\n")
}
if( EFUN_DBG ) {
dbgfile.printf("tgt_trough = [\n")
for i = 0,tgtAPtrough.size()-1 {
dbgfile.printf("\t%g\t%g\n",tgtAPtrough.x[i],tgtAPtrough_V.x[i])
}
dbgfile.printf("];\n")
dbgfile.printf("mdl_trough = [\n")
for i = 0,modelAPtrough.size()-1 {
dbgfile.printf("\t%g\t%g\n",modelAPtrough.x[i],modelAPtrough_V.x[i])
}
dbgfile.printf("];\n")
}
if( EFUN_DBG ) {
printf("Looked again for model APs from %g - (cut %g ) - %g; ",boundary.x[0],bound_2exp,boundary.x[1])
printf("Found %g to %g to %g; ", $o2.x[modelAPtimes.x[mod_APleft]], $o2.x[modelAPtimes.x[mod_APmid]],$o2.x[modelAPtimes.x[mod_APright]])
printf("%g spikes, %g windows\n",n_modAP,n_modWin)
}
frptr = new Pointer(&modFRmean)
cvptr = new Pointer(&modCVmean)
frptrsd = new Pointer(&modFRsd)
cvptrsd = new Pointer(&modCVsd)
findFR($o2.ind(modelAPtimes),modFR,frptr,frptrsd,cvptr,cvptrsd)
e = 0
nstr = 1
nstr2 = 1
nlbl = 1
if( EFUN_DBG ) { print_Mfile_header($o1,$o2,e) }
//printf("[1] e %g\t",e)
// Evaluate each component of the error function.
//
//
n_both = n_AP
if( n_modAP > n_AP ) {
n_both = n_AP
n_big = n_modAP
} else {
n_both = n_modAP
n_big = n_AP
}
if( n_modWin > n_expWin ) {
n_bigWin = n_modWin
n_bothWin = n_expWin
} else {
n_bothWin = n_modWin
n_bigWin = n_expWin
}
if( EFUN_DBG ) { dbgfile.printf("n_both = %d; n_big=%d;\n",n_both,n_big) }
//printf("[3] e %g\t",e)
// Penalize ISI's which are less than a user-specified threshold brst_thold,
// by a user specified amount brst_pnlty
//
// This will penalize doublets and bursting.
ebrst = 0
if( brst_pnlty > 0 ) {
ebrst = efun_penburst($o2)
//printf("ebrst = %g\t",ebrst)
e += ebrst
} else {
if( EFUN_DBG ) {
//dbgfile.printf("str(%d) = {'ISI burst penalty (0 times) 0'};\n", nstr)
//nstr += 1
dbgfile.printf("lbl(%d) = {'ISI burst penalty (0 times) 0'};\n", nlbl)
nlbl += 1
printf("ISI burst penalty (0 times) 0\n")
}
}
//printf("[4] e %g\t",e)
// Penalize the difference between the target and model delay to first spike.
// Scale factor on this term is delay_pnlty; for description of this term see
// the func description below
if( delay_pnlty > 0 ) {
e_dly = efun_pendelay($o1,$o2)
//printf("edly = %g\t",e_dly)
e += e_dly
} else {
if( EFUN_DBG ) {
//dbgfile.printf("str(%d) = {'Spike delay penalty (0 times) 0'};\n", nstr)
//nstr += 1
dbgfile.printf("lbl(%d) = {'Spike delay penalty (0 times) 0'};\n", nlbl)
nlbl += 1
printf("Spike delay penalty (0 times) 0\n")
}
}
//printf("[5] e %g\t",e)
//
// 1. for each model spike, add Yfit_i, making sure that it is scaled
// appropriately as a fraction of the window size (e.g. multiply scale factor
// from 0 to 1)
//
// If there are no model APs, just calculate the error from
// the peak model voltage. Otherwise, calculate the average yfitness
// error.
ey = 0
mnind = $o2.indwhere(">=",boundary.x[0])
if( mnind < 0 ) {
//printf("First boundary never reached; therefore NO SHAPE ERROR.\n")
return
}
mxind = $o2.indwhere(">=",boundary.x[1])
if( mxind < 0 ) {
mxind = $o2.size-1
//printf("Boundary not reached; set end to %g\n",$o2.x[mxind])
}
//printf("mnind %g, mxind %g\n",mnind,mxind)
// if match_all == 1, Match each model AP sequentially to the
// corresponding experimental AP.
// else, match each model AP to the each of the experimental APs
// selected.
//
//printf("[6] e %g\t",e)
calc_shape_error($o1,$o2)
//printf("ey = %g\t",ey)
//printf("[7] e %g\t",e)
if( EFUN_DBG ) {
//dbgfile.printf("str(%d) = {'Shape Error (mV-ms) %g'};\n",nstr,ey )
//nstr += 1
dbgfile.printf("lbl(%d) = {'Shape Error (mV-ms) %g'};\n",nlbl,ey )
nlbl += 1
if( shape_scale > 0 ) {
print "Shape Error (mV-ms)\t",ey/shape_scale," --> ",ey
} else { print "Shape Error (mV-ms)\t0 --> 0" }
}
e += ey
//printf("[8] e %g\t",e)
//
// For both the model & experimental data, fit a straight line to the
// instantaneous FR. Use the intercept & slope to calculate an
// appropriate error. Note: there must be at least three target APs
// (i.e. two FR points) to fit the line.
//
//if( n_AP > 3 ) { // try fitting even when there is no firing.
efr = 0
//if( mod_APmid > mod_APleft ) {
// early phase of adaptation
if( EFUN_DBG ) {
printf("------About to fit early phase SFA\n")
}
efr = efun_expfit($o1,$o2,exp_APleft,exp_APmid,mod_APleft,mod_APmid,intcp_scale,slope_scale,boundary.x[0])
//printf("efr1 = %g\n",efr)
if( EFUN_DBG ) {
printf("------DONE early phase SFA\n")
}
if( VERBOSE ) {
printf("Total Early FR Fit Error: %g = ",efr)
}
//printf("[9] e %g\t",e)
//} else { if( EFUN_DBG ) print "Skipping early phase SFA" }
e += efr
//if( mod_APmid < mod_APright ) {
//printf("[10] e %g\t",e)
if( EFUN_DBG ) {
printf("------About to fit late phase SFA\n")
}
// late phase of adaptation
efr = efun_expfit($o1,$o2,exp_APmid,exp_APright,mod_APmid,mod_APright,intcp2_scale,slope2_scale,bound_2exp)
//printf("efr1 = %g\n",efr)
if( EFUN_DBG ) {
printf("------DONE late phase SFA\n")
}
if( VERBOSE ) {
printf("Total Late FR Fit Error: %g = ",efr)
}
//printf("[9prime] e %g\t",e)
//} else { if( EFUN_DBG ) print "Skipping late phase SFA" }
e += efr // if there is no second phase of adaptation, add the error twice.
//printf("[10prime] e %g\t",e)
//}
//printf("[11] e %g\t",e)
if( EFUN_DBG ) {
print_Mfile_tail($o1,$o2,e)
}
//printf("[13] e %g\t",e)
if( EFUN_DBG ) {
dbgfile.printf("%% now e=%g\n",e)
dbgfile.close()
}
//if( VERBOSE) {
//printf("EFUN (%g,%g,%g,%g,%g,%g,%g,%g) [%g %g %g %g %g] %g %g %g %g: Returning e = %g\n",gnabar_fn,gkbar_fn,gbar_kca,Kp_cad,Rca_cad,gbar_ka,gbar_cahi,gbar_nap,shape_scale,tot_frscale,nf_frscale,delay_pnlty,slope_scale,e_dly,ey,efr,e_slope,e)
//}
return e
}
proc calc_shape_error() { local mnind, mxind, i, exind, loc_ind, mtchall, iloop
// $o1, $o2 same as for efun()
mtchall = 1
// mnind : first index of model AP window
// mxind : last index of model AP window
mnind = $o2.indwhere(">=",boundary.x[0])
mxind = $o2.indwhere(">=",boundary.x[1])
if( EFUN_DBG ) printf("In calc_shape_error, time indices are [%g, %g]\n",$o2.x[mnind],$o2.x[mxind])
if( n_modWin == 0 ) {
// there's no APs; just match the first model peak to the first (or
// each selected) AP peak
if( EFUN_DBG ) printf("peak of subset = %g, values %g - %g - %g\n",$o1.c(mnind,mxind).firstmax,$o1.x[mnind+$o1.c(mnind,mxind).firstmax-1],$o1.x[mnind+$o1.c(mnind,mxind).firstmax],$o1.x[mnind+$o1.c(mnind,mxind).firstmax+1])
peak = $o2.x[mnind+$o1.c(mnind,mxind).firstmax]
xpeak = xdat.x[ydat.firstpeak]
//printf("Time of first model peak is at %g, value %g\n",peak,$o1.x[mnind+$o1.c(mnind,mxind).firstmax])
if( (shape_scale > 0) && (apwin.x[1]-apwin.x[0] > 0) ) {
if( mtchall ) {
if( EFUN_DBG ) printf("Model shape error from peak:\n")
if( EFUN_DBG ) printf("No model win, but Match ALL: before efun_yfit, time indices are [%g, %g]\n",$o2.x[mnind],$o2.x[mxind])
ey += efun_yfit($o1,$o2,0,exp_APleft,0,$o2.size-1)
if( EFUN_DBG ) printf("(4) Return from this AP.\n")
} else {
for i = exp_APleft, exp_APright {
if( EFUN_DBG ) printf("Shape error from model peak against exp AP %d:\n",i)
if( EFUN_DBG ) printf("No model win, match expAP %d: before efun_yfit, time indices are [%g, %g]\n",i,$o2.x[mnind],$o2.x[mxind])
ey += efun_yfit($o1,$o2,0,i,0,$o2.size-1)
if( EFUN_DBG ) printf("(3) Return from this AP.\n")
if( EFUN_DBG && VERBOSE ) {
dbgfile.printf("%% No model APs\n")
dbgfile.printf("ey(%d) = %g;\t",i,ey)
}
}
}
}
if( EFUN_DBG && VERBOSE ) {
dbgfile.printf("\n%% No model APs\n")
dbgfile.printf("peak(1) = %g; xpeak(1) = %g;\n",peak,xpeak)
dbgfile.printf("ey(1) = %g;\n",ey)
}
} else {
// for each model AP, either match to each corresponding AP
// (if mtchall), or else match each model AP to each chosen
// experimental AP.
//printf("modelAPtimes, from %d to %d\n",mod_APleft,mod_APleft+n_modWin-1)
//for iloop = mod_APleft,mod_APright {
//printf("\t[%d] = %g\n",modelAPtimes.x[iloop],$o2.x[modelAPtimes.x[iloop]])
//}
for iloop = 0, n_modWin-1 {
//for i = mod_APleft, mod_APright {
i = mod_APleft + iloop
//printf("i= %d = %d + %d\n",i,mod_APleft,iloop)
peak = $o2.x[$o1.nextpeak(modelAPtimes.x[i])]
//printf("Found peak %g which is the time of the %d th model AP, index %d\n",peak,i,modelAPtimes.x[i])
// If we don't want to limit the AP Win boundaries
// for neighboring APs, comment out this section
if( i == 0 ) {
mnind = $o2.indwhere(">=",boundary.x[0])
if( mnind <= 0 ) mnind = 0
//printf("Set mnind = %d\n",mnind)
} else {
// this is where the last peak was, plus twice the time
// diff b/w the last peak and its threshold crossing.
// if this occurs later than the prescribed AP window,
// use this cutoff instead. (distance wasn't doubled
// before; this was added 12/14/04).
//printf("Look for x >= %g = 3*%g - 2*%g\n",3*$o2.x[$o1.nextpeak(modelAPtimes.x[i-1])]-2*$o2.x[modelAPtimes.x[i-1]],$o2.x[$o1.nextpeak(modelAPtimes.x[i-1])],$o2.x[modelAPtimes.x[i-1]])
mnind = $o2.indwhere(">=",3*$o2.x[$o1.nextpeak(modelAPtimes.x[i-1])]-2*$o2.x[modelAPtimes.x[i-1]])
//printf("**Set mnind = %d\n",mnind)
}
if( i == mod_APright ) {
mxind = $o2.indwhere(">=",boundary.x[1])
if( mxind <= 0 ) mxind = $o2.size-1
//printf("Set mxind = %d\n",mxind)
} else {
// upper bound of window is the threshold crossing of
// the next spike, minus the distance from that
// next threshold crossing and its AP peak.
// (this was changed 12/14/04; was = to next AP thold)
//printf("Look for x >= %g = 2*%g -%g\n",2*$o2.x[modelAPtimes.x[i+1]]-$o2.x[$o1.nextpeak(modelAPtimes.x[i+1])],$o2.x[modelAPtimes.x[i+1]],$o2.x[$o1.nextpeak(modelAPtimes.x[i+1])])
mxind = $o2.indwhere(">=",2*$o2.x[modelAPtimes.x[i+1]]-$o2.x[$o1.nextpeak(modelAPtimes.x[i+1])])
if( mxind < 0 ) {
mxind = $o2.size-1
//printf("\tNot found; use mxind = %d instead\n",mxind)
}
//printf("**Set mxind = %d\n",mxind)
}
if( mtchall ) {
// match model spike i with experimental spike i;
// if there are more model spikes, compare the extras
// against the last experimental spike
if( i-mod_APleft > n_expWin-2 ) {
exind = exp_APleft+n_expWin-1
} else {
exind = exp_APleft + (i - mod_APleft)
}
xpeak = xdat.x[ydat.nextpeak(xyAPtimes.x[exind])]
loc_ind = exind - exp_APleft
if( EFUN_DBG && VERBOSE ) {
dbgfile.printf("%% There are %d experimental APs\n",n_expAP)
dbgfile.printf("%% and %d experimental AP windows\n",n_expWin)
dbgfile.printf("%% Comparing target %d (peak %g) with model %d (peak %g)\n",exind,xpeak,i,peak)
dbgfile.printf("%% before efun_yfit\n")
dbgfile.printf("mnind(%d) = %d; mxind(%d) = %d;\n",i+1,mnind,i+1,mxind)
dbgfile.printf("exind(%d) = %d;\n",i+1,exind)
//print "Match model AP ",i," to experiment AP ",exind
if( i>0 ) {
dbgfile.printf("t_th(%d)=%g; t_pk(%d)=%g; t_apk(%d)=%g;\n",i+1,$o2.x[modelAPtimes.x[i-1]],i+1,$o2.x[$o1.nextpeak(modelAPtimes.x[i-1])],i+1,$o2.x[mnind])
}
dbgfile.printf("%% y direction: Y %d X %d\n",ydat_rel_list.object(loc_ind).size,\
xdat_rel_list.object(loc_ind).size)
dbgfile.printf("%% now y:\n ydat_rel%d = [",exind+1)
ydat_rel_list.object(loc_ind).printf(dbgfile)
dbgfile.printf("];\n")
dbgfile.printf("%% now x:\nxdat_rel%d = [",exind+1)
xdat_rel_list.object(loc_ind).printf(dbgfile)
dbgfile.printf("];\n")
dbgfile.printf("%% here xpeak = %g;\n",xpeak)
}
if( (shape_scale > 0) && (apwin.x[1]-apwin.x[0] > 0 )) {
//printf("MULT model win, Match ALL: before efun_yfit, time indices are [%g, %g]\n",$o2.x[mnind],$o2.x[mxind])
ey += efun_yfit($o1,$o2,i,exind,mnind,mxind,0)
//printf("(2) Return from this AP, ey now %g.\n",ey)
}
if( EFUN_DBG && VERBOSE ) {
dbgfile.printf("%% first pass\n")
dbgfile.printf("peak(%d) = %g; xpeak(%d) = %g;\n",i+1,peak,i-mod_APleft+exp_APleft+1,xpeak)
dbgfile.printf("ey(%d) = %g;\n",i+1,ey)
}
} else {
// loop over all chosen experimental spikes
for exind = exp_APleft, exp_APright {
loc_ind = exind - exp_APleft
//printf("Chose specific APs. Compare model AP %d to exp AP %d",i,exind)
xpeak = xdat.x[ydat.nextpeak(xyAPtimes.x[exind])]
if( EFUN_DBG && VERBOSE ) {
dbgfile.printf("%% before efun_yfit\n")
dbgfile.printf("mnind(%d) = %d; mxind(%d) = %d;\n",i+1,mnind,i+1,mxind)
dbgfile.printf("exind(%d) = %d;\n",i+1,exind)
print "Match model AP ",i," to experiment AP ",exind
if( i>0 ) {
dbgfile.printf("t_th(%d)=%g; t_pk(%d)=%g; t_apk(%d)=%g;\n",i+1,$o2.x[modelAPtimes.x[i-1]],i+1,$o2.x[$o1.nextpeak(modelAPtimes.x[i-1])],i+1,$o2.x[mnind])
}
dbgfile.printf("%% y direction: Y %d X %d\n",ydat_rel_list.object(loc_ind).size,\
xdat_rel_list.object(loc_ind).size)
dbgfile.printf("%% now y:\n ydat_rel%d = [",exind+1)
ydat_rel_list.object(loc_ind).printf(dbgfile)
dbgfile.printf("];\n")
dbgfile.printf("%% now x:\nxdat_rel%d = [",exind+1)
xdat_rel_list.object(loc_ind).printf(dbgfile)
dbgfile.printf("];\n")
}
if( (shape_scale > 0) && (apwin.x[1]-apwin.x[0] > 0 )) {
//printf("MULT model win, match AP %d: before efun_yfit, time indices are [%g, %g]\n",i,$o2.x[mind],$o2.x[mxind])
ey += efun_yfit($o1,$o2,i,exind,mnind,mxind,0)
//printf("(1) Return from this AP.\n")
}
if( EFUN_DBG && VERBOSE ) {
dbgfile.printf("%% second pass\n")
dbgfile.printf("peak(%d) = %g; xpeak(%d) = %g;\n",i+1,peak,i-mod_APleft+exp_APleft+1,xpeak)
dbgfile.printf("ey(%d) = %g;\n",i+1,ey)
}
} // end of exp spike loop
} // end else
} // end of modelAP loop
if( shape_scale > 0 ) {
//printf("Now scaling the shape error %g by %g / %d\n",ey,shape_scale,n_modWin)
ey *= shape_scale / n_modWin
//ey *= shape_scale / n_modAP
} else { ey = 0 }
//if( shape_scale > 0 ) { ey *= shape_scale / n_modAP } else { ey = 0 }
if( EFUN_DBG && VERBOSE ) {
dbgfile.printf("ey_mean = %g;\n",ey)
}
}
//printf("Returning from calc_shape_error()\n")
}
// Penalize ISI's which are less than a user-specified threshold brst_thold,
// by a user specified amount brst_pnlty
//
// This will penalize doublets and bursting.
//
func efun_penburst() { local etot, eisi, cnt
etot = 0
cnt = 0
for i = mod_APleft, mod_APright-1 {
eisi = $o1.x[modelAPtimes.x[i+1]] - $o1.x[modelAPtimes.x[i]]
if( eisi <= brst_thold ) {
etot += brst_pnlty
cnt += 1
if( EFUN_DBG && VERBOSE ) {
if( i >= 0 ) {
dbgfile.printf("eISI(%d)=%g;\n",i+1,brst_pnlty)
} else {
dbgfile.printf("fprintf('Trouble! Found ISI index %d\\n');\n",i+1)
dbgfile.printf("eISI(1)=%g;\n",brst_pnlty)
}
}
} else {
if( EFUN_DBG && VERBOSE ) {
if( i >= 0 ) {
dbgfile.printf("eISI(%d)=0;\n",i+1)
} else {
dbgfile.printf("fprintf('Trouble! Found ISI index %d\n');\n",i+1)
dbgfile.printf("eISI(1)=0;\n")
}
}
}
}
if( EFUN_DBG ) {
//dbgfile.printf("str(%d) = {'ISI burst penalty (%d times) %g'};\n", nstr,cnt,etot)
//nstr += 1
dbgfile.printf("lbl(%d) = {'ISI burst penalty (%d times) %g'};\n", nlbl,cnt,etot)
nlbl += 1
printf("ISI burst penalty (%d times) %g\n", cnt,etot)
}
return etot
}
/***********************************************************************
Penalize the difference between the target and model delay to first spike.
Let
S0 = first target spike between the error bounds
s0 = first model spike between the error bounds
(this equals bounds.x[1]-bounds.x[0] if there is no model spike)
tISI = the mean target ISI during the error period,
i.e. reciprocal of mean FR (in ms).
Then the first spike penalty is given by
(S0 - s0)^2 / tISI^2,
that is, the mean squared difference between the first spike times,
normalized as a multiple of the mean target ISI.
This penalty is scaled by the factor DELAY_PNLTY.
INPUT $o1 y values of model output
$o2 t values of model output
***********************************************************************/
func efun_pendelay() { local errval, mISI, tspk, mspk
if( n_expAP < 1 ) return 0
//
// determine first spike times for target and model
//
tspk = xdat.x[xyAPtimes.x[exp_APleft]]
if( n_modAP == 0 ) {
mspk = boundary.x[1]
} else {
mspk = $o2.x[modelAPtimes.x[mod_APleft]]
}
//printf("tspk %g - mspk %g = %g\n",tspk,mspk,tspk-mspk)
//
// determine the mean ISI (in ms) during the error bounds
//
if( n_expAP == 1 ) {
mISI = boundary.x[1] - tspk
} else {
/****
printf("calculating the mean FR between spike %d to %d\n",exp_APleft,exp_APright)
printf("all target spike times : ")
xdat.ind(xyAPtimes).printf
printf("all of datFR = ")
datFR.printf()
printf("we care about indices from %d to %d: ",exp_APleft,exp_APright-1)
datFR.printf("%g ",exp_APleft,exp_APright-1)
****/
if( exp_APright-exp_APleft > 2 ) {
mISI = 1000 / datFR.mean(exp_APleft,exp_APright-1)
} else {
mISI = 1000 / datFR.x[exp_APleft]
}
//printf("found mean ISI %g\n",mISI)
}
errval = ( (tspk - mspk) / mISI ) ^ 2
//printf("found error %g; tspk %g, mspk %g, mISI %g, times %g\n",errval,tspk,mspk,mISI,delay_pnlty)
if( EFUN_DBG ) {
//dbgfile.printf("str(%d) = {'First spike delay penalty (times %g) %g'};\n", nstr,delay_pnlty,errval)
//nstr += 1
dbgfile.printf("lbl(%d) = {'First spike delay penalty (times %g) %g'};\n", nlbl,delay_pnlty,errval)
nlbl += 1
printf("First spike delay penalty (times %g) %g\n", delay_pnlty, errval)
if( VERBOSE ) {
dbgfile.printf("e_dly = %g;\n",errval)
}
}
return delay_pnlty * errval
}
// the least squares error function in y
func efun_yfit() { local e, xtmp_modind, tmp_expind, e_old, mnind, mxind, loc_ind, eval, loc_mdl,idx
// only use e_old for debug
// INPUT:
// $o1 y-values of model output (e.g soma.v(.5) )
// $o2 t-values of model output (time, ms)
mod_ind = $3
exp_ind = $4
loc_ind = exp_ind - exp_APleft
loc_mdl = mod_ind - mod_APleft
mnind = $5
mxind = $6
objref xtmp, ytmp
// Shift voltage trace of target data to align target and model AP troughs along the V axis.
if( VERBOSE ) {
printf("loc_ind = %d; loc_mdl = %d; tgt size = %d %d; mod size = %d, %d\n",\
loc_ind,loc_mdl,tgtAPtrough.size,tgtAPtrough_V.size,modelAPtrough.size,\
modelAPtrough_V.size)
//printf("tgt troughs : ")
//tgtAPtrough_V.printf
//printf("model troughs : ")
//modelAPtrough_V.printf
}
xtmp = xdat_rel_list.object(loc_ind).c
if( modelAPtrough.size() > 0 && modelAPtrough.size() > loc_mdl ) {
if( VERBOSE ) {
printf("\tTarget AP %d, add - %g + %g = %g\n",loc_ind,\
tgtAPtrough_V.x[loc_ind],modelAPtrough_V.x[loc_mdl],\
-1*tgtAPtrough_V.x[loc_ind]+modelAPtrough_V.x[loc_mdl])
}
ytmp = ydat_rel_list.object(loc_ind).c.add(-1*tgtAPtrough_V.x[loc_ind]+modelAPtrough_V.x[loc_mdl])
} else {
if( VERBOSE ) {
printf("No model APs, so no V shift.\n")
}
ytmp = ydat_rel_list.object(loc_ind).c
}
//printf("In EFUN_YFIT(): the model data go from t = %g to %g\n",$o2.x[mnind],$o2.x[mxind])
e = 0
if( EFUN_DBG ) {
dbgfile.printf("%% e = %g; peak=%g; xpeak=%g;\n",e,peak,xpeak)
}
//printf("mod_ind = %d, n_mod %d; modsz %d\n",mod_ind,n_modAP,modelAPtimes.size())
if( n_modAP > 0 ) {
tmp_modind = $o1.nextpeak(modelAPtimes.x[mod_ind])
} else {
tmp_modind = $o1.firstmax
if( tmp_modind == 0 ) tmp_modind += 1
}
//printf("set tmp_modind to %g\n",tmp_modind)
if ( xdat_rel_list.object(loc_ind).size > 0 ) {
if( EFUN_DBG ) {
e_old = e
tmp_expind = ytmp.firstpeak
//printf("tmp_expind = %g\n",tmp_expind)
dbgfile.printf("%% Before yfitness, model peak = %g (%g, %g - %g - %g), ",\
tmp_modind,$o2.x[tmp_modind],$o1.x[tmp_modind-1],$o1.x[tmp_modind],$o1.x[tmp_modind+1])
dbgfile.printf("exp peak = %g (%g, %g - %g - %g )\n",tmp_expind,xtmp.x[tmp_expind],\
ytmp.x[tmp_expind]-1,ytmp.x[tmp_expind],ytmp.x[tmp_expind+1])
}
//printf("Calling ywnscl_fitness MDx.sz %d, MDy.sz %d; TGx.sz %d TGy.sz %d\n",$o1.size,$o2.size,ytmp.size,xtmp.size)
//printf("\tmodel pkind = %d [%g, %g], exp pk ind %d [%g, %g];",tmp_modind,$o2.x[tmp_modind],$o1.x[tmp_modind],ytmp.firstpeak,xtmp.x[ytmp.firstpeak],ytmp.x[ytmp.firstpeak])
//printf(" time indices [%g, %g]\n",$o2.x[mnind],$o2.x[mxind])
eval = $o1.ywnscl_fitness($o2, peak, ytmp, xtmp,ytmp.firstpeak,tmp_modind,mnind,mxind)
if( eval >= 0 ) {
if( EFUN_DBG ) printf("\tModel AP %d (time %.2f) shape error = %g\n",loc_mdl,$o2.x[modelAPtimes.x[loc_mdl]],eval)
e += eval
} else {
printf("Error! No model window was found; shape error for this AP is -1.\n")
}
if( EFUN_DBG ) {
dbgfile.printf("%% yfitness %s Y %d X %d\n",$o2.label,peak,ytmp.size,xtmp.size)
dbgfile.printf("ey_scale(%d) = %g;\n",mod_ind+1,e-e_old)
dbgfile.printf("tgt_AP{%d} = [\n",loc_ind+1)
for idx = 0, ytmp.size()-1 {
dbgfile.printf("\t%g\t%g\n",xtmp.x[idx],ytmp.x[idx])
}
dbgfile.printf("];\n")
}
}
if( EFUN_DBG ) { dbgfile.printf("e_1AP(%d)=%g;\n",mod_ind+1,e) }
return e
}
func efun_expfit() { local mnind, mxind, lastfr, err, i, tauerr, Aerr
// $o1 voltage
// $o2 time vec
// $3 leftmost AP from exp data
// $4 rightmost AP from exp data
// $5 leftmost AP from model
// $6 rightmost AP from model
// $7 weight for intercept error
// $8 weight for time constant error
// $9 time offset (start of region for exp'l fitting)
ntag = ""
mnind = $3
mxind = $4 - 1
//
// First, perform a linear fit to the target data.
//
//printf("Start of target fit: mnind = %d, mxind = %d.\n",mnind,mxind)
if( mxind < 0 || mnind < 0 || (mxind-mnind < 2) ) {
// Not enough experimental spikes to fit a line
if( EFUN_DBG ) printf("Not enough exp data for fit.\n")
pen_silent = 0
exp_intcpt = exp_slope = 0
expsa = expsb = expchi = 0
} else {
// do the linear fit
ptra = new Pointer(&exp_intcpt)
ptrb = new Pointer(&exp_slope)
ptrsiga = new Pointer(&expsa)
ptrsigb = new Pointer(&expsb)
ptrchi = new Pointer(&expchi)
ytmp = new Vector()
ytmp = datFR.c
for i = 0, ytmp.size-1 {
ytmp.x[i] = log(datFR.x[i])
}
if( EFUN_DBG ) {
printf("TARGET DATA: Exponential fit, Fitting times (minus %g ms): \n",$9)
xdat.ind(xyAPtimes.c(mnind,mxind)).printf
printf("And firing rates: \n")
datFR.c(mnind,mxind).printf
printf("Converted to log:\n")
ytmp.c(mnind,mxind).printf
printf("Done converting target data, now fit.\n")
}
linfit(xdat.ind(xyAPtimes).sub($9),ytmp,mnind,mxind,ptra,ptrb,ptrsiga,ptrsigb,ptrchi)
//printf("Done the fit, now convert parameters.\n")
//printf("Lin fit complete, %g and %g\n",ptra.val,ptrb.val)
exp_intcpt = exp(ptra.val)
if( exp_slope == 0 || ptrb.val == 0 ) {
exp_slope = -999999.99
} else {
exp_slope = -1 / ptrb.val
}
expsa = ptrsiga.val
expsb = ptrsigb.val
expchi = ptrchi.val
// target data *is* firing APs; add extra penalty for non-firing models
pen_silent = 1
//printf("Finished doing the fit; now check it.\n")
}
//
// Now, perform a linear fit to the model data.
//
mmnind = $5
mmxind = $6 - 1
//printf("Now for model data: Before starting fit.\n")
if( mmxind < 0 || mmnind < 0 || mmxind-mmnind < 2 ) {
if( EFUN_DBG ) printf("Not enough model data for fit.\n")
modsa = modsb = modchi = 0
mod_slope = -1e5
mod_intcpt = 0
//mod_slope = -999999.99 // 7/8/11: penalty too high
//mod_intcpt = 1
} else {
// do the linear fit
ptra = new Pointer(&mod_intcpt)
ptrb = new Pointer(&mod_slope)
ptrsiga = new Pointer(&modsa)
ptrsigb = new Pointer(&modsb)
ptrchi = new Pointer(&modchi)
ytmp = new Vector()
ytmp = modFR.c
for i = 0, ytmp.size-1 {
ytmp.x[i] = log(modFR.x[i])
}
if( EFUN_DBG ) {
printf("MODEL DATA: Exponential fit, Fitting times (minus %g ms): \n",$9)
$o2.ind(modelAPtimes.c(mmnind,mmxind)).printf
printf("And firing rates: \n")
modFR.c(mmnind,mmxind).printf
printf("Converted to log:\n")
ytmp.c(mmnind,mmxind).printf
printf("Copied data, now do lin fit.\n")
}
linfit($o2.ind(modelAPtimes).sub($9),ytmp,mmnind,mmxind,ptra,ptrb,ptrsiga,ptrsigb,ptrchi)
//printf("Lin fit complete, %g and %g\n",ptra.val,ptrb.val)
mod_intcpt = exp(ptra.val)
if( mod_slope == 0 || ptrb.val == 0 ) {
mod_slope = -999999.99
mod_slope = -999999.99
} else {
mod_slope = -1 / ptrb.val
}
modsa = ptrsiga.val
modsb = ptrsigb.val
modchi = ptrchi.val
}
//printf("Done fitting model data.\n")
// use the exponential fit to estimate firing rates at each of the TARGET AP times.
modFit = new Vector()
for i = mnind, mxind {
modFit.append(mod_intcpt*exp(-1*xdat.x[xyAPtimes.x[i]]/mod_slope))
}
//printf("Target fit estimated.\n")
if( EFUN_DBG ) {
printf("Fitting an exponential to the FR data\n")
printf("Experiment LSFIT tau %.8f +- %.8f\n",exp_slope,expsb)
printf("Experiment A0 %.8f +- %.8f\n",exp_intcpt,expsa)
printf("Experiment chi_sq %.8f\n\n",expchi)
//sprintf(title,"y = %.2f + %.2f * x\0",a,b)
printf("Model LSFIT tau %.8f +- %.8f\n",mod_slope,modsb)
printf("Model A0 %.8f +- %.8f\n",mod_intcpt,modsa)
printf("Model chi_sq %.8f\n\n",modchi)
nstr2 = 1
dbgfile.printf("str2(%d) = {'Exponential FR fit'};\n",nstr2)
nstr2 += 1
dbgfile.printf("str2(%d) = {'Exper LSFIT tau %g +- %g'};\n",nstr2,exp_slope,expsb)
nstr2 += 1
dbgfile.printf("str2(%d) = {'Exper A0 %.8f +- %.8f'};\n",nstr2,exp_intcpt,expsa)
nstr2 += 1
dbgfile.printf("str2(%d) = {'Exper chi_sq %.8f'};\n\n",nstr2,expchi)
nstr2 += 1
dbgfile.printf("str2(%d) = {'Model LSFIT tau %.8f +- %.8f'};\n",nstr2,mod_slope,modsb)
nstr2 += 1
dbgfile.printf("str2(%d) = {'Model A0 %.8f +- %.8f'};\n",nstr2,mod_intcpt,modsa)
nstr2 += 1
dbgfile.printf("str2(%d) = {'Model chi_sq %.8f'};\n\n",nstr2,modchi)
nstr2 += 1
}
// Add the error between the two components of the fit: The initial offset, and the time constant.
err = 0
if( EFUN_DBG ) {
}
Aerr = $7 * abs(exp_intcpt - mod_intcpt)
tauerr = $8 * abs(exp_slope-mod_slope)
err = Aerr + tauerr
//printf("err %g = Aerr %g + tauerr %g = %g * abs(%g - %g) + %g * abs(%g - %g) \n",err,Aerr,tauerr,intcp_scale,exp_intcpt,mod_intcpt,slope_scale,exp_slope,mod_slope)
if( EFUN_DBG ) {
dbgfile.printf("str2(%d) = {'Error in initial FR = %.2f * %.2f = %2f'};\n\n",nstr2, $7,abs(exp_intcpt-mod_intcpt),Aerr)
dbgfile.printf("str2(%d) = {'Error in time constants = %.2f * %.2f = %2f'};\n\n",nstr2, $8,abs(exp_slope-mod_slope),tauerr)
print "Error in initial FR = ", $7, " * ", abs(exp_intcpt-mod_intcpt)," = ",Aerr
print "Error in time constants = ", $8, " * ", abs(exp_slope-mod_slope)," = ",tauerr
print "Total exponential fit error = ",err
nstr2 += 1
}
return err
}
/*
* From NUMERICAL RECIPES IN C
*
* Given a set of points x[strt ... end], y[strt ... end], with standard
* deviations sig[strt ... end], fit them to
*
* y = a + bx (!!!!!!!!!)
*
* by minimizing chi_sq.
* Returned are a,b, siga, sigb, chi_sq (chi2), and the goodness of fit
* probability q (that the fit would have chi_sq this large or larger).
* If mwt = 0 * on input, then the standard deviations are assumed
* unavailable, q is returned as 1.0, and the normalization of chi2
* is to unit standard deviation on all points.
*/
func SQR() { return $1*$1 }
proc linfit() { local strt, endfit, a,b,siga,sigb,chi2, i,wt,tfit,sxoss,ss, sigdat, sx, sy, st2, ndata
// x = $o1
// y = $o2
strt = $3
endfit = $4
ndata = endfit-strt+1
if( EFUN_DBG ) {
printf("start of linfit(): start at %d, end at %d\n",strt,endfit)
$o1.printf
$o2.printf
printf("1st %d=(%g,%g), last %d=(%g,%g)\n",strt,$o1.x[strt],$o2.x[strt],endfit,$o1.x[endfit],$o2.x[endfit])
}
sx = 0.0
sy = 0.0
st2 = 0.0
b = 0.0
for i = strt, endfit {
sx += $o1.x[i] sy += $o2.x[i]
}
ss = endfit+1-strt
//printf("sx=%g\tsy=%g\tss=%g\n",sx,sy,ss)
sxoss = sx/ss
//printf("sxoss = %g\n",sxoss)
for i = strt, endfit {
tfit = $o1.x[i] - sxoss
st2 += tfit*tfit
b += tfit*$o2.x[i]
}
//printf("st2=%g\tb=%g\n",st2,b)
b /= st2
a = (sy-sx*b)/ss
//printf("a=%g\n",a)
siga = sqrt((1.0 + sx*sx/(ss*st2))/ss)
sigb = sqrt(1.0/st2)
chi2 = 0.0
for i = strt, endfit { chi2 += SQR($o2.x[i]-a-b*$o1.x[i]) }
if( ndata > 2 ) { sigdat = sqrt(chi2/ndata)
} else { sigdat = 0 }
siga *= sigdat
sigb *= sigdat
if( EFUN_DBG ) {
printf("LSFIT slope %.8f +- %.8f\n",b,sigb)
printf(" int %.8f +- %.8f\n",a,siga)
printf(" chi_sq %.8f\n\n",chi2)
//sprintf(title,"y = %.2f + %.2f * x\0",a,b)
}
$o5.val = a
$o6.val = b
$o7.val = siga
$o8.val = sigb
$o9.val = chi2
}
/*******************************
Return RMS residual between target instantaneous firing rates, and
firing rates at those same target times, as estimated by a linear fit
of the model instantaneous firing rates.
*******************************/
func efun_linfit() { local mnind, mxind, err, pen_silent, mmnind, mmxind
// $1 start place
// $2 end place
ntag = ""
mnind = exp_APleft
mxind = exp_APright - 1
//
// First, perform a linear fit to the target data.
//
if( mxind < 0 || mnind < 0 || (mxind-mnind < 2) ) {
// Not enough experimental spikes to fit a line
pen_silent = 0
exp_intcpt = exp_slope = 0
expsa = expsb = expchi = 0
} else {
// do the linear fit
ptra = new Pointer(&exp_intcpt)
ptrb = new Pointer(&exp_slope)
ptrsiga = new Pointer(&expsa)
ptrsigb = new Pointer(&expsb)
ptrchi = new Pointer(&expchi)
//printf("TARGET DATA: Linear fit, Fitting times: \n")
//xdat.ind(xyAPtimes).printf
//printf("And firing rates: \n")
//datFR.printf
linfit(xdat.ind(xyAPtimes),datFR,mnind,mxind,ptra,ptrb,ptrsiga,ptrsigb,ptrchi)
exp_intcpt = ptra.val
exp_slope = ptrb.val
expsa = ptrsiga.val
expsb = ptrsigb.val
expchi = ptrchi.val
// target data *is* firing APs; add extra penalty for non-firing models
pen_silent = 1
}
datFit = new Vector()
datFit = xdat.ind(xyAPtimes).c(mnind,mxind)
// datFit = xyAPtimes.c(mnind,mxind)
datFit.mul(exp_slope)
datFit.add(exp_intcpt)
//
// Now, perform a linear fit to the model data.
//
mmnind = mod_APleft
mmxind = mod_APright - 1
if( mmxind < 0 || mmnind < 0 || mmxind-mmnind < 2 ) {
mod_intcpt = mod_slope = modsa = modsb = modchi = 0
} else {
// do the linear fit
ptra = new Pointer(&mod_intcpt)
ptrb = new Pointer(&mod_slope)
ptrsiga = new Pointer(&modsa)
ptrsigb = new Pointer(&modsb)
ptrchi = new Pointer(&modchi)
//printf("MODEL DATA: Linear fit, Fitting times: \n")
//$o2.ind(modelAPtimes).printf
//printf("And firing rates: \n")
//modFR.printf
linfit($o2.ind(modelAPtimes),modFR,mmnind,mmxind,ptra,ptrb,ptrsiga,ptrsigb,ptrchi)
mod_intcpt = ptra.val
mod_slope = ptrb.val
modsa = ptrsiga.val
modsb = ptrsigb.val
modchi = ptrchi.val
}
// use the linear fit to estimate firing rates at each of the TARGET AP times.
modFit = new Vector()
modFit = xdat.ind(xyAPtimes).c(mnind,mxind)
//modFit = xyAPtimes.c(mnind,mxind)
modFit.mul(mod_slope)
modFit.add(mod_intcpt)
if( EFUN_DBG ) {
printf("Fitting a straight line to the FR data\n")
printf("Experiment LSFIT slope %.8f +- %.8f\n",exp_slope,expsb)
printf("Experiment int %.8f +- %.8f\n",exp_intcpt,expsa)
printf("Experiment chi_sq %.8f\n\n",expchi)
//sprintf(title,"y = %.2f + %.2f * x\0",a,b)
printf("Model LSFIT slope %.8f +- %.8f\n",mod_slope,modsb)
printf("Model int %.8f +- %.8f\n",mod_intcpt,modsa)
printf("Model chi_sq %.8f\n\n",modchi)
nstr2 = 1
dbgfile.printf("str2(%d) = {'Straight line FR fit'};\n",nstr2)
nstr2 += 1
dbgfile.printf("str2(%d) = {'Exper LSFIT slope %g +- %g'};\n",nstr2,exp_slope,expsb)
nstr2 += 1
dbgfile.printf("str2(%d) = {'Exper int %.8f +- %.8f'};\n",nstr2,exp_intcpt,expsa)
nstr2 += 1
dbgfile.printf("str2(%d) = {'Exper chi_sq %.8f'};\n\n",nstr2,expchi)
nstr2 += 1
dbgfile.printf("str2(%d) = {'Model LSFIT slope %.8f +- %.8f'};\n",nstr2,mod_slope,modsb)
nstr2 += 1
dbgfile.printf("str2(%d) = {'Model int %.8f +- %.8f'};\n",nstr2,mod_intcpt,modsa)
nstr2 += 1
dbgfile.printf("str2(%d) = {'Model chi_sq %.8f'};\n\n",nstr2,modchi)
nstr2 += 1
}
//
// Calculate the RMS error between the model fit and the target fit,
// evaluated at each target AP time.
//
err = 0
for i = 0, modFit.size-1 {
if( EFUN_DBG ) {
printf("Err was %f; Add SQR( E %f - M %f)\n",err,datFit.x[i],modFit.x[i])
dbgfile.printf("str2(%d) = {'Exp AP %d, error %f'};\n",nstr2,i,SQR(datFit.x[i]-modFit.x[i]))
nstr2 += 1
}
err += SQR(datFit.x[i]-modFit.x[i])
}
if( modFit.size > 0 ) {
err = sqrt( err / modFit.size )
}
//
// Scale the error further as an extra penalty for models which do not
// fire. The scaling factor is defined by the user in scale_factors() below.
// This extra penalty is only applied when the target does exhibit APs.
//
if( pen_silent && (mmxind-mmnind < 2) ) {
err *= nf_frscale
sprint(ntag,"(silent - extra scaled by %f)",nf_frscale)
}
if( EFUN_DBG ) {
dbgfile.printf("str2(%d) = {'Total linear fit error = %.2f * %.2f %s'};\n\n",nstr2,err,tot_frscale,ntag)
print "Total linear fit error = ",err," * ",tot_frscale, ntag
nstr2 += 1
}
return err
}
proc init() {local i
print "initializing APFRfit_2exp, last modified 7 Jul 11"
EFUN_DBG = 0
VERBOSE = 0
GBAR_SMRY = 0
sprint(tag, "%s", this)
sscanf(tag, "%[^[]", tag)
use_x = 0
use_gui = 0
have_data = 0
//mintstop_ = 0
xdat = new Vector(0)
ydat = new Vector(0)
idx_rel_list = new List()
xdat_rel_list = new List()
ydat_rel_list = new List()
boundary = new Vector(2)
bound_2exp = 0
apwin = new Vector(2)
apwin.x[0] = apwin.x[1] = 0
outname = "err_data"
xyAPtimes = new Vector(0)
xpeak = 0
t_hat = 10
pow1 = 0.4
pow2 = 1
frscale=1
cvscale=1
time_scale = miss_scale = 0.1
shape_scale = 5
frm_fac = 10
frs_fac = 5
cvm_fac = 100
cvs_fac = 1
frcv_int = 1000
sprint(scalelabel, "scale shape=%g FR=%g CV=%g",\
shape_scale,frscale,cvscale)
brst_pnlty = brst_thold = 0
delay_pnlty = 1
nf_frscale = 1
allList = new SectionList()
// these are to measure AHP
MS_AHP = 0
neg_p = 0
trgh_sz = 0
trgh_ps = 0
thr_ac = 0
//print "\tdone initializing AP_instFR_Fitness"
/******************
NOTE: From build():
the following functions are needed to set up the data properly
g.menu_tool("Adjust", "adjust_region")
xbutton("Paste data from Clipboard", "clipboard_data()")
xbutton("Select overall time window", "region_panel()") *** use set_r()
xbutton("Shape Error: AP window and scale factor", "set_apwin()")
xbutton("Firing Rate Error: Scale factors for function fit", "scale_factors()")
xbutton("Firing Rate Error: Overall scale factors", "FR_GenScale_factors()")
xbutton("Error for Bursting","set_burstpen()")
xbutton("Error for Delay in Firing","set_delaypen()")
xbutton("Set info for file output", "output_panel()")
********************/
}
proc clone() {
$o1 = new APFRfit_2exp()
$o1.have_data = have_data
if (have_data) {
$o1.set_data(xdat, ydat)
}
$o1.boundary = boundary.c
$o1.bound_2exp = bound_2exp
$o1.apwin = apwin.c
$o1.xyAPtimes = xyAPtimes.c
$o1.set_modelx(idx_rel_list)
$o1.t_hat = t_hat
$o1.pow1 = pow1
$o1.pow2 = pow2
$o1.shape_scale = shape_scale
$o1.frscale = frscale
$o1.cvscale = cvscale
$o1.frm_fac = frm_fac
$o1.frs_fac = frs_fac
$o1.cvm_fac = cvm_fac
$o1.cvs_fac = cvs_fac
$o1.brst_thold = brst_thold
$o1.brst_pnlty = brst_pnlty
$o1.delay_pnlty = delay_pnlty
sprint(scalelabel, "scale shape=%g FR=%g CV=%g",\
shape_scale,frscale,cvscale)
}
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[xdat.size-1]
bound_2exp = boundary.mean
have_data = 1
xpeak = xdat.x[ydat.max_ind]
}
proc set_modelx() {local i
objref idx_rel_list
idx_rel_list = new List()
for i = exp_APleft, exp_APright {
idx_rel_list.append($o1.object(i))
}
set_allxy()
}
func mintstop() {
return boundary.x[1]
}
proc setxy() {local i, ap_ind, glob_ind
// $1 index of AP in the master list
// $2 index of AP, within the boundary.
if( PASTE_DBG ) {
print "SETXY START:"
for i = 0, ydat_rel_list.count-1 {
print "\ty[",i,"] peak ",ydat_rel_list.object(i).x[ydat_rel_list.object(i).firstpeak]
}
dbgfile = new File()
dbgfile.aopen("paste_curve.m")
dbgfile.printf("%% in setxy() \n")
}
mintstop_ = 0
glob_ind = $1
ap_ind = $2
if( PASTE_DBG ) { printf("SETXY: appending AP %d global %d\n",ap_ind,glob_ind) }
xpeak = xdat.x[ydat.nextpeak(xyAPtimes.x[glob_ind])]
if( PASTE_DBG ) {
print "\t",ap_ind,".\tIn setxy, index size ",idx_rel_list.object(ap_ind).size
dbgfile.printf("idx(:,%d) = [",ap_ind+1)
idx_rel_list.object(ap_ind).printf(dbgfile)
print "xdat size ",xdat.size
dbgfile.printf("];\n")
}
xtmp.index(xdat, idx_rel_list.object(ap_ind))
xtmp.sub(xpeak)
xdat_rel_list.append(xtmp.c)
if( PASTE_DBG ) {
print "xdat size ",xdat_rel_list.object(ap_ind).size
dbgfile.printf("%% xdat_rel\nxdat_rel_sub(:,%d) = [",ap_ind+1)
xdat_rel_list.object(ap_ind).printf(dbgfile)
dbgfile.printf("];\n")
}
ytmp.index(ydat, idx_rel_list.object(ap_ind))
ydat_rel_list.append(ytmp.c)
if( PASTE_DBG ) {
print "ydat size ",ydat_rel_list.object(ap_ind).size
dbgfile.printf("%% ydat_rel\nydat_rel_set(:,%d) = [",ap_ind+1)
ydat_rel_list.object(ap_ind).printf(dbgfile)
dbgfile.printf("];\n")
}
if ( xtmp.x[xtmp.size-1] < boundary.x[1] && \
mintstop_ < xtmp.x[xtmp.size-1]) {
mintstop_ = xtmp.x[xtmp.size-1]
}
if( PASTE_DBG ) {
dbgfile.printf("mintstop = %g;\n",mintstop_)
dbgfile.close()
print "SETXY END:"
for i = 0, ydat_rel_list.count-1 {
print "\ty[",i,"] peak ",ydat_rel_list.object(i).x[ydat_rel_list.object(i).firstpeak]
}
}
}
proc set_allxy() { local iloop
//print "Here in set_allxy()"
left_ptr = new Pointer(&exp_APleft)
mid_ptr = new Pointer(&exp_APmid)
right_ptr = new Pointer(&exp_APright)
n_ptr = new Pointer(&n_expAP)
nw_ptr = new Pointer(&n_expWin)
find_APtimes(xdat,ydat,xyAPtimes,"Exp AP",1,left_ptr,right_ptr,n_ptr,nw_ptr,mid_ptr)
set_relative_idx()
for iloop = exp_APleft, exp_APright {
//print "\t",iloop,". setting xy"
setxy(iloop,iloop-exp_APleft)
}
//print "end of set_allxy(), mintstop = ",mintstop_
}
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")
}
set_allxy()
}
proc build() {
/**** the following functions are needed to set up the data properly
g.menu_tool("Adjust", "adjust_region")
xbutton("Paste data from Clipboard", "clipboard_data()")
xbutton("Select overall time window", "region_panel()")
xbutton("Shape Error: AP window and scale factor", "set_apwin()")
xcheckbox("Match all APs",&match_all,"match_allAP()")
xbutton("Firing Rate Error: Class of functions for FR fit","choose_function_panel()")
xbutton("Firing Rate Error: Scale factors for function fit", "scale_factors()")
xbutton("Firing Rate Error: Overall scale factors", "FR_GenScale_factors()")
xbutton("Error for Bursting","set_burstpen()")
xbutton("Error for Delay in Firing","set_delaypen()")
xbutton("Set info for file output", "output_panel()")
****/
}
/**************************************************************
modified from APtv_instFR_Fitness factors() proc, which needs
the GUI.
input $1 t_hat
$2 pow1
$3 pow2
**************************************************************/
proc factors() {local th, pw, pw2
//"Enter space separated t_hat, pow1, pow2 parameters"
t_hat = $1
pow1 = $2
pow2 = $3
}
proc FR_GenScale_factors() {local sls, ics
//"Enter scale factors for (1) TOTAL WEIGHT of firing rate error and
// (2) EXTRA PENALTY FOR NON-FIRING MODELS"
tot_frscale = $1
nf_frscale = $2
}
proc exp_scale_factors() {local sls, ics
//"Enter scale factors for the (1) TIME CONSTANT and
// (2) AMPLITUDE of the exponential fit of FR.",tmpstr)){
slope_scale = $1
intcp_scale = $2
slope2_scale = $3
intcp2_scale = $4
nf_frscale = $2
}
proc scale_factors() {local shs, bs
exp_scale_factors()
}
proc set_burstpen() {local th, pw
//"For each ISI < BURST_THOLD (ms), add BURST_PENALTY.\n
//Enter space separated BURST_THOLD and BURST_PENALTY parameters")
brst_thold = $1
brst_pnlty = $2
}
proc set_delaypen() { local pn
//sprint(tmpname,"Enter scale factor for First Spike Delay Penalty")
delay_pnlty = $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() {
// 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("AP_instFR_Fitness xdat ydat boundary apwin (lines=%d) ",\
4 + 2*xdat.size + 2 +2 )
$o1.printf(" %g %g %g %g %g %g %g\n",\
shape_scale,tot_frscale, nf_frscale, slope_scale,\
brst_thold, brst_pnlty, delay_pnlty)
$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\n%g\n", boundary.x[0],bound_2exp,boundary.x[1])
$o1.printf("%d\n", apwin.size)
apwin.printf($o1)
}
proc rfile() {local i, n
shape_scale = $o1.scanvar
tot_frscale = $o1.scanvar
nf_frscale = $o1.scanvar
intcp_scale = nf_frscale // for exponential fits - reuse this variable but give it a more meaningful name.
slope_scale = $o1.scanvar
printf("INIT intcp_scale %g, slope_scale %g\n",intcp_scale,slope_scale)
intcp2_scale = $o1.scanvar
slope2_scale = $o1.scanvar
brst_thold = $o1.scanvar
brst_pnlty = $o1.scanvar
delay_pnlty = $o1.scanvar
sprint(scalelabel, "shape=%g",shape_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
apwin.resize(n)
apwin.scanf($o1, n)
//print "apwin = ",apwin.x[0],apwin.x[1]
//print "here in rfile()"
set_relative_idx()
set_allxy()
//print "rfile done set_allxy() "
//printf("rfile(): \n")
//print shape_scale, tot_frscale, nf_frscale,slope_scale,brst_thold,brst_pnlty,delay_pnlty
//apwin.printf
}
/****************************************************
Count APs which occur during the specified time boundary.
****************************************************/
proc find_APtimes() { local i, firing, thresh, idx
// $o1 x component, $o2 y component of a voltage trace
// $o3 stores AP times
// $s4 stores descriptive string
// $5 1 or 0, print result to screen?
// $o6, points to AP after left boundary
// $o7, points to AP before right boundary
// $o8, points to count of APs within boundary
// $o9, pointer to the number of AP shape windows within
// boundary (which may differ from the number of APs)
// $o10, points to AP after middle boundary for 2-exponential fit
thresh = -20
$o3 = new Vector()
firing = 0
check($o1,$o2,$o3,0,thresh,firing)
for i = 0, $o1.size-1 {
firing = check($o1,$o2,$o3,i,thresh,firing)
}
/***
if( $5 ) {
print $s4, "Found ",$o3.size," APs:"
for i = 0, $o3.size-1 {
printf("%g.\t(x,y)[%d] = (%g,%g)\n",i,$o3.x[i],$o1.x[$o3.x[i]],$o2.x[$o3.x[i]])
}
}
***/
// find first AP after, and last AP before, the time bounds
$o6.val = $o1.ind($o3).indwhere(">=", boundary.x[0])
// probably, keep left boundary less than zero.
// if( $o6.val < 0 ) $o6.val = $o3.size
$o7.val = $o1.ind($o3).indwhere(">", boundary.x[1]) - 1
if( $o7.val < 0 ) $o7.val = $o3.size-1
// find first AP after the cutoff for the first exponential
$o10.val = $o1.ind($o3).indwhere(">=", bound_2exp)
$o8.val = $o7.val-$o6.val+1
if( $o7.val < 0 || $o6.val < 0 ) {
$o8.val = 0
}
// if the AP window begins AFTER the AP peak, it is possible
// that the last AP will have no window associated with it.
if( ($o7.val >= 0) && ($o3.size > 0)) {
if ($o1.x[$o3.x[$o7.val]]+apwin.x[0] < boundary.x[1]) {
$o9.val = $o8.val
} else { $o9.val = $o8.val - 1 }
} else $o9.val = 0
if( VERBOSE ) {
//printf("left = %g, right = %g, count = %g\n",$o6.val,$o7.val,$o8.val)
if( $o6.val > 0 ) {
printf("first %s after bound %g has time %g\n",$s4,boundary.x[0],$o1.x[$o3.x[$o6.val]])
}
if( $o10.val > 0 ) {
printf("first %s after 2-exp cutoff %g has time %g\n",$s4,bound_2exp,$o1.x[$o3.x[$o10.val]])
}
if( $o7.val > 0 ) {
printf("last %s before bound %g has time %g\n",$s4,boundary.x[1],$o1.x[$o3.x[$o7.val]])
}
printf("Total valid APs %d\n",$o8.val)
printf("Total valid AP windows %d\n",$o9.val)
}
}
func check() { local idx, thresh, firing, time
// $o1 x component, $o2 y component of a voltage trace
// $o3 stores AP times
idx = $4
thresh = $5
firing = $6
if( $o2.x[idx] >= thresh && !firing ) {
firing = 1
time = $o1.x[idx]
$o3.append(idx)
}
if(firing && $o2.x[idx] < thresh && $o1.x[idx] > time) { firing = 0 }
return firing
}
/*******************************************************************************
Find the mean & variance of FR. Report (mean) CV. For multiple model
runs, also report the variance of the CV.
The user has specified an interval, frcv_int, over which the mean and
std values of FR and CV will be calculated. For each non-overlapping
interval of this length contained within the matching bounds, values
of FR and CV will be calculated.
*******************************************************************************/
proc findFR() { local i, nint, cvmean
// $o1 AP times
// $o2 instantaneous FRs
// $o3 mean FR vector
// $o4 stdv FR vector
// $o5 mean CV vector
// $o6 stdv CV vector
//// CHECK THIS FOR BUGS! JAN 05 05
$o2 = new Vector()
// vector $o1 has all of the APtimes throughout the simulation;
// calculate associated firing rates.
frwhere = $o1.indwhere(">=",boundary.x[1])
if( frwhere == -1 ) frwhere = $o1.size
if( EFUN_DBG && VERBOSE ) dbgfile.printf("%% last index should be %d\n",frwhere)
//for i = 0, $o1.size-2 {
for i = 0, frwhere-2 {
$o2.append(1000/($o1.x[i+1]-$o1.x[i]))
}
// how many intervals will be calculated?
//
// if match region is less than specified interval, use the entire match
// region to calculate FR & CV
//
nint = int( (boundary.x[1]-boundary.x[0]) / frcv_int )
if( nint <= 0 ) { nint = 1 }
if( EFUN_DBG && VERBOSE ) {
dbgfile.printf("\n%% Looking at APs: ")
$o1.printf(dbgfile,"%g\t")
dbgfile.printf("%% FRs are : ")
$o2.printf(dbgfile,"%g\t")
}
// frindx stores the indices of the start/endpoints of each interval.
// not needed outside this procedure.
frindx = new Vector()
for( i = boundary.x[0]; i <= boundary.x[1]; i += frcv_int ) {
frwhere = $o1.indwhere(">=",i)
if( VERBOSE ) {
dbgfile.printf("%% %s > %g at index %d\n",$o1,i,frwhere)
}
if( frwhere == -1 || frwhere > $o2.size-1 ) {
// all firing occurs before this time boundary; the last FR
// value is the appropriate index.
frindx.append($o2.size)
} else {
frindx.append(frwhere)
}
}
for i = frindx.size, nint+1 {
frindx.append($o2.size)
}
if( EFUN_DBG && VERBOSE ) {
dbgfile.printf("%% interval bound indices are : ")
frindx.printf(dbgfile,"%g\t")
}
if( $o2.size == 0 ) {
$o3.val = $o4.val = $o5.val = $o6.val = 0
return
}
frtmp = new Vector(nint)
frtmpsd = new Vector(nint)
cvtmp = new Vector(nint)
for i = 0, nint-1 {
if( frindx.x[i] == frindx.x[i+1] ) {
frtmp.x[i] = 0
} else {
if( (frindx.x[i+1] - frindx.x[i]) < 2 ) {
// BUG! found 15 mar
// frtmp.x[i] = $o2.x[i]
frtmp.x[i] = $o2.x[frindx.x[i]]
} else {
frtmp.x[i] = $o2.mean(frindx.x[i],frindx.x[i+1]-1)
}
}
if( (frindx.x[i+1] - frindx.x[i]) < 2 ) {
frtmpsd.x[i] = 0
cvtmp.x[i] = 0
} else {
frtmpsd.x[i] = $o2.stdev(frindx.x[i],frindx.x[i+1]-1)
cvtmp.x[i] = frtmpsd.x[i] / frtmp.x[i]
}
if( EFUN_DBG && VERBOSE ) {
dbgfile.printf("%% Interval %d firing rates: ",i)
if( frindx.x[i]==frindx.x[i+1] ) {
dbgfile.printf("None here.\n")
} else { $o2.printf(dbgfile,"%g\t",frindx.x[i],frindx.x[i+1]-1) }
dbgfile.printf("%% Mean FR = %g, FR SD = %g, CV = %g\n",frtmp.x[i],frtmpsd.x[i],cvtmp.x[i])
}
}
$o3.val = frtmp.mean
$o4.val = frtmpsd.mean
$o5.val = cvtmp.mean
if( cvtmp.size < 2 ) {
$o6.val = 0
} else {
$o6.val = cvtmp.stdev
}
if( EFUN_DBG && VERBOSE ) {
dbgfile.printf("%% CV info: ")
cvtmp.printf(dbgfile,"%g\t")
dbgfile.printf("%% CV mean = %g, CV std = %g\n\n",$o5.val,$o6.val)
}
}
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
}
}
if( bound_2exp == 0 ) bound_2exp = boundary.mean
}
set_allxy()
}
proc ask_apwin() { local ws, we, w_scl
//if (string_dialog("For AP shape error, enter (1) START and
//(2) END of time window, relative to AP peak, and
//(3) TOTAL WEIGHT for shape error",tmpstr)){
apwin.x[0] = $1
apwin.x[1] = $2
shape_scale = $3
}
proc set_apwin() {
ask_apwin()
set_relative_idx()
set_allxy()
if( mintstop_ < boundary.x[1] ) { mintstop_ = boundary.x[1] }
//print "Set AP window to [",apwin.x[0],",",apwin.x[1],"] -- check [",xtmp.min,",",xtmp.max,"]"
//print "mintstop = ",mintstop_
}
proc set_apwin_noGUI() {
apwin.x[0] = $1
apwin.x[1] = $2
set_relative_idx()
set_allxy()
if( mintstop_ < boundary.x[1] ) { mintstop_ = boundary.x[1] }
//print "Set AP window to [",apwin.x[0],",",apwin.x[1],"] -- check [",xtmp.min,",",xtmp.max,"]"
}
proc set_relative_idx() { local tmin, tmax, i,j,k, ws, we
//print "Entering set_relative_idx() "
left_ptr = new Pointer(&exp_APleft)
mid_ptr = new Pointer(&exp_APmid)
right_ptr = new Pointer(&exp_APright)
n_ptr = new Pointer(&n_expAP)
nw_ptr = new Pointer(&n_expWin)
find_APtimes(xdat,ydat,xyAPtimes,"Exp AP",1,left_ptr,right_ptr,n_ptr,nw_ptr,mid_ptr)
//objref tmp_idx
if( PASTE_DBG ) {
print "in set_relative_idx()"
dbgfile=new File()
dbgfile.wopen("set_win.m")
dbgfile.printf("apwin = [")
apwin.printf(dbgfile)
dbgfile.printf("];\n")
}
idx_rel_list = new List()
xdat_rel_list = new List()
ydat_rel_list = new List()
for k = exp_APleft, exp_APright {
//for iloop = 0, n_expWin-1 {
//k = exp_APleft + iloop
ws = apwin.x[0]
we = apwin.x[1]
if( apwin.x[0] + xdat.x[ydat.nextpeak(xyAPtimes.x[k])] < boundary.x[0] ) {
ws = boundary.x[0] - xdat.x[ydat.nextpeak(xyAPtimes.x[k])]
printf("Reset ws = %g - %g = %g.\t",boundary.x[0],xdat.x[ydat.nextpeak(xyAPtimes.x[k])],ws)
}
//if( apwin.x[0] + xdat.x[ydat.nextpeak(xyAPtimes.x[k])] > boundary.x[1] ) {
// ws = boundary.x[1]
//printf("Reset ws = %g.\t",boundary.x[1])
//}
if( apwin.x[1] + xdat.x[ydat.nextpeak(xyAPtimes.x[k])] > boundary.x[1] ) {
we = boundary.x[1] - xdat.x[ydat.nextpeak(xyAPtimes.x[k])] - dt
//printf("Reset we = %g - %g -%g = %g.\t",boundary.x[1],xdat.x[ydat.nextpeak(xyAPtimes.x[k])],dt,we)
}
//printf("Set APWin %d to [%d, %d]\n",k,ws,we)
tmin = xdat.indwhere(">=",ws+xdat.x[ydat.nextpeak(xyAPtimes.x[k])])
//if( tmin < 0 ) tmin = xdat.size-1
tmax = xdat.indwhere(">=",we+xdat.x[ydat.nextpeak(xyAPtimes.x[k])])
if( tmax < 0 ) tmax = xdat.size-1
if( tmin < 0 ) tmin = tmax
//printf("Setting window %d indices to [%d, %d]\n",k,tmin,tmax)
if( tmin > tmax ) tmin = tmax
//printf("Had tmin > tmax; reset them to [%d, %d]\n",tmin,tmax)
xtmp = xdat.c(tmin,tmax)
ytmp = ydat.c(tmin,tmax)
if( PASTE_DBG ) {
print "apwin[",k,"] = (",ws,", ",we,")"
dbgfile.printf("xdat(:,%d)=[",k+1)
xdat.printf(dbgfile)
dbgfile.printf("];\nydat(:,%d)=[",k+1)
ydat.printf(dbgfile)
dbgfile.printf("];\n")
}
tmp_idx = new Vector(0)
for i=0,xtmp.size-1 {
if( PASTE_DBG ) {
dbgfile.printf("%%adding %d (%g, %g)\n",i,xtmp.x[i],ytmp.x[i])
}
j = pick_x(xtmp.x[i])
tmp_idx.append(j)
if( PASTE_DBG ) { if(i%500 == 0) { print "...done" } }
}
if( PASTE_DBG ) {print "done the loop, now sort()" }
tmp_idx.sort()
idx_rel_list.append(tmp_idx.c)
}
if( PASTE_DBG ) {
print "done set_relative_idx()"
dbgfile.close()
print "Leaving set_relative_idx() "
print "SET_REL_IDX:"
if( ydat_rel_list.count > 0 ) {
for i = 0, ydat_rel_list.count-1 {
print "\ty[",i,"] peak ",ydat_rel_list.object(i).x[ydat_rel_list.object(i).firstpeak]
}
} else { print "\tydat not set yet" }
}
}
proc print_Mfile_header() { local e
e = $3
dbgfile.printf("mod_APleft = %g; mod_APright = %g;\n",mod_APleft,mod_APright)
dbgfile.printf("exp_APleft = %g; exp_APright = %g;\n",exp_APleft,exp_APright)
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)
if( ismembrane("itGHK" )) {
dbgfile.printf("pcabar_CaT=%g;\n",pcabar_itGHK)
}
if( ismembrane("cat" )) {
dbgfile.printf("gbar_cat=%g;\n",gbar_cat)
}
}
/****
dbgfile.printf("lbl(%d) = {'Target FR %.1f Hz, CV %g'};\n",nlbl,datFRmean,datCVmean)
nlbl += 1
dbgfile.printf("lbl(%d) = {'Model FR %.1f Hz, CV %g'};\n",nlbl,modFRmean,modCVmean)
nlbl += 1
dbgfile.printf("lbl(%d) = {'F_{shape} = %g, F_{SubShape} = %g'};\n",\
nlbl,shape_scale, subshape_scale)
nlbl += 1
****/
dbgfile.printf("lbl(%d) = {'F_{shape} = %g, F_{fr} = %g'};\n",nlbl,shape_scale,tot_frscale)
nlbl += 1
dbgfile.printf("lbl(%d) = {'F_{slope} = %g'};\n",nlbl,slope_scale)
nlbl += 1
dbgfile.printf("lbl(%d) = {'Burst ISI threshold = %g ms, penalty %g;'};\n", nlbl,brst_thold,brst_pnlty)
nlbl += 1
dbgfile.printf("lbl(%d) = {'1st spike delay penalty = %g;'};\n", nlbl,delay_pnlty)
nlbl += 1
dbgfile.printf("lbl(%d) = {'AP Window [%g, %g]'};\n",nlbl,apwin.x[0],apwin.x[1])
nlbl += 1
dbgfile.printf("apwin = [%g %g];\n",apwin.x[0],apwin.x[1])
dbgfile.printf("lbl(%d) = {' '};\n",nlbl)
nlbl += 1
dbgfile.printf("lbl(%d) = {'Extra silent model FR penalty = %g;'};\n", nlbl,nf_frscale)
//dbgfile.printf("tgt_val = [0.01 0.004 0.001 0.05 0.0125 0 0.00025 5e-5];\n")
//dbgfile.printf("mdl_val = [%g %g %g %g %g %g %g %g];\n",gnabar_fn,gkbar_fn,\
// gbar_kca,Kp_cad,Rca_cad,gbar_ka,gbar_cahi,gbar_nap)
//dbgfile.printf("lblstr = {' g_Na','g_KDR','g_KCa',' K_p',' R_Ca',' g_KA',' g_Ca','g_NaP'};\n")
if( GBAR_SMRY ) {
dbgfile.printf("tgt_val = [0.01 0.004 0.001 0.05 0.0125 0 0.00025 5e-5 3e-5 0];\n")
dbgfile.printf("mdl_val = [%g %g %g %g %g %g %g %g %g ",gnabar_fn,gkbar_fn,\
gbar_kca,Kp_cad,Rca_cad,gbar_ka,gbar_cahi,gbar_nap,g_pas)
if( ismembrane("itGHK") ) {
dbgfile.printf("%g",pcabar_itGHK)
}
if( ismembrane("cat") ) {
dbgfile.printf("%g",gbar_cat)
}
dbgfile.printf(" ];\n")
dbgfile.printf("lblstr = {' g_Na','g_KDR','g_KCa',' K_p',' R_Ca',' g_KA',' g_Ca','g_NaP','g_L','g_CaT'};\n")
}
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
if( ismembrane("itGHK" )) {
dbgfile.printf("lbl(%d) = {'p_{CaT} = %g'};\n",nlbl,pcabar_itGHK)
nlbl += 1
}
if( ismembrane("cat" )) {
dbgfile.printf("lbl(%d) = {'g_{CaT} = %g'};\n",nlbl,gbar_cat)
nlbl += 1
}
}
******/
dbgfile.printf("boundary = [%g %g];\n",boundary.x[0],boundary.x[1])
if( VERBOSE ) {
dbgfile.printf("%% Exp APs during interval = %d\n",n_AP)
dbgfile.printf("xmodel = [")
$o2.printf(dbgfile)
dbgfile.printf("];\n")
dbgfile.printf("ymodel = [")
$o1.printf(dbgfile)
dbgfile.printf("];\n")
dbgfile.printf("xexpt = [")
xdat.printf(dbgfile)
dbgfile.printf("];\n")
dbgfile.printf("yexpt = [")
ydat.printf(dbgfile)
dbgfile.printf("];\n")
dbgfile.printf("expAPtimes = [")
xyAPtimes.printf(dbgfile)
dbgfile.printf("];\n")
dbgfile.printf("modAPtimes = [")
modelAPtimes.printf(dbgfile)
dbgfile.printf("];\n")
}
nstr = 1
}
proc print_Mfile_tail() { local e
e = $3
print "Total Error\t",e
dbgfile.printf("e_final = %g;\n",e)
//dbgfile.printf("str(%d) = {'Total Error %g'};\n",nstr,e)
//nstr += 1
dbgfile.printf("lbl(%d) = {'Total Error %g'};\n",nlbl,e)
nlbl += 1
print "Raw voltage trace"
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("y1=[min(yexpt) max(yexpt)];\n")
dbgfile.printf("\nsubplot(2,1,1);\n")
dbgfile.printf("pl = plot(xexpt,yexpt,'--',xmodel,ymodel,'-',x1,y1,'r',x2,y1,'r');\n")
dbgfile.printf("set(pl(1),'LineWidth',2); set(pl(2),'LineWidth',2); \n")
dbgfile.printf("xlim([boundary(1)-10 boundary(2)+10]);\n")
dbgfile.printf("legend('target','model');\n")
dbgfile.printf("title(ttl);\n")
dbgfile.printf("\nsubplot(2,2,3);\n")
if( GBAR_SMRY ) {
dbgfile.printf("xmx=plot_deviation_from_tgt(tgt_val(end:-1:1),mdl_val(end:-1:1),lblstr(end:-1:1),'');\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")
dbgfile.printf("text(1.2*xmx,length(tgt_val)/2,lbl,'FontSize',10);\n")
if( VERBOSE ) { dbgfile.printf("orient(gcf,'tall');saveas(gcf,'%s_V.fig');\n\n",outname) }
//
// print out time-modulated FR fits
//
print "Now, time-modulated FR fits"
dbgfile.printf("figure(2);\n")
dbgfile.printf("h2 = axes('Position',[0 0 1 1],'Visible','off');\n")
dbgfile.printf("ttl2 = '%s: Total time-modulated FR Error %g';\n",outname,efr)
if( VERBOSE ) {
if( n_AP > 3 ) {
dbgfile.printf("datFit=[ ")
if( datFit.size > 0 ) { datFit.printf(dbgfile) }
dbgfile.printf("];\n")
dbgfile.printf("modFit=[ ")
if( modFit.size > 0 ) { modFit.printf(dbgfile) }
dbgfile.printf("];\n")
if( modFR.size > 0 ) {
dbgfile.printf("modFR=[ ")
modFR.printf(dbgfile)
dbgfile.printf("];\n")
} else { dbgfile.printf("modFR=[];\n") }
if( datFR.size > 0 ) {
dbgfile.printf("datFR=[ ")
datFR.printf(dbgfile)
dbgfile.printf("];\n")
} else { dbgfile.printf("datFR=[];\n") }
dbgfile.printf("axes('Position',[.1 .5 .8 .4]);\n")
dbgfile.printf("f1=[min([datFR ; modFR]) max([datFR ; modFR])];\n")
if( modFR.size < 1 ) {
// recall, Matlab starts their indices at 1, not 0
dbgfile.printf("frpl = plot(xexpt(expAPtimes(1:%d)),datFR,'b.',xexpt(expAPtimes(%d:%d)),datFit,'b-',[%f %f], [0 0],'g+',xexpt(expAPtimes(%d:%d)),modFit,'g--',x1,f1,'r',x2,f1,'r');\n",\
exp_APright,exp_APleft+1,exp_APright,boundary.x[0],boundary.x[1],\
exp_APleft+1, exp_APright)
} else {
dbgfile.printf("frpl = plot(xexpt(expAPtimes(1:%d)),datFR,'b.',xexpt(expAPtimes(%d:%d)),datFit,'b-',xmodel(modAPtimes(1:%d)),modFR,'g+',xexpt(expAPtimes(%d:%d)),modFit,'g--',x1,f1,'r',x2,f1,'r');\n",\
exp_APright, exp_APleft+1, exp_APright, mod_APright, exp_APleft+1, exp_APright)
}
dbgfile.printf("set(frpl(2),'LineWidth',2); set(frpl(4),'LineWidth',2); \n")
dbgfile.printf("legend('target','target fit','model','model fit');\n")
dbgfile.printf("title(ttl2);\n")
}
}
dbgfile.printf("set(gcf,'CurrentAxes',h2);\n")
dbgfile.printf("text(.05,.25,lbl,'FontSize',12);\n")
dbgfile.printf("text(.5,.25,str2,'FontSize',12);\n")
if( VERBOSE ) { dbgfile.printf("orient(gcf,'tall');saveas(gcf,'%s_FR.fig');",outname) }
}
/**************************************************
Determine whether there are 2 components to the
AHP of each action potential. Between 2 APs,
the sign changes in dV/dt will give how
many components there are to the AP.
Sign change 1 = AP peak
Sign change 2 = bottom of spike
if there are 2 more changes,
Sign change 3 = peak between quick & slow AHPs
Sign change 4 = trough of slower AHP
For each spike with 2 AHP components, record
the depths of the quick & slow AHP troughs, and
the ratio between them.
input $o1 time vector
$o2 voltage vector
$o3 AP times
$4 vector with start, end bounds
output $o5 time of each trough
$o6 voltage of each trough
**************************************************/
proc find_APtrough() { local i, k, start_ind,end_ind, cnt
// Identify boundaries of the AHP regions.
dvdt = new Vector()
dvdt.deriv($o2,$o1.x[1]-$o1.x[0],2)
// calculate AHP data for each spike separately
for k=0,$o3.size-1 {
// set start & end points for AHP checks: start when this AP crosses threshold; stop at next AP's threshold, or end of the data trace.
start_ind = $o3.x[k]
if( k == $o3.size-1 ) {
end_ind = $o1.indwhere(">=",$o4.x[1])
} else { end_ind = $o3.x[k+1] - 1 }
//if( VERBOSE ) printf("Finding trough of AP %d, from %g to %g\n",k,$o1.x[start_ind],$o1.x[end_ind])
cnt = 0
for i=start_ind,end_ind-1 {
if( (dvdt.x[i] >= 0 && dvdt.x[i+1] <= 0) || (dvdt.x[i] <= 0 && dvdt.x[i+1] >= 0)){
cnt += 1
//printf("\tSign change %d at %d (%g,%g)",cnt,i,$o1.x[i],$o2.x[i])
if( cnt == 2 ) {
$o5.append(i)
$o6.append($o2.x[i])
//printf("\tTROUGH")
}
//printf("\n")
}
if( cnt == 2 ) break
}
//printf("cnt = %d; on to next AP\n\n",cnt)
}
if( VERBOSE ) {
printf("\tFound trough times: \t")
$o5.printf
printf("\tFound trough voltages:\t")
$o6.printf
}
}
endtemplate TrajDens_fit2