if (name_declared("use_traubexact") == 0) {
use_traubexact = 0
}
begintemplate TraubExactInfo
public tci, traub_parent
objref tci[14], traub_parent[14]
endtemplate TraubExactInfo
objref traubExactInfo
if (use_traubexact) { // for network context traub exact connections
traubExactInfo = new TraubExactInfo()
}
objref tci, traub_parent
strdef tstr
proc traub_connect() {local x, ip, ic localobj child
if (numarg() == 1) {
tci = new Matrix($1, $1, 2)
traub_parent = new Vector($1)
return
}
if (use_traubexact) {
if (traubExactInfo.tci[$o1.type] == nil) {
traubExactInfo.tci[$o1.type] = tci
traubExactInfo.traub_parent[$o1.type] = traub_parent
}
}
if (numarg() == 5) {x = $5} else { x = 1 }
ip = $2 ic = $3
// which is the child and which is the parent?
// assume a child is connected to the parent before it
// makes a delta connection (to child connected to same parent)
// if $3 is already connected assume the child is $2
$o1.comp[$3] child = new SectionRef()
if (child.has_parent) {
$o1.comp[$2] child = new SectionRef()
ip = $3 ic = $2
}
// standard NEURON connection only if child not already connected
if (!child.has_parent) {
$o1.comp[ip] connect child.sec(0), x
traub_parent.x[ic] = ip
}
// save info for optional traub style connections
tci.x[$2][$3] = $4
tci.x[$3][$2] = $4
// note that the above matrix has the topology of a tree unless
// child to child connections form delta loops.
}
// exact() transforms the delta and pyramidal axial resistance branch points to Y's and X's.
// single_cell points at current cell object, .e.g TCR[0]
proc exact(){
print "exact() has been called."
execute("traubexact(single_cell, tci)")
}
proc traubexact() {
//$o1 is the cell, $o2 is the connection matrix
// assume the names of sections range from
// $o1.comp[1] to $o1.comp[$o2.nrow-1]
// the connection matrix
// the traub connection rule is that a childs 0 end is connected
// to the .5 location of the parent and the resistance is defined
// by the connection matrix (uS) with one exception. Sometimes
// pairs or triples of children are also connected to each other
// and these children are moved to the 1 end of the parent and
// the delta network is replaced by the equivalent wye network.
traubexact_topology($o1, $o2, traub_parent)
doNotify()
traubexact_coef($o1, $o2, traub_parent)
}
proc traubexact_topology() { local j, jc, ip, ic, gc, x
// move all children to .5 location of parent except if
// a group of children have a delta topology in which case move
// them to the 1 location of the parent.
for ic = 2, $o2.nrow-1 $o1.comp[ic] {
// parent must be the least index
ip = $o3.x[ic]
x = .5
for j = 0, $o2.sprowlen(ic)-1 {
gc = $o2.spgetrowval(ic, j, &jc)
// the question is whether jc is also connected
// to ic's parent.
if ($o2.getval(ip, jc)) { // yes
x = 1
break
}
}
if (x != parent_connection()) {
disconnect()
$o1.comp[ip] connect $o1.comp[ic](0), x
}
}
}
proc traubexact_coef() { local i, j, id, ip, ic, gp, gc, jc, x localobj dtopol, dmat, dindx, wye
// now that the topology is stable we can safely change the
// connection coefficients so they reflect the desired axial
// resistance
dtopol = new Vector($o2.nrow)
for ic = 2, $o2.nrow-1 $o1.comp[ic] {
ip = $o3.x[ic]
gp = $o2.getval(ip,ic)
if (parent_connection() == 0.5) {
rold = ri(0.5)
scale_connection_coef(0.5, rold*gp)
}else{ // this is delta-topology child.
// so save info for further analysis
dtopol.x[ic] = ip
}
}
// examine the dtopol info. For each delta topology, construct a
// conductance matrix and a section index vector (the parent is always
// the first index) and then find the equivalent wye conductances and
// use those to specify the connection coefficients
for (id = dtopol.indwhere(">",0); id != -1; id = dtopol.indwhere(">",0)){
// a delta exists and the parent is
ip = dtopol.x[id]
// note: error if ip != $o3.x[id]
// so all the children are
dindx = dtopol.c.indvwhere("==", ip)
// prepare for next iteration over id
for i=0, dindx.size-1 { dtopol.x[dindx.x[i]] = 0 }
// and the complete delta index vector is
dindx.insrt(0, ip)
// the delta conductance matrix
dmat = new Matrix(dindx.size, dindx.size)
if (dindx.size < 3) {
execerror("not a delta", 0)
}
for i = 0, dindx.size-1 {
ic = dindx.x[i]
for j = i, dindx.size-1 {
jc = dindx.x[j]
dmat.x[i][j] = $o2.getval(ic,jc)
dmat.x[j][i] = dmat.x[i][j]
}
}
// transform into an equivalent wye
wye = delta2wye(dmat)
// modify the connection coefficients
for i = 0, dindx.size-1 $o1.comp[dindx.x[i]] {
if (i == 0) { x = 1 } else { x = .5 }
rold = ri(x)
scale_connection_coef(x, rold/wye.x[i])
}
}
}
obfunc delta2wye() {localobj g, x, row0
// $o1 is the conductance matrix. return vector is the equivalent
// wye vector of resistances.
// construct a full Gij*(Vi - Vj) = Ii matrix equation
x = new Vector($o1.nrow)
x.fill(1)
$o1.mulv(x, x)
g = new Matrix($o1.nrow, $o1.ncol)
g.setdiag(0, x)
g.add($o1.muls(-1))
//presently indeterminate. save first row and replace by ground equation
row0 = g.getrow(0)
g.setrow(0, 0)
g.x[0][0] = 1
g.inverse($o1)
$o1.getdiag(0, x) // 1 to n-1 are the r0 + ri resistances and x.x[0]=1
// need one more equation (not involving rp) and assume
// numerical accuracy allows an arbitrary pair of nodes for
// current injection and ground.
// so put the row back, select node 1 as the ground point,
// and inject into node 2.
g.setrow(0, row0)
g.setrow(1, 0)
g.x[1][1] = 1
//wasteful but...
g.inverse($o1)
x.x[0] = $o1.x[2][2]
// now the 0 element is r1+r2
// set up the matrix equation to solve for ri
g.zero
g.setcol(0, 1)
g.setdiag(0, 1)
// but the first equation is r1 + r2 = x0
g.x[0][0] = 0
g.x[0][1] = 1
g.x[0][2] = 1
x = g.solv(x)
return x
}
