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
}