begintemplate ModelViewSeclistAnalysis

public all, analyse, mkdisplay, selected, dmlist, parmsubsets, parmsecs, parminhomos, allseclist
public selsub, selsec, disp, grph
public subsets, prependname, mview, unlink
external hoc_sf_, nil, lambda_f
objref this, all, tobj, tobj1, allseclist, mview
objref mt, ms, subsets, subset, dislist_, distinctvalues
objref parmsubsets, parmsecs, parminhomos, dmech, dmlist
objref pmat, segmat, sortx, psub, ionstyle
objref nsegvec, lvec, lamvec, disp, msubstr, msub, disp1
strdef tstr, mname, pname, pfullname

proc init() {local i
	analysed = 0
	mview = $o1
	distinctvalues = mview.distinctvalues
	allseclist = $o2
	all = new List()
	forsec allseclist all.append(new SectionRef())
}

proc unlink() {local i
	if (dmlist != nil) for i=0, dmlist.count - 1 {
		dmlist.object(i).unlink()
	}
	objref all, tobj, tobj1, allseclist, mview
	objref mt, ms, subsets, subset, dislist_
	objref parmsubsets, parmsecs, parminhomos, dmech, dmlist
	objref pmat, segmat, sortx, psub, ionstyle
	objref nsegvec, lvec, lamvec, disp, msubstr, msub, disp1
}

proc prependname() {
	all.object($1).sec sprint($s2, "%s.%s", secname(), $s2)
}

proc analyse() {
	analysed = 1
	if (all.count == 0) { return }
	parmsubsets = new ModelViewSubsetList()
	parmsubsets.subsets.append(new ModelViewParmSubset(new Vector(all.count), this))
	parmsubsets.subsets.object(0).subset.indgen
	mk_existence_subsets()
	parmsecs = new List()
	parminhomos = new List()
	dislist_ = new List()
	mk_parm_subsets()
}

proc mkdisplay() {local i, j, k
	disp = $o1

	if (!analysed) {
		analyse()
	}

	mk_nseg()
	msub = new TreeViewItem(disp, "", this, -1)
	j = 0
	for i=0, dmlist.count - 1 {
		dmech = dmlist.object(i)
		if (object_id(dmech.location) == 0) { continue }
		dmech.display(msub, tstr)
		j += 1
	}
	sprint(tstr, "%d inserted mechanisms", j)
	msub.str(tstr)

	if (parmsubsets.subsets.count > 1) {
		sprint(tstr, "%d subsets with constant parameters", parmsubsets.subsets.count)
		disp1 = new TreeViewItem(disp, tstr, this, -1)
	}else{
		disp1 = disp
	}
	for i=0, parmsubsets.subsets.count-1 {
		tobj = parmsubsets.subsets.object(i)
		sprint(tstr, "%s (%d sections)", tobj, tobj.subset.size)
		msub = new TreeViewItem(disp1, tstr, tobj, -1)
		prsub(tobj.subset)
		for j=0, tobj.parmlist.count-1 {
			tobj1 = tobj.parmlist.object(j)
			sprint(tstr, "%s = %g", tobj1.parmelm.name, tobj1.value)
			msub.append(tstr, tobj, j)
		}
	}

	if (parmsecs.count + parminhomos.count > 1) {
		sprint(tstr, "%d sections with unique parameters", parmsecs.count + parminhomos.count)
		disp1 = new TreeViewItem(disp, tstr)
	}else{
		disp1 = disp
	}
	for i=0, parmsecs.count-1 {
		tobj = parmsecs.object(i)
		all.object(tobj.isec).sec sprint(tstr, "%s {", secname())
		msub = new TreeViewItem(disp1, tstr, tobj, -1)
		for j=0, tobj.parmlist.count-1 {
			tobj1 = tobj.parmlist.object(j)
			sprint(tstr, "    %s = %g", tobj1.parmelm.name, tobj1.value)
			msub.append(tstr, tobj, j)
		}
		msub.append("}", tobj, j)
	}

	for i=0, parminhomos.count-1 {
		tobj = parminhomos.object(i)
		all.object(tobj.isec).sec sprint(tstr, "%s {", secname())
		msub = new TreeViewItem(disp1, tstr, tobj, -1)
		for j=0, tobj.parmlist.count-1 {
			tobj1 = tobj.parmlist.object(j)
			sprint(tstr, "    %s =", tobj1.parmelm.name)
			for k=0, tobj1.vec.size-1 {
				sprint(tstr, " %s %g", tstr, tobj1.vec.x[k])
			}
			msub.append(tstr, tobj, j)
		}
		msub.append("}", tobj, j)
	}
}

proc selected() {
	if ($4 == -2) {
		mview.gui.selsub()
	}else if ($4 == -1) {
		mview.gui.selsub(allseclist)
	}else if ($4 > -1) {
		seclist(tobj, subsets.object($4))
		mview.gui.selsub(tobj)
	}
}

proc selsub() {
	seclist(tobj, $o1)
	mview.gui.selsub(tobj)
}

proc selsec() {
	mview.gui.selsec(all.object($1))
}

proc seclist() {local i
	$o1 = new SectionList()
	for i=0, $o2.size-1 {
		all.object($o2.x[i]).sec $o1.append()
	}
}

// subsets of density mechanism existence
// subsets is a list of distinct subsets
// the objects of subsets are vectors whose elements are the indices of
//  the all list of SectionRef defining the cell.
// typevec,subsetvec are parallel vectors. typevec elements are
//   the mechanism types of actually existing density mechanisms and
//   subsetvec elements are the subsets list index of their locations
//   i.e. multiple mechanisms can share same subset.
proc mk_existence_subsets() { local b, it, i
	dmlist = new List()
	mt = new MechanismType(0)
	dmlist.append(new ModelViewDMech("Ra")) // so 0 is Ra and rest are type
	dmlist.object(0).location = parmsubsets.subsets.object(0)
	for it = 1, mt.count-1 { // skip morphology
		mt.select(it)
		mt.selected(mname)
		dmlist.append(new ModelViewDMech(mname))
		tobj = new Vector()
		for i=0, all.count-1 all.object(i).sec if (ismembrane(mname)) {
			tobj.append(i)
		}
		if (tobj.size) {
			i = parmsubsets.isubset(tobj, this)
			tobj = parmsubsets.subsets.object(i)
			dmlist.object(dmlist.count-1).location = tobj
		}
	}
}

// If a parameter in an existence subset is not uniform then there
// is a problem about how to succinctly express the variation.
// Uniform is easy, just rely on the existence subset.
// If there are only a few distinct values then perhaps they should define
//   subsubsets.

proc prsub() {local i
return
	for i=0, $o1.size-1 {
		all.object($o1.x[i]).sec printf(" %s", secname())
	}
	printf("\n")
}

proc mk_parm_subsets() {local nrow, ncol, irow, icol, itype, iname, isec, is, x, size, b, doit

	// deal with Ra
	tobj = new Vector(all.count)
	for i=0, all.count-1 {
		tobj.x[i] = all.object(i).sec.Ra
	}
	dmech = dmlist.object(0)
	subset = dmech.location.subset
	// just remember the following are 0 and 1
	dmech.parms.append(new ModelViewDParm(nil, 0, "Ra"))
	dmech.parms.append(new ModelViewDParm(nil, 0, "rallbranch"))
	parm_assess(tobj.c, 0, 0)

	// only deal with rallbranch != 1
	tobj = new Vector(all.count)
	for i=0, all.count-1 {
		x = all.object(i).sec.rallbranch
		if (x != 1) {
			tobj.x[i] = x
		}else{
			tobj.x[i] = -1e50
		}
	}
	parm_assess(tobj.c, 1, 0)

	// deal with all density mechanisms
	for itype=1, dmlist.count-1 {
		dmech = dmlist.object(itype)
		if (object_id(dmech.location) == 0) { continue }
		mt.select(itype)
		mt.selected(mname)
		if (hoc_sf_.tail(mname, "_ion", tstr) != -1) {
			// ions need special handling. i.e. use ena or nai,nao
			is_ion = 1
			ms = new MechanismStandard(mname, 0)
		}else{
			// parameters only
			is_ion = 0
			ms = new MechanismStandard(mname, 1)
		}
		dmech.fill(ms)
		subset = dmech.location.subset
		// make a matrix of all params(columns) for all sections(rows)
		// also make a list of inhomogeneous sections
		// how many params
		ncol = 0
		for iname=0, ms.count-1 {
			ncol += ms.name(pname, iname)
		}
		if (ncol == 0) { continue }
		// how many sections
		nrow = subset.size
		if (is_ion) {
			ionstyle = new Vector(subset.size)
		}
		pmat = new Matrix(nrow, ncol)
		// fill by rows
		irow = 0
		for isec=0, subset.size-1 all.object(subset.x[isec]).sec {
			segmat = new Matrix(nseg, ncol)
			if (is_ion) {
				x = ion_style(mname)
				// style=1 means concentrations are parameters
				// 2 means reversal potential is a parameter
				// 0 means neither are parameters
				x = ((x%4) == 1) + 2*((int(x/8)%4) == 1)
				ionstyle.x[isec] = x
			}
			iseg = 0
			for (x, 0) {
				ms.in(x)
				icol = 0
				for iname=0, ms.count-1 {
					size = ms.name(pname, iname)
					for is=0, size-1 {
						segmat.x[iseg][icol] = ms.get(pname, is)
						icol += 1
					}
				}
				iseg += 1
			}
			// now look at within section uniformity
			icol = 0
			b = 1
			for iname=0, ms.count-1 {
				size = ms.name(pname, iname)
				for is=0, size-1 {
					if(distinctvalues.ndistinct(segmat.getcol(icol)) > 1) {
						doit = 1
						if (is_ion) {
							doit = 0
if (ionstyle.x[isec] == 2 && icol == 0) { doit = 1} //ena
if (ionstyle.x[isec] == 1 && (icol == 1 || icol == 2)) { doit = 1 } // nai and nao
						}
						if (doit == 0) { continue }
						addparm2inhomo(isec, segmat.getcol(icol), iname, is)
						// sentinal to avoid adding below
						segmat.x[0][icol] = -1e50
					}
					icol+=1
				}
			}
			pmat.setrow(isec, segmat.getrow(0))
		}
		// now look at the columns (parameter vs position)
		icol = 0
		if (is_ion) {
			handle_ion(ms, pmat, ionstyle)
		}else{
			for iname=0, ms.count-1 {
				size = ms.name(pname, iname)
				for is=0, size-1 {
					parm_assess(pmat.getcol(icol), iname, is)
					icol += 1
				}
			}
		}
	}
}

proc handle_ion() {local i
	if (distinctvalues.ndistinct($o3) > 1) {
		printf("ions with multiple styles not handled")
	}else{
		if ($o3.x[0] == 2) { //ena
			parm_assess($o2.getcol(0), 0, 0)
		}
		if ($o3.x[0] == 1) { // nai and nao
			for i=1,2 {
				parm_assess($o2.getcol(i), i, 0)
			}
		}
	}
}

func parm_assess() {local n
	n = distinctvalues.ndistinct($o1)
	if (n == 0) { return n }
//	if (n > 20) {
//		sprint(tstr, "%d distinct constant values ranging from %g to %g", \
//			n, distinctvalues.min, distinctvalues.max)
//		dmech.parms.object($2).parmelms.object($3).append( \
//			new ModelViewParmStr(tstr, dmech.location))
//		return n
//	}
	distinctvalues.dislist(dislist_)
	if (n > 1 || $o1.size == 1) {
		for i=0, dislist_.count-1 {
			tobj1 = dislist_.object(i)
			val = $o1.x[tobj1.x[0]]
			if (val == -1e50) { continue }
			if (tobj1.size > 1) {
				addparm2subset(tobj1, val, $2, $3)
			}else{
				addparm2sec(tobj1.x[0], val, $2, $3)
			}
		}
	}else if (n == 1 && $o1.x[0] != -1e50) {
		addparm2subset(dislist_.object(0), $o1.x[0], $2, $3)
	}
	if (n > 1) {
		dmech.grph = this
		dmech.parms.object($2).parmelms.object($3).dmech = dmech
	}
	return n
}

// args are subset indices, val, name index, array index
proc addparm2subset() {local i, val
	tobj = $o1.c.index(subset, $o1).sort
	i = parmsubsets.isubset(tobj, this)
	tobj = parmsubsets.subsets.object(i)
	tobj1=new ModelViewParmVal($2, tobj)
	dmech.append($3, $4, tobj1)
	tobj.parmlist.append(tobj1)
}

proc addparm2sec() {local i, isec, val
	isec = subset.x[$1]
	for i=0, parmsecs.count-1 {
		tobj = parmsecs.object(i)
		if (isec == tobj.isec) {
			tobj1 = new ModelViewParmVal($2,tobj)
			dmech.append($3, $4, tobj1)
			tobj.parmlist.append(tobj1)
			return
		}
	}
	tobj = new ModelViewParmSec(isec, this)
	tobj1=new ModelViewParmVal($2, tobj)
	dmech.append($3, $4, tobj1)
	tobj.parmlist.append(tobj1)
	parmsecs.append(tobj)
}

proc addparm2inhomo() {local i, isec
	isec = subset.x[$1]
	for i=0, parminhomos.count-1 {
		tobj = parminhomos.object(i)
		if (isec == tobj.isec) {
			tobj1 = new ModelViewParmVec($o2,tobj)
			dmech.append($3, $4, tobj1)
			dmech.grph = this
			dmech.parms.object($3).parmelms.object($4).dmech = dmech
			tobj.parmlist.append(tobj1)
			return
		}
	}
	tobj = new ModelViewParmInhomo(isec, this)
	tobj1=new ModelViewParmVec($o2, tobj)
	dmech.append($3, $4, tobj1)
	dmech.grph = this
	dmech.parms.object($3).parmelms.object($4).dmech = dmech
	tobj.parmlist.append(tobj1)
	parminhomos.append(tobj)
}

proc mk_nseg() {local i, n_nseg, n_dx, dx, dlam
	nseg_strategy = 0
	nsegvec = new Vector(all.count)
	lvec = new Vector(all.count)
	lamvec = new Vector(all.count)
	for i=0, all.count-1 all.object(i).sec {
		nsegvec.x[i] = nseg
		lvec.x[i] = L
		lamvec.x[i] = L/lambda_f(100)
	}
	sprint(tstr, "%d sections; %d segments", all.count, nsegvec.sum)
	disp.append(tstr)
	n_nseg = distinctvalues.ndistinct(nsegvec)
	distinctvalues.dislist(dislist_)
	sprint(tstr, "%d distinct values of nseg", n_nseg)
	msub = new TreeViewItem(disp, tstr)
//	for i=0, dislist_.count-1 { printf(" %d of %g\n", dislist_.object(i).size, nsegvec.x[dislist_.object(i).x[0]])}

	i = lvec.c.div(nsegvec).max_ind
	all.object(i).sec {
		dx = L/nseg
		sprint(tstr, "Longest dx is %g at %s with nseg=%d", dx, secname(), nseg)
		msub.append(tstr)
	}
	dx = int(dx/10+.999)*10
	tobj = lvec.c.div(dx).floor.add(1).sub(nsegvec)
	if (tobj.dot(tobj) == 0) {
		sprint(tstr, "all nseg consistent under assumption of dx=%g and nseg=L/dx+1", dx)
		msub.append(tstr)
		nseg_strategy=1
	}else{
		tobj = lvec.c.div(2*dx).add(.45).floor.mul(2).add(1).sub(nsegvec)
		if (tobj.dot(tobj) == 0) {
			sprint(tstr, "nseg consistent under assumption of dx=%g and nseg=2*int(L/dx/2)+1", dx)
			msub.append(tstr)
			nseg_strategy=2
		}else{
			sprint(tstr, "nseg not consistent under constant dx=%g assumption", dx)
			msub.append(tstr)
		}
	}

	
	i = lamvec.c.div(nsegvec).max_ind
	all.object(i).sec {
		dlam = lamvec.x[i]/nseg
		sprint(tstr, "Greatest dlambda is %g at %s with nseg=%d", dlam, secname(), nseg)
		msub.append(tstr)
	}
	if (dlam > .01) {
		if (dlam > .1) {rscl = .1} else { rscl = .01 }
		dlam = int(dlam/rscl+.45)*rscl
		tobj = lamvec.c.div(2*dlam).add(.45).floor.mul(2).add(1).sub(nsegvec)
		if (tobj.dot(tobj) == 0) {
			sprint(tstr, "nseg consistent under assumption of dlambda=%g", dlam)
			msub.append(tstr)
			nseg_strategy=3
		}else{
			sprint(tstr, "nseg not consistent under constant dlambda=%g assumption", dlam)
			msub.append(tstr)
		}
	}
	if (all.count == 1) {
		if (nsegvec.x[0] > 1) {
//			subset = parmsubsets.subsets.object(0).subset
//			addparm2sec(0,"nseg",nsegvec.x[0], -1)
			sprint(tstr, "nseg = %d", nsegvec.x[0])
			msub.append(tstr)
			nseg_strategy = 0
		}
		return
	}
//	if (nseg_strategy == 0) {
//		subset = subsets.object(0)
//		parm_assess(nsegvec, "nseg")
//	}			
}
		
proc grph() {
	mview.gui.grph($o1, $o2)
}

endtemplate ModelViewSeclistAnalysis
