collapse a dendritic tree

Managing anatomically complex model cells with the CellBuilder. Importing morphometric data with NEURON's Import3D tool or Robert Cannon's CVAPP. Where to find detailed morphometric data.
mikey
Posts: 37
Joined: Mon Jun 19, 2006 3:06 pm

collapse a dendritic tree

Post by mikey »

Hoc code to collapse a dendritic tree
Collapses a dendritic tree into 3 compartments "proximal", "middle" and "distal".

http://cns.iaf.cnrs-gif.fr/files/COLLAPSE/

It collapses the dendritic tree into 3 compartments. I wonder if anyone really good with .hoc knows how to manipulate the code to instead collapse the tree into 20 compartments (or any other number of compartments specified by the user)? This is a great tool - but I think this manipulation could make it even better. I would be very grateful to anyone that thinks they can do it and takes the time to look.

--------------
bit of BACKGROUND to the algorithm used if you interested:
it does this: "collapsing the dendritic structure into a few compartments based on the conservation of axial resistance, a theme previously proposed by Bush and Sejnowski (1993).Bush–Sejnowski algorithm: Collapsing technique based on conserving Ra rather than the membrane surface area. This is done by making the cross-sectional area of the equivalent cylinder equal to the sum of the cross-sectional areas of all the dendrites represented by that equivalent cylinder. The length of the equivalent cylinder is just the average length of all the dendrites represented by the equivalent cylinder.
-----------------
mikey
Posts: 37
Joined: Mon Jun 19, 2006 3:06 pm

Post by mikey »

Still really struggling with this. So would really really appreciate some help. Below I have cut and paste 1) the original code and 2) my best stab at modifying it (which is not working code - there must be something wrong with it, although I cannot see what) - to collapse a dendritic tree into 4 compartments instead of 3. Although I actually want to collapse it into 40 compartments. Help!! I would be so, so grateful to anyone that can help as this is holding me back quite a lot.

1) THE ORIGINAL
(note that the orignial can also be found here:
http://cns.iaf.cnrs-gif.fr/files/COLLAPSE/
)

Code: Select all



/* ---------------------------------------------------------------------

	HOC CODE TO COLLAPSE A DENDRITIC TREE
	=====================================

   Collapse a dendritic tree into 3 compartments "proximal", "middle"
   and "distal".  

   The collapse is made such as the collapsed dendritic structure
   preserves the axial resistance of the original structure. 
   The algorithm works by merging successive pairs of dendritic 
   branches into an equivalent branch (a branch that preserves the 
   axial resistance of the two original branches).
   This merging of branches can be done according to different methods
   selectable in the present code:

   AREABYLONG

	The code rejects the smallest original branch if it is more than
	5 times smaller than the long orginal branch.  If AREABYLONG is
	set to 1, then the small branch is rescaled to the longest one 
	(a new diameter is assigned such that the small branch conserves
	its axial resistance).  This minimizes errors for dendritic 
	structures where branches have a large range of possible lengths.
	AREABYLONG = 0 is the default.


   AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF

	During the merging of two branches, the length of the equivalent 
	branch is calculated by one out of three possible methods.  

	First, a simple average of the two branch length:

		L = (L1+L2)/2

	This is the method originally used by Bush and Sejnowski
	(J. Neurosci. Methods 46: 159-166, 1993).
	It is selected by setting AVGLENGTH=1.

	Second, by weighting this average with the diameters of each branch

		L = (diam1*L1 + diam2*L2)/(diam1+diam2)

	This algorithm produces more fair merging in the case two branches of
	very different diameters are to be merged.  It is selected by setting
	AVGLENGTHWDIAM=1.

	Third, by weighting the average with membrane area:

		L = (area1*L1 + area2*L2)/(area1+area2)

	This algorithm is better in the case two branches of very different
	areas are to be merged.  It is selected by setting AVGLENGTHWSURF=1.


   CUTOFF1, CUTOFF2

	The dendritic branches are assigned to one of the three compartments
	("proximal", "middle" or "distal") based on their path axial 
	resistance (in Meg-Ohms) to the soma.  The cutoff values (CUTOFF1
	and CUTOFF2) determine these cut-offs.  For example, CUTOFF1=10 means
	that all branches laying within 10 Meg-Ohm of axial resistance from 
	the soma will be merged together into the same ("proximal") section.
	

   Written by M. Neubig & A. Destexhe
   Laval University, 1997; CNRS 2000

-----------------------------------------------------------------------*/

CUTOFF1 = 7.3	// define the cutoff between "proximal" and "middle"
CUTOFF2 = 71	// define the cutoff between "middle" and "distal"

AVGLENGTH = 0		// if 1, merge based on ageraged length
AVGLENGTHWDIAM = 1	// if 1, merge based on avg length, weighted by diam
AVGLENGTHWSURF = 0	// if 1, merge based on avg length, weighted by area

AREABYLONG = 0		// if 1, rescale the length before merging



xopen("tc-cell.geo")		// read geometry file



Lbysurf_areabyasis = 0
Lbysurf_areabylong = 0
Lbydiam_areabyasis = 0
Lbydiam_areabylong = 0
Lbyavg_areabyasis = 0
Lbyavg_areabylong = 0

Lbylong_areabylong = 0

if(AREABYLONG==1) {
  print "Merge by rescaling length"
  if(AVGLENGTH == 1) {
	Lbyavg__areabylong = 1
	print "Merging method based on averaged length"
  } else if(AVGLENGTHWDIAM == 1) {
	Lbydiam_areabylong = 1
	print "Merging method based on averaged length, weigthed by diameters"
  } else if(AVGLENGTHWSURF == 1) {
	Lbysurf_areabylong = 1
  }
} else {
  print "Do not rescale length"
  if(AVGLENGTH == 1) {
	Lbyavg__areabyasis = 1
	print "Merging method based on averaged length"
  } else if(AVGLENGTHWDIAM == 1) {
	Lbydiam_areabyasis = 1
	print "Merging method based on averaged length, weigthed by diameters"
  } else if(AVGLENGTHWSURF == 1) {
	Lbysurf_areabyasis = 1
  }
}

forall nseg=10
forall Ra=260


objectvar secaddress		
objectvar paraddress		
objectvar order			
objectvar parsecndx		
				
objectvar secri			
objectvar secpathri		
objectvar secL			
objectvar secRa			
objectvar secCm			
objectvar mrgL		
objectvar mrgdiam	
objectvar mrgri			
				
objectvar tosec			
objectvar tocyc			
				
objectvar rvp			
objectvar slsoma		
objectvar slroot		

objectvar cycL
objectvar cycdiam
objectvar cycnin





proc initsec(){
	secaddress=new Vector(NSEC,0)		//set/reset
	paraddress=new Vector(NSEC,0)		//set/reset
	order=new Vector(NSEC,0)		//set/reset
	parsecndx=new Vector(NSEC,0)		//set/reset

	secri=new Vector(NSEC,0)		//set/reset
	secpathri=new Vector(NSEC,0.0)		//set/reset
	secL=new Vector(NSEC,0)   		//set/reset
	secRa=new Vector(NSEC,0)   		//set/reset
	secCm=new Vector(NSEC,0)		//set/reset
}
proc initmrg(){
	mrgL=new Vector(NSEC,0)   		//set/reset
	mrgdiam=new Vector(NSEC,0)		//set/reset
	mrgri=new Vector(NSEC,0)		//set/reset

}
proc initcyc(){
	cycL=new Vector(4,0)
	cycdiam=new Vector(4,0)
	cycnin=new Vector(4,0)
}
proc initto(){
	tosec=new Vector(NSEC,-999)		//set/reset
	tocyc=new Vector(NSEC,-999)		//set/reset
}
proc initVectors(){
	initto()
	initsec()
	initmrg()
}
													
proc arborcharacterize(){local sec,sec1,sec2, aaa, x

	{NSEC=0  forall NSEC=NSEC+1}		
	initVectors()
	initcyc()

	sec=0
	forall {
		if(sec==0){
	          secaddress.set(0, this_section())
	          paraddress.set(0, parent_section(0))
	          secRa.set(0,Ra)
	          secL.set(0,L)

	          secCm.set(0,cm)
	          mrgL.set(0,L)
	          {for(x) if(x>0) aaa=aaa+area(x)   mrgdiam.set(0,aaa/PI/L)}
	          mrgri.set(0,secri.x(0))
		}
		if(sec!=0){
		  secaddress.set(sec, this_section())
		  paraddress.set(sec, parent_section(0))
		  {slsoma=new SectionList() rvp=new RangeVarPlot("v")
		   {soma  rvp.begin(.5)}  rvp.end(.5)  rvp.list(slsoma)}
		  {slroot=new SectionList() rvp=new RangeVarPlot("v")
		   {ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)}  ss=ss+1}}
		   rvp.end(.5)  rvp.list(slroot)}
	  	  forsec slroot order.set(sec, order.x(sec)+1)
		  for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x))
		  if(secri.x(sec)>9999) secri.set(sec, -9999)
		  forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x))
	}sec=sec+1}
	setsec()
	setmrg()
	for sec1=1,NSEC-1 {
		for sec2=0,NSEC-1 {
			if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)}
			}
		}
}
proc setsec(){
	sec=0
	forall {
		if(sec!=0){
		  secRa.set(sec,Ra)
		  secL.set(sec,L)
		  secCm.set(sec,cm)
	}sec=sec+1}
}
proc setmrg(){
	sec=0
	forall {
		if(sec!=0){
		  mrgL.set(sec,L)
		 
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
		  mrgri.set(sec,secri.x(sec))
	}sec=sec+1}
}

proc getnextABP(){local extent, sec
	extent=$1
	secA=0			
	for sec=1,NSEC-1 {if(tosec.x(sec)==-999){			
			    if(order.x(sec)>=order.x(secA)){	
			      if(abs(tocyc.x(sec))==extent){		
				secA=sec
 	}}}}
	secB=0
	for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){	
			    if(sec!=secA){
			      if(tosec.x(sec)==-999){			
				if(abs(tocyc.x(sec))==extent){		
			          secB=sec
	}}}}}
	secP=0
	if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA)
}	
proc shortnms(){
	{AmRa=0 AmL=0 Amdiam=0 Amri=0}
	{BmRa=0 BmL=0 Bmdiam=0 Bmri=0}
	{PmRa=0 PmL=0 Pmdiam=0 Pmri=0}
	if(secA>0) {
		AmRa=secRa.x(secA)
		AmL=mrgL.x(secA)  
		Amdiam=mrgdiam.x(secA) 
		Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
	}
	if(secB>0) {
		BmRa=secRa.x(secB)
		BmL=mrgL.x(secB)  
		Bmdiam=mrgdiam.x(secB) 
		Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
	}
	if(secP>0) {
		PmRa=secRa.x(secP)
		PmL=mrgL.x(secP)  
		Pmdiam=mrgdiam.x(secP) 
		Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
	}
	if(AREABYLONG==1){
	  if(secA>0 && secB>0){
	    if(AmL>BmL){
	      {BmL=AmL  mrgL.set(secB,BmL)}
	      {Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)}
	    }
	    if(BmL>AmL){
	      {AmL=BmL  mrgL.set(secA,AmL)}
	      {Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)}
	    }
	  }
	}
}










proc joinAP(){local newRa, newL, newri, newdiam
		newRa=(AmRa+PmRa)/2
		newL=AmL+PmL
		newri=Amri+Pmri
		newdiam=sqrt(newRa*newL*4/newri/PI/100)

		mrgL.set(secP, newL)
		mrgdiam.set(secP, newdiam)
		mrgri.set(secP, newri)
	tosec.set(secA, secP)
}
proc joinBP(){local newRa, newL, newri, newdiam
		newRa=(BmRa+PmRa)/2
		newL=BmL+PmL
		newri=Bmri+Pmri
		newdiam=sqrt(newRa*newL*4/newri/PI/100)

		mrgL.set(secP, newL)
		mrgdiam.set(secP, newdiam)
		mrgri.set(secP, newri)
	tosec.set(secB, secP)
}
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam

	AmL2BmL=AmL/BmL
	BmL2AmL=BmL/AmL
	if(AmL2BmL>=5) {
		joinAP()
		tosec.set(secB, -997)
	}
	if(BmL2AmL>=5) {
		joinBP()
		tosec.set(secA, -997)
	}
	if(AmL2BmL<5 && BmL2AmL<5){
		if(Lbysurf_areabyasis==1 || Lbysurf_areabylong==1) {
			//Lbysurf_areabyasis BSC_Lbysurf_areabylong
			Asurf=AmL*PI*Amdiam
			Bsurf=BmL*PI*Bmdiam
			AmBmL=(AmL*Asurf+BmL*Bsurf)/(Asurf+Bsurf)
			AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)

			AmBmcross=PI*(AmBmdiam^2)/4	
			AmBmRa=(AmRa+BmRa)/2
			AmBmri=AmBmRa*AmBmL/AmBmcross/100	

			newRa=(AmBmRa+PmRa)/2
			newL=AmBmL+PmL
			newri=AmBmri+Pmri
			newdiam=sqrt(newRa*newL*4/newri/PI/100)
		}
		if(Lbydiam_areabyasis==1 || Lbydiam_areabylong==1){
			//Lbydiam_areabyasis  Lbydiam_areabylong
			AmBmL=(AmL*Amdiam+BmL*Bmdiam)/(Amdiam+Bmdiam)
			AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
			AmBmcross=PI*(AmBmdiam^2)/4	

			AmBmRa=(AmRa+BmRa)/2
			AmBmri=AmBmRa*AmBmL/AmBmcross/100	

			newRa=(AmBmRa+PmRa)/2
			newL=AmBmL+PmL
			newri=AmBmri+Pmri
			newdiam=sqrt(newRa*newL*4/newri/PI/100)
		}
		if(Lbyavg__areabyasis==1 || Lbyavg__areabylong==1 || Lbylong_areabylong==1) {
			//Lbyavg__areabyasis Lbyavg__areabylong Lbylong_areabylong
			AmBmL=(AmL+BmL)/2
			AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
			AmBmcross=PI*(AmBmdiam^2)/4	

			AmBmRa=(AmRa+BmRa)/2
			AmBmri=AmBmRa*AmBmL/AmBmcross/100	

			newRa=(AmBmRa+PmRa)/2
			newL=AmBmL+PmL
			newri=AmBmri+Pmri
			newdiam=sqrt(newRa*newL*4/newri/PI/100)
		}
		mrgL.set(secP, newL)
		mrgri.set(secP, newri)
		mrgdiam.set(secP, newdiam)

		tosec.set(secA, secP)
		tosec.set(secB, secP)
	}
	
}
proc tagAcyc(){
	tosec.set(secA, -$1)
	tocyc.set(secA, $1)
}
proc tagBcyc(){
	tosec.set(secB, -$1)
	tocyc.set(secB, $1)
}
proc tagABcyc(){
	
	tagAcyc($1)
	tagBcyc($1)
}
proc cycmake(){local sec, extent
	initto()
	initcyc()
	setmrg()

	for sec=1,NSEC-1{
		if(secpathri.x(sec)>=CUTOFF2 && secpathri.x(sec)< 1e10) {
		  tocyc.set(sec, -1)
		}
	}
	for sec=1,NSEC-1{
		if(secpathri.x(sec)>=CUTOFF1 && secpathri.x(sec)< CUTOFF2) {
		  tocyc.set(sec, -2)
		}
	}
	for sec=1,NSEC-1{
		if(secpathri.x(sec)>=0 && secpathri.x(sec)< CUTOFF1) {
		  tocyc.set(sec, -3)
		}
	}
	for extent=1,3 {
		getnextABP(extent)
		shortnms()
		while(secA>0) {	
			if (secB==0 && secP!=0) joinAP()
			if (secB!=0 && secP!=0) joinABP()
			if (secB==0 && secP==0) tagAcyc(extent)
			if (secB!=0 && secP==0) tagABcyc(extent)

			getnextABP(extent)
			shortnms()
		}
	}

	if(Lbylong_areabylong==1){
		//Lbylong_areabylong
	  for extent=1,3 {
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, maxL)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbysurf_areabyasis==1){
		//Lbysurf_areabyasis
	  for extent=1,3 {
		totsurf=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  totsurf=totsurf+surf
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbysurf_areabylong==1){
		//Lbysurf_areabylong
	  for extent=1,3 {
		totsurf=0
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  totsurf=totsurf+surf
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbydiam_areabyasis==1){
		//Lbydiam_areabyasis
	  for extent=1,3 {
		totdiam=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  totdiam=totdiam+mrgdiam.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbydiam_areabylong==1){
		//Lbydiam_areabylong
	  for extent=1,3 {
		maxL=0
		totdiam=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  totdiam=totdiam+mrgdiam.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbyavg_areabyasis==1){
		//Lbyavg__areabyasis
	  for extent=1,3 {
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbyavg_areabylong==1){
		//Lbyavg__areabylong
	  for extent=1,3 {
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}

}

arborcharacterize()

cycmake()

print "Section\tLength (um)\tDiameter (um)"

for ii=1,3 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)


2) MY VERSION (not working)

Code: Select all



/* ---------------------------------------------------------------------

	HOC CODE TO COLLAPSE A DENDRITIC TREE
	=====================================

   Collapse a dendritic tree into 3 compartments "proximal", "middle"
   and "distal".  

   The collapse is made such as the collapsed dendritic structure
   preserves the axial resistance of the original structure. 
   The algorithm works by merging successive pairs of dendritic 
   branches into an equivalent branch (a branch that preserves the 
   axial resistance of the two original branches).
   This merging of branches can be done according to different methods
   selectable in the present code:

   AREABYLONG

	The code rejects the smallest original branch if it is more than
	5 times smaller than the long orginal branch.  If AREABYLONG is
	set to 1, then the small branch is rescaled to the longest one 
	(a new diameter is assigned such that the small branch conserves
	its axial resistance).  This minimizes errors for dendritic 
	structures where branches have a large range of possible lengths.
	AREABYLONG = 0 is the default.


   AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF

	During the merging of two branches, the length of the equivalent 
	branch is calculated by one out of three possible methods.  

	First, a simple average of the two branch length:

		L = (L1+L2)/2

	This is the method originally used by Bush and Sejnowski
	(J. Neurosci. Methods 46: 159-166, 1993).
	It is selected by setting AVGLENGTH=1.

	Second, by weighting this average with the diameters of each branch

		L = (diam1*L1 + diam2*L2)/(diam1+diam2)

	This algorithm produces more fair merging in the case two branches of
	very different diameters are to be merged.  It is selected by setting
	AVGLENGTHWDIAM=1.

	Third, by weighting the average with membrane area:

		L = (area1*L1 + area2*L2)/(area1+area2)

	This algorithm is better in the case two branches of very different
	areas are to be merged.  It is selected by setting AVGLENGTHWSURF=1.


   CUTOFF1, CUTOFF2

	The dendritic branches are assigned to one of the three compartments
	("proximal", "middle" or "distal") based on their path axial 
	resistance (in Meg-Ohms) to the soma.  The cutoff values (CUTOFF1
	and CUTOFF2) determine these cut-offs.  For example, CUTOFF1=10 means
	that all branches laying within 10 Meg-Ohm of axial resistance from 
	the soma will be merged together into the same ("proximal") section.
	

   Written by M. Neubig & A. Destexhe
   Laval University, 1997; CNRS 2000

-----------------------------------------------------------------------*/

CUTOFF1 = 7.3	// define the cutoff between "proximal" and "middle"
CUTOFF2 = 30	// define the cutoff between "middle" and "distal"
CUTOFF3 = 71

AVGLENGTH = 0		// if 1, merge based on ageraged length
AVGLENGTHWDIAM = 1	// if 1, merge based on avg length, weighted by diam
AVGLENGTHWSURF = 0	// if 1, merge based on avg length, weighted by area

AREABYLONG = 0		// if 1, rescale the length before merging



xopen("tc-cell.geo")		// read geometry file



Lbysurf_areabyasis = 0
Lbysurf_areabylong = 0
Lbydiam_areabyasis = 0
Lbydiam_areabylong = 0
Lbyavg_areabyasis = 0
Lbyavg_areabylong = 0

Lbylong_areabylong = 0

if(AREABYLONG==1) {
  print "Merge by rescaling length"
  if(AVGLENGTH == 1) {
	Lbyavg__areabylong = 1
	print "Merging method based on averaged length"
  } else if(AVGLENGTHWDIAM == 1) {
	Lbydiam_areabylong = 1
	print "Merging method based on averaged length, weigthed by diameters"
  } else if(AVGLENGTHWSURF == 1) {
	Lbysurf_areabylong = 1
  }
} else {
  print "Do not rescale length"
  if(AVGLENGTH == 1) {
	Lbyavg__areabyasis = 1
	print "Merging method based on averaged length"
  } else if(AVGLENGTHWDIAM == 1) {
	Lbydiam_areabyasis = 1
	print "Merging method based on averaged length, weigthed by diameters"
  } else if(AVGLENGTHWSURF == 1) {
	Lbysurf_areabyasis = 1
  }
}

forall nseg=10
forall Ra=260


objectvar secaddress		
objectvar paraddress		
objectvar order			
objectvar parsecndx		
				
objectvar secri			
objectvar secpathri		
objectvar secL			
objectvar secRa			
objectvar secCm			
objectvar mrgL		
objectvar mrgdiam	
objectvar mrgri			
				
objectvar tosec			
objectvar tocyc			
				
objectvar rvp			
objectvar slsoma		
objectvar slroot		

objectvar cycL
objectvar cycdiam
objectvar cycnin





proc initsec(){
	secaddress=new Vector(NSEC,0)		//set/reset
	paraddress=new Vector(NSEC,0)		//set/reset
	order=new Vector(NSEC,0)		//set/reset
	parsecndx=new Vector(NSEC,0)		//set/reset

	secri=new Vector(NSEC,0)		//set/reset
	secpathri=new Vector(NSEC,0.0)		//set/reset
	secL=new Vector(NSEC,0)   		//set/reset
	secRa=new Vector(NSEC,0)   		//set/reset
	secCm=new Vector(NSEC,0)		//set/reset
}
proc initmrg(){
	mrgL=new Vector(NSEC,0)   		//set/reset
	mrgdiam=new Vector(NSEC,0)		//set/reset
	mrgri=new Vector(NSEC,0)		//set/reset

}
proc initcyc(){
	cycL=new Vector(4,0)
	cycdiam=new Vector(4,0)
	cycnin=new Vector(4,0)
}
proc initto(){
	tosec=new Vector(NSEC,-999)		//set/reset
	tocyc=new Vector(NSEC,-999)		//set/reset
}
proc initVectors(){
	initto()
	initsec()
	initmrg()
}
													
proc arborcharacterize(){local sec,sec1,sec2, aaa, x

	{NSEC=0  forall NSEC=NSEC+1}		
	initVectors()
	initcyc()

	sec=0
	forall {
		if(sec==0){
	          secaddress.set(0, this_section())
	          paraddress.set(0, parent_section(0))
	          secRa.set(0,Ra)
	          secL.set(0,L)

	          secCm.set(0,cm)
	          mrgL.set(0,L)
	          {for(x) if(x>0) aaa=aaa+area(x)   mrgdiam.set(0,aaa/PI/L)}
	          mrgri.set(0,secri.x(0))
		}
		if(sec!=0){
		  secaddress.set(sec, this_section())
		  paraddress.set(sec, parent_section(0))
		  {slsoma=new SectionList() rvp=new RangeVarPlot("v")
		   {soma  rvp.begin(.5)}  rvp.end(.5)  rvp.list(slsoma)}
		  {slroot=new SectionList() rvp=new RangeVarPlot("v")
		   {ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)}  ss=ss+1}}
		   rvp.end(.5)  rvp.list(slroot)}
	  	  forsec slroot order.set(sec, order.x(sec)+1)
		  for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x))
		  if(secri.x(sec)>9999) secri.set(sec, -9999)
		  forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x))
	}sec=sec+1}
	setsec()
	setmrg()
	for sec1=1,NSEC-1 {
		for sec2=0,NSEC-1 {
			if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)}
			}
		}
}
proc setsec(){
	sec=0
	forall {
		if(sec!=0){
		  secRa.set(sec,Ra)
		  secL.set(sec,L)
		  secCm.set(sec,cm)
	}sec=sec+1}
}
proc setmrg(){
	sec=0
	forall {
		if(sec!=0){
		  mrgL.set(sec,L)
		 
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
		  mrgri.set(sec,secri.x(sec))
	}sec=sec+1}
}

proc getnextABP(){local extent, sec
	extent=$1
	secA=0			
	for sec=1,NSEC-1 {if(tosec.x(sec)==-999){			
			    if(order.x(sec)>=order.x(secA)){	
			      if(abs(tocyc.x(sec))==extent){		
				secA=sec
 	}}}}
	secB=0
	for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){	
			    if(sec!=secA){
			      if(tosec.x(sec)==-999){			
				if(abs(tocyc.x(sec))==extent){		
			          secB=sec
	}}}}}
	secP=0
	if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA)
}	
proc shortnms(){
	{AmRa=0 AmL=0 Amdiam=0 Amri=0}
	{BmRa=0 BmL=0 Bmdiam=0 Bmri=0}
	{PmRa=0 PmL=0 Pmdiam=0 Pmri=0}
	if(secA>0) {
		AmRa=secRa.x(secA)
		AmL=mrgL.x(secA)  
		Amdiam=mrgdiam.x(secA) 
		Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
	}
	if(secB>0) {
		BmRa=secRa.x(secB)
		BmL=mrgL.x(secB)  
		Bmdiam=mrgdiam.x(secB) 
		Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
	}
	if(secP>0) {
		PmRa=secRa.x(secP)
		PmL=mrgL.x(secP)  
		Pmdiam=mrgdiam.x(secP) 
		Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
	}
	if(AREABYLONG==1){
	  if(secA>0 && secB>0){
	    if(AmL>BmL){
	      {BmL=AmL  mrgL.set(secB,BmL)}
	      {Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)}
	    }
	    if(BmL>AmL){
	      {AmL=BmL  mrgL.set(secA,AmL)}
	      {Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)}
	    }
	  }
	}
}










proc joinAP(){local newRa, newL, newri, newdiam
		newRa=(AmRa+PmRa)/2
		newL=AmL+PmL
		newri=Amri+Pmri
		newdiam=sqrt(newRa*newL*4/newri/PI/100)

		mrgL.set(secP, newL)
		mrgdiam.set(secP, newdiam)
		mrgri.set(secP, newri)
	tosec.set(secA, secP)
}
proc joinBP(){local newRa, newL, newri, newdiam
		newRa=(BmRa+PmRa)/2
		newL=BmL+PmL
		newri=Bmri+Pmri
		newdiam=sqrt(newRa*newL*4/newri/PI/100)

		mrgL.set(secP, newL)
		mrgdiam.set(secP, newdiam)
		mrgri.set(secP, newri)
	tosec.set(secB, secP)
}
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam

	AmL2BmL=AmL/BmL
	BmL2AmL=BmL/AmL
	if(AmL2BmL>=5) {
		joinAP()
		tosec.set(secB, -997)
	}
	if(BmL2AmL>=5) {
		joinBP()
		tosec.set(secA, -997)
	}
	if(AmL2BmL<5 && BmL2AmL<5){
		if(Lbysurf_areabyasis==1 || Lbysurf_areabylong==1) {
			//Lbysurf_areabyasis BSC_Lbysurf_areabylong
			Asurf=AmL*PI*Amdiam
			Bsurf=BmL*PI*Bmdiam
			AmBmL=(AmL*Asurf+BmL*Bsurf)/(Asurf+Bsurf)
			AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)

			AmBmcross=PI*(AmBmdiam^2)/4	
			AmBmRa=(AmRa+BmRa)/2
			AmBmri=AmBmRa*AmBmL/AmBmcross/100	

			newRa=(AmBmRa+PmRa)/2
			newL=AmBmL+PmL
			newri=AmBmri+Pmri
			newdiam=sqrt(newRa*newL*4/newri/PI/100)
		}
		if(Lbydiam_areabyasis==1 || Lbydiam_areabylong==1){
			//Lbydiam_areabyasis  Lbydiam_areabylong
			AmBmL=(AmL*Amdiam+BmL*Bmdiam)/(Amdiam+Bmdiam)
			AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
			AmBmcross=PI*(AmBmdiam^2)/4	

			AmBmRa=(AmRa+BmRa)/2
			AmBmri=AmBmRa*AmBmL/AmBmcross/100	

			newRa=(AmBmRa+PmRa)/2
			newL=AmBmL+PmL
			newri=AmBmri+Pmri
			newdiam=sqrt(newRa*newL*4/newri/PI/100)
		}
		if(Lbyavg__areabyasis==1 || Lbyavg__areabylong==1 || Lbylong_areabylong==1) {
			//Lbyavg__areabyasis Lbyavg__areabylong Lbylong_areabylong
			AmBmL=(AmL+BmL)/2
			AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
			AmBmcross=PI*(AmBmdiam^2)/4	

			AmBmRa=(AmRa+BmRa)/2
			AmBmri=AmBmRa*AmBmL/AmBmcross/100	

			newRa=(AmBmRa+PmRa)/2
			newL=AmBmL+PmL
			newri=AmBmri+Pmri
			newdiam=sqrt(newRa*newL*4/newri/PI/100)
		}
		mrgL.set(secP, newL)
		mrgri.set(secP, newri)
		mrgdiam.set(secP, newdiam)

		tosec.set(secA, secP)
		tosec.set(secB, secP)
	}
	
}
proc tagAcyc(){
	tosec.set(secA, -$1)
	tocyc.set(secA, $1)
}
proc tagBcyc(){
	tosec.set(secB, -$1)
	tocyc.set(secB, $1)
}
proc tagABcyc(){
	
	tagAcyc($1)
	tagBcyc($1)
}
proc cycmake(){local sec, extent
	initto()
	initcyc()
	setmrg()

	for sec=1,NSEC-1{
		if(secpathri.x(sec)>=CUTOFF3 && secpathri.x(sec)< 1e10) {
		  tocyc.set(sec, -1)
		}
	}
	for sec=1,NSEC-1{
		if(secpathri.x(sec)>=CUTOFF2 && secpathri.x(sec)< 1e10) {
		  tocyc.set(sec, -2)
		}
	}
	for sec=1,NSEC-1{
		if(secpathri.x(sec)>=CUTOFF1 && secpathri.x(sec)< CUTOFF2) {
		  tocyc.set(sec, -3)
		}
	}
	for sec=1,NSEC-1{
		if(secpathri.x(sec)>=0 && secpathri.x(sec)< CUTOFF1) {
		  tocyc.set(sec, -4)
		}
	}
	for extent=1,4 {
		getnextABP(extent)
		shortnms()
		while(secA>0) {	
			if (secB==0 && secP!=0) joinAP()
			if (secB!=0 && secP!=0) joinABP()
			if (secB==0 && secP==0) tagAcyc(extent)
			if (secB!=0 && secP==0) tagABcyc(extent)

			getnextABP(extent)
			shortnms()
		}
	}

	if(Lbylong_areabylong==1){
		//Lbylong_areabylong
	  for extent=1,4 {
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, maxL)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbysurf_areabyasis==1){
		//Lbysurf_areabyasis
	  for extent=1,4 {
		totsurf=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  totsurf=totsurf+surf
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbysurf_areabylong==1){
		//Lbysurf_areabylong
	  for extent=1,4 {
		totsurf=0
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  totsurf=totsurf+surf
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbydiam_areabyasis==1){
		//Lbydiam_areabyasis
	  for extent=1,4 {
		totdiam=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  totdiam=totdiam+mrgdiam.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbydiam_areabylong==1){
		//Lbydiam_areabylong
	  for extent=1,4 {
		maxL=0
		totdiam=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  totdiam=totdiam+mrgdiam.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbyavg_areabyasis==1){
		//Lbyavg__areabyasis
	  for extent=1,4 {
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbyavg_areabylong==1){
		//Lbyavg__areabylong
	  for extent=1,4 {
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}

}

arborcharacterize()

cycmake()

print "Section\tLength (um)\tDiameter (um)"

for ii=1,4 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)

Mike Neubig
Posts: 2
Joined: Tue May 23, 2006 8:39 pm

Post by Mike Neubig »

in cycmake(), this line looks suspicious:

if(secpathri.x(sec)>=CUTOFF2 && secpathri.x(sec)< 1e10) {

it probably needs to be something like:

if(secpathri.x(sec)>=CUTOFF2 && secpathri.x(sec)< CUTOFF3) {

==========

...... but i would suggest approaching your revisions in a new way:

1) make sure original code works okay for a variety of parameters

2) refactor the original code to use a Vector for CUTOFFs and use
the length of that vector when rewriting lines such as:

for extent=1,4 {

and in other places as called for.

3) in each small step of the refactoring process, make sure that, when n = 3, the modified code produces results that are identical to the original code.

4) after all refactoring is completed, debug for n != 3

http://en.wikipedia.org/wiki/Refactoring

best,
mike
mikey
Posts: 37
Joined: Mon Jun 19, 2006 3:06 pm

Post by mikey »

thanks mike so much for your kind reply. I have changed the code to incorporate the line change you suggested. Unfortunately the code still doesn't work. dam!

I understand the alternative approach you suggested is superiour but I'm afraid I don't think I have the programming capability to do it. I am really only just hanging on by my finger tips doing these line changes as we are now. If I can get this code working for 4 as opposed to 3 - I think I will know how to just flesh it out manually (using cuts and pastes and little edits etc.) for 40 compartments. Then I will be very, very happy. Would you be kind enough to take a look at the newest, and still non-working, version of the code I have cut and paste below? I would be so, so grateful.

Code: Select all


/* ---------------------------------------------------------------------

	HOC CODE TO COLLAPSE A DENDRITIC TREE
	=====================================

   Collapse a dendritic tree into 3 compartments "proximal", "middle"
   and "distal".  

   The collapse is made such as the collapsed dendritic structure
   preserves the axial resistance of the original structure. 
   The algorithm works by merging successive pairs of dendritic 
   branches into an equivalent branch (a branch that preserves the 
   axial resistance of the two original branches).
   This merging of branches can be done according to different methods
   selectable in the present code:

   AREABYLONG

	The code rejects the smallest original branch if it is more than
	5 times smaller than the long orginal branch.  If AREABYLONG is
	set to 1, then the small branch is rescaled to the longest one 
	(a new diameter is assigned such that the small branch conserves
	its axial resistance).  This minimizes errors for dendritic 
	structures where branches have a large range of possible lengths.
	AREABYLONG = 0 is the default.


   AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF

	During the merging of two branches, the length of the equivalent 
	branch is calculated by one out of three possible methods.  

	First, a simple average of the two branch length:

		L = (L1+L2)/2

	This is the method originally used by Bush and Sejnowski
	(J. Neurosci. Methods 46: 159-166, 1993).
	It is selected by setting AVGLENGTH=1.

	Second, by weighting this average with the diameters of each branch

		L = (diam1*L1 + diam2*L2)/(diam1+diam2)

	This algorithm produces more fair merging in the case two branches of
	very different diameters are to be merged.  It is selected by setting
	AVGLENGTHWDIAM=1.

	Third, by weighting the average with membrane area:

		L = (area1*L1 + area2*L2)/(area1+area2)

	This algorithm is better in the case two branches of very different
	areas are to be merged.  It is selected by setting AVGLENGTHWSURF=1.


   CUTOFF1, CUTOFF2

	The dendritic branches are assigned to one of the three compartments
	("proximal", "middle" or "distal") based on their path axial 
	resistance (in Meg-Ohms) to the soma.  The cutoff values (CUTOFF1
	and CUTOFF2) determine these cut-offs.  For example, CUTOFF1=10 means
	that all branches laying within 10 Meg-Ohm of axial resistance from 
	the soma will be merged together into the same ("proximal") section.
	

   Written by M. Neubig & A. Destexhe
   Laval University, 1997; CNRS 2000

-----------------------------------------------------------------------*/

CUTOFF1 = 7.3	// define the cutoff between "proximal" and "middle"
CUTOFF2 = 30	// define the cutoff between "middle" and "distal"
CUTOFF3 = 71

AVGLENGTH = 0		// if 1, merge based on ageraged length
AVGLENGTHWDIAM = 1	// if 1, merge based on avg length, weighted by diam
AVGLENGTHWSURF = 0	// if 1, merge based on avg length, weighted by area

AREABYLONG = 0		// if 1, rescale the length before merging



xopen("tc-cell.geo")		// read geometry file



Lbysurf_areabyasis = 0
Lbysurf_areabylong = 0
Lbydiam_areabyasis = 0
Lbydiam_areabylong = 0
Lbyavg_areabyasis = 0
Lbyavg_areabylong = 0

Lbylong_areabylong = 0

if(AREABYLONG==1) {
  print "Merge by rescaling length"
  if(AVGLENGTH == 1) {
	Lbyavg__areabylong = 1
	print "Merging method based on averaged length"
  } else if(AVGLENGTHWDIAM == 1) {
	Lbydiam_areabylong = 1
	print "Merging method based on averaged length, weigthed by diameters"
  } else if(AVGLENGTHWSURF == 1) {
	Lbysurf_areabylong = 1
  }
} else {
  print "Do not rescale length"
  if(AVGLENGTH == 1) {
	Lbyavg__areabyasis = 1
	print "Merging method based on averaged length"
  } else if(AVGLENGTHWDIAM == 1) {
	Lbydiam_areabyasis = 1
	print "Merging method based on averaged length, weigthed by diameters"
  } else if(AVGLENGTHWSURF == 1) {
	Lbysurf_areabyasis = 1
  }
}

forall nseg=10
forall Ra=260


objectvar secaddress		
objectvar paraddress		
objectvar order			
objectvar parsecndx		
				
objectvar secri			
objectvar secpathri		
objectvar secL			
objectvar secRa			
objectvar secCm			
objectvar mrgL		
objectvar mrgdiam	
objectvar mrgri			
				
objectvar tosec			
objectvar tocyc			
				
objectvar rvp			
objectvar slsoma		
objectvar slroot		

objectvar cycL
objectvar cycdiam
objectvar cycnin





proc initsec(){
	secaddress=new Vector(NSEC,0)		//set/reset
	paraddress=new Vector(NSEC,0)		//set/reset
	order=new Vector(NSEC,0)		//set/reset
	parsecndx=new Vector(NSEC,0)		//set/reset

	secri=new Vector(NSEC,0)		//set/reset
	secpathri=new Vector(NSEC,0.0)		//set/reset
	secL=new Vector(NSEC,0)   		//set/reset
	secRa=new Vector(NSEC,0)   		//set/reset
	secCm=new Vector(NSEC,0)		//set/reset
}
proc initmrg(){
	mrgL=new Vector(NSEC,0)   		//set/reset
	mrgdiam=new Vector(NSEC,0)		//set/reset
	mrgri=new Vector(NSEC,0)		//set/reset

}
proc initcyc(){
	cycL=new Vector(4,0)
	cycdiam=new Vector(4,0)
	cycnin=new Vector(4,0)
}
proc initto(){
	tosec=new Vector(NSEC,-999)		//set/reset
	tocyc=new Vector(NSEC,-999)		//set/reset
}
proc initVectors(){
	initto()
	initsec()
	initmrg()
}
													
proc arborcharacterize(){local sec,sec1,sec2, aaa, x

	{NSEC=0  forall NSEC=NSEC+1}		
	initVectors()
	initcyc()

	sec=0
	forall {
		if(sec==0){
	          secaddress.set(0, this_section())
	          paraddress.set(0, parent_section(0))
	          secRa.set(0,Ra)
	          secL.set(0,L)

	          secCm.set(0,cm)
	          mrgL.set(0,L)
	          {for(x) if(x>0) aaa=aaa+area(x)   mrgdiam.set(0,aaa/PI/L)}
	          mrgri.set(0,secri.x(0))
		}
		if(sec!=0){
		  secaddress.set(sec, this_section())
		  paraddress.set(sec, parent_section(0))
		  {slsoma=new SectionList() rvp=new RangeVarPlot("v")
		   {soma  rvp.begin(.5)}  rvp.end(.5)  rvp.list(slsoma)}
		  {slroot=new SectionList() rvp=new RangeVarPlot("v")
		   {ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)}  ss=ss+1}}
		   rvp.end(.5)  rvp.list(slroot)}
	  	  forsec slroot order.set(sec, order.x(sec)+1)
		  for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x))
		  if(secri.x(sec)>9999) secri.set(sec, -9999)
		  forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x))
	}sec=sec+1}
	setsec()
	setmrg()
	for sec1=1,NSEC-1 {
		for sec2=0,NSEC-1 {
			if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)}
			}
		}
}
proc setsec(){
	sec=0
	forall {
		if(sec!=0){
		  secRa.set(sec,Ra)
		  secL.set(sec,L)
		  secCm.set(sec,cm)
	}sec=sec+1}
}
proc setmrg(){
	sec=0
	forall {
		if(sec!=0){
		  mrgL.set(sec,L)
		 
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
		  mrgri.set(sec,secri.x(sec))
	}sec=sec+1}
}

proc getnextABP(){local extent, sec
	extent=$1
	secA=0			
	for sec=1,NSEC-1 {if(tosec.x(sec)==-999){			
			    if(order.x(sec)>=order.x(secA)){	
			      if(abs(tocyc.x(sec))==extent){		
				secA=sec
 	}}}}
	secB=0
	for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){	
			    if(sec!=secA){
			      if(tosec.x(sec)==-999){			
				if(abs(tocyc.x(sec))==extent){		
			          secB=sec
	}}}}}
	secP=0
	if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA)
}	
proc shortnms(){
	{AmRa=0 AmL=0 Amdiam=0 Amri=0}
	{BmRa=0 BmL=0 Bmdiam=0 Bmri=0}
	{PmRa=0 PmL=0 Pmdiam=0 Pmri=0}
	if(secA>0) {
		AmRa=secRa.x(secA)
		AmL=mrgL.x(secA)  
		Amdiam=mrgdiam.x(secA) 
		Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
	}
	if(secB>0) {
		BmRa=secRa.x(secB)
		BmL=mrgL.x(secB)  
		Bmdiam=mrgdiam.x(secB) 
		Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
	}
	if(secP>0) {
		PmRa=secRa.x(secP)
		PmL=mrgL.x(secP)  
		Pmdiam=mrgdiam.x(secP) 
		Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
	}
	if(AREABYLONG==1){
	  if(secA>0 && secB>0){
	    if(AmL>BmL){
	      {BmL=AmL  mrgL.set(secB,BmL)}
	      {Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)}
	    }
	    if(BmL>AmL){
	      {AmL=BmL  mrgL.set(secA,AmL)}
	      {Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)}
	    }
	  }
	}
}










proc joinAP(){local newRa, newL, newri, newdiam
		newRa=(AmRa+PmRa)/2
		newL=AmL+PmL
		newri=Amri+Pmri
		newdiam=sqrt(newRa*newL*4/newri/PI/100)

		mrgL.set(secP, newL)
		mrgdiam.set(secP, newdiam)
		mrgri.set(secP, newri)
	tosec.set(secA, secP)
}
proc joinBP(){local newRa, newL, newri, newdiam
		newRa=(BmRa+PmRa)/2
		newL=BmL+PmL
		newri=Bmri+Pmri
		newdiam=sqrt(newRa*newL*4/newri/PI/100)

		mrgL.set(secP, newL)
		mrgdiam.set(secP, newdiam)
		mrgri.set(secP, newri)
	tosec.set(secB, secP)
}
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam

	AmL2BmL=AmL/BmL
	BmL2AmL=BmL/AmL
	if(AmL2BmL>=5) {
		joinAP()
		tosec.set(secB, -997)
	}
	if(BmL2AmL>=5) {
		joinBP()
		tosec.set(secA, -997)
	}
	if(AmL2BmL<5 && BmL2AmL<5){
		if(Lbysurf_areabyasis==1 || Lbysurf_areabylong==1) {
			//Lbysurf_areabyasis BSC_Lbysurf_areabylong
			Asurf=AmL*PI*Amdiam
			Bsurf=BmL*PI*Bmdiam
			AmBmL=(AmL*Asurf+BmL*Bsurf)/(Asurf+Bsurf)
			AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)

			AmBmcross=PI*(AmBmdiam^2)/4	
			AmBmRa=(AmRa+BmRa)/2
			AmBmri=AmBmRa*AmBmL/AmBmcross/100	

			newRa=(AmBmRa+PmRa)/2
			newL=AmBmL+PmL
			newri=AmBmri+Pmri
			newdiam=sqrt(newRa*newL*4/newri/PI/100)
		}
		if(Lbydiam_areabyasis==1 || Lbydiam_areabylong==1){
			//Lbydiam_areabyasis  Lbydiam_areabylong
			AmBmL=(AmL*Amdiam+BmL*Bmdiam)/(Amdiam+Bmdiam)
			AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
			AmBmcross=PI*(AmBmdiam^2)/4	

			AmBmRa=(AmRa+BmRa)/2
			AmBmri=AmBmRa*AmBmL/AmBmcross/100	

			newRa=(AmBmRa+PmRa)/2
			newL=AmBmL+PmL
			newri=AmBmri+Pmri
			newdiam=sqrt(newRa*newL*4/newri/PI/100)
		}
		if(Lbyavg__areabyasis==1 || Lbyavg__areabylong==1 || Lbylong_areabylong==1) {
			//Lbyavg__areabyasis Lbyavg__areabylong Lbylong_areabylong
			AmBmL=(AmL+BmL)/2
			AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
			AmBmcross=PI*(AmBmdiam^2)/4	

			AmBmRa=(AmRa+BmRa)/2
			AmBmri=AmBmRa*AmBmL/AmBmcross/100	

			newRa=(AmBmRa+PmRa)/2
			newL=AmBmL+PmL
			newri=AmBmri+Pmri
			newdiam=sqrt(newRa*newL*4/newri/PI/100)
		}
		mrgL.set(secP, newL)
		mrgri.set(secP, newri)
		mrgdiam.set(secP, newdiam)

		tosec.set(secA, secP)
		tosec.set(secB, secP)
	}
	
}
proc tagAcyc(){
	tosec.set(secA, -$1)
	tocyc.set(secA, $1)
}
proc tagBcyc(){
	tosec.set(secB, -$1)
	tocyc.set(secB, $1)
}
proc tagABcyc(){
	
	tagAcyc($1)
	tagBcyc($1)
}
proc cycmake(){local sec, extent
	initto()
	initcyc()
	setmrg()
		
	for sec=1,NSEC-1{
		if(secpathri.x(sec)>=CUTOFF3 && secpathri.x(sec)< 1e10) {
		  tocyc.set(sec, -1)
		}
	}
	for sec=1,NSEC-1{
		if(secpathri.x(sec)>=CUTOFF2 && secpathri.x(sec)< CUTOFF3) {
		  tocyc.set(sec, -2)
		}
	}
	for sec=1,NSEC-1{
		if(secpathri.x(sec)>=CUTOFF1 && secpathri.x(sec)< CUTOFF2) {
		  tocyc.set(sec, -3)
		}
	}
	for sec=1,NSEC-1{
		if(secpathri.x(sec)>=0 && secpathri.x(sec)< CUTOFF1) {
		  tocyc.set(sec, -4)
		}
	}
	for extent=1,4 {
		getnextABP(extent)
		shortnms()
		while(secA>0) {	
			if (secB==0 && secP!=0) joinAP()
			if (secB!=0 && secP!=0) joinABP()
			if (secB==0 && secP==0) tagAcyc(extent)
			if (secB!=0 && secP==0) tagABcyc(extent)

			getnextABP(extent)
			shortnms()
		}
	}

	if(Lbylong_areabylong==1){
		//Lbylong_areabylong
	  for extent=1,4 {
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, maxL)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbysurf_areabyasis==1){
		//Lbysurf_areabyasis
	  for extent=1,4 {
		totsurf=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  totsurf=totsurf+surf
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbysurf_areabylong==1){
		//Lbysurf_areabylong
	  for extent=1,4 {
		totsurf=0
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  totsurf=totsurf+surf
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbydiam_areabyasis==1){
		//Lbydiam_areabyasis
	  for extent=1,4 {
		totdiam=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  totdiam=totdiam+mrgdiam.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbydiam_areabylong==1){
		//Lbydiam_areabylong
	  for extent=1,4 {
		maxL=0
		totdiam=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  totdiam=totdiam+mrgdiam.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbyavg_areabyasis==1){
		//Lbyavg__areabyasis
	  for extent=1,4 {
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbyavg_areabylong==1){
		//Lbyavg__areabylong
	  for extent=1,4 {
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}

}

arborcharacterize()

cycmake()

print "Section\tLength (um)\tDiameter (um)"

for ii=1,4 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)


Mike Neubig
Posts: 2
Joined: Tue May 23, 2006 8:39 pm

Post by Mike Neubig »

the distributed nature of a dendritic arbor creates an uncoupled environment for synapses and membrane channels. collapsing a dendritic arbor into one with fewer branches, or a single tree, or a long cyclinder couples them.
mikey
Posts: 37
Joined: Mon Jun 19, 2006 3:06 pm

Post by mikey »

Mike Neubig wrote:the distributed nature of a dendritic arbor creates an uncoupled environment for synapses and membrane channels. collapsing a dendritic arbor into one with fewer branches, or a single tree, or a long cyclinder couples them.
Im not putting synapses onto my dendritic tree so this isn't so much of a worry for me. Am still really, really stuck with the coding. Haven't got past the non-working coding iteration in my last post. Any help would really make my weekend and would be so gratefully received.
mikey
Posts: 37
Joined: Mon Jun 19, 2006 3:06 pm

Post by mikey »

Thanks so much to Dimitris who gave me this tip:
---------
Are you sure these vectors have the correct length?
proc initcyc(){
cycL=new Vector(4,0)
cycdiam=new Vector(4,0)
cycnin=new Vector(4,0)
}
---------
This got me to a version of the code working for 4 compartments as opposed to 3. I have cut and paste this working code below. HOWEVER, when I then tried to further modify the code to yield 10 compartments as opposed to 4 it didnt work. I have cut and paste this code below as well (in the next post). Stumped again. Dam! Would be so, so grateful for another pair of eyes. There must be a stupid mistake in there.

1) WORKING CODE FOR 4 COMPARTMENTS

Code: Select all


/* ---------------------------------------------------------------------

	HOC CODE TO COLLAPSE A DENDRITIC TREE
	=====================================

   Collapse a dendritic tree into 3 compartments "proximal", "middle"
   and "distal".  

   The collapse is made such as the collapsed dendritic structure
   preserves the axial resistance of the original structure. 
   The algorithm works by merging successive pairs of dendritic 
   branches into an equivalent branch (a branch that preserves the 
   axial resistance of the two original branches).
   This merging of branches can be done according to different methods
   selectable in the present code:

   AREABYLONG

	The code rejects the smallest original branch if it is more than
	5 times smaller than the long orginal branch.  If AREABYLONG is
	set to 1, then the small branch is rescaled to the longest one 
	(a new diameter is assigned such that the small branch conserves
	its axial resistance).  This minimizes errors for dendritic 
	structures where branches have a large range of possible lengths.
	AREABYLONG = 0 is the default.


   AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF

	During the merging of two branches, the length of the equivalent 
	branch is calculated by one out of three possible methods.  

	First, a simple average of the two branch length:

		L = (L1+L2)/2

	This is the method originally used by Bush and Sejnowski
	(J. Neurosci. Methods 46: 159-166, 1993).
	It is selected by setting AVGLENGTH=1.

	Second, by weighting this average with the diameters of each branch

		L = (diam1*L1 + diam2*L2)/(diam1+diam2)

	This algorithm produces more fair merging in the case two branches of
	very different diameters are to be merged.  It is selected by setting
	AVGLENGTHWDIAM=1.

	Third, by weighting the average with membrane area:

		L = (area1*L1 + area2*L2)/(area1+area2)

	This algorithm is better in the case two branches of very different
	areas are to be merged.  It is selected by setting AVGLENGTHWSURF=1.


   CUTOFF1, CUTOFF2

	The dendritic branches are assigned to one of the three compartments
	("proximal", "middle" or "distal") based on their path axial 
	resistance (in Meg-Ohms) to the soma.  The cutoff values (CUTOFF1
	and CUTOFF2) determine these cut-offs.  For example, CUTOFF1=10 means
	that all branches laying within 10 Meg-Ohm of axial resistance from 
	the soma will be merged together into the same ("proximal") section.
	

   Written by M. Neubig & A. Destexhe
   Laval University, 1997; CNRS 2000

-----------------------------------------------------------------------*/

CUTOFF1 = 7.3	// define the cutoff between "proximal" and "middle"
CUTOFF2 = 30	// define the cutoff between "middle" and "distal"
CUTOFF3 = 71

AVGLENGTH = 0		// if 1, merge based on ageraged length
AVGLENGTHWDIAM = 1	// if 1, merge based on avg length, weighted by diam
AVGLENGTHWSURF = 0	// if 1, merge based on avg length, weighted by area

AREABYLONG = 0		// if 1, rescale the length before merging



xopen("tc-cell.geo")		// read geometry file



Lbysurf_areabyasis = 0
Lbysurf_areabylong = 0
Lbydiam_areabyasis = 0
Lbydiam_areabylong = 0
Lbyavg_areabyasis = 0
Lbyavg_areabylong = 0

Lbylong_areabylong = 0

if(AREABYLONG==1) {
  print "Merge by rescaling length"
  if(AVGLENGTH == 1) {
	Lbyavg__areabylong = 1
	print "Merging method based on averaged length"
  } else if(AVGLENGTHWDIAM == 1) {
	Lbydiam_areabylong = 1
	print "Merging method based on averaged length, weigthed by diameters"
  } else if(AVGLENGTHWSURF == 1) {
	Lbysurf_areabylong = 1
  }
} else {
  print "Do not rescale length"
  if(AVGLENGTH == 1) {
	Lbyavg__areabyasis = 1
	print "Merging method based on averaged length"
  } else if(AVGLENGTHWDIAM == 1) {
	Lbydiam_areabyasis = 1
	print "Merging method based on averaged length, weigthed by diameters"
  } else if(AVGLENGTHWSURF == 1) {
	Lbysurf_areabyasis = 1
  }
}

forall nseg=10
forall Ra=260


objectvar secaddress		
objectvar paraddress		
objectvar order			
objectvar parsecndx		
				
objectvar secri			
objectvar secpathri		
objectvar secL			
objectvar secRa			
objectvar secCm			
objectvar mrgL		
objectvar mrgdiam	
objectvar mrgri			
				
objectvar tosec			
objectvar tocyc			
				
objectvar rvp			
objectvar slsoma		
objectvar slroot		

objectvar cycL
objectvar cycdiam
objectvar cycnin





proc initsec(){
	secaddress=new Vector(NSEC,0)		//set/reset
	paraddress=new Vector(NSEC,0)		//set/reset
	order=new Vector(NSEC,0)		//set/reset
	parsecndx=new Vector(NSEC,0)		//set/reset

	secri=new Vector(NSEC,0)		//set/reset
	secpathri=new Vector(NSEC,0.0)		//set/reset
	secL=new Vector(NSEC,0)   		//set/reset
	secRa=new Vector(NSEC,0)   		//set/reset
	secCm=new Vector(NSEC,0)		//set/reset
}
proc initmrg(){
	mrgL=new Vector(NSEC,0)   		//set/reset
	mrgdiam=new Vector(NSEC,0)		//set/reset
	mrgri=new Vector(NSEC,0)		//set/reset

}
proc initcyc(){
	cycL=new Vector(5,0)
	cycdiam=new Vector(5,0)
	cycnin=new Vector(5,0)
}
proc initto(){
	tosec=new Vector(NSEC,-999)		//set/reset
	tocyc=new Vector(NSEC,-999)		//set/reset
}
proc initVectors(){
	initto()
	initsec()
	initmrg()
}
													
proc arborcharacterize(){local sec,sec1,sec2, aaa, x

	{NSEC=0  forall NSEC=NSEC+1}		
	initVectors()
	initcyc()

	sec=0
	forall {
		if(sec==0){
	          secaddress.set(0, this_section())
	          paraddress.set(0, parent_section(0))
	          secRa.set(0,Ra)
	          secL.set(0,L)

	          secCm.set(0,cm)
	          mrgL.set(0,L)
	          {for(x) if(x>0) aaa=aaa+area(x)   mrgdiam.set(0,aaa/PI/L)}
	          mrgri.set(0,secri.x(0))
		}
		if(sec!=0){
		  secaddress.set(sec, this_section())
		  paraddress.set(sec, parent_section(0))
		  {slsoma=new SectionList() rvp=new RangeVarPlot("v")
		   {soma  rvp.begin(.5)}  rvp.end(.5)  rvp.list(slsoma)}
		  {slroot=new SectionList() rvp=new RangeVarPlot("v")
		   {ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)}  ss=ss+1}}
		   rvp.end(.5)  rvp.list(slroot)}
	  	  forsec slroot order.set(sec, order.x(sec)+1)
		  for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x))
		  if(secri.x(sec)>9999) secri.set(sec, -9999)
		  forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x))
	}sec=sec+1}
	setsec()
	setmrg()
	for sec1=1,NSEC-1 {
		for sec2=0,NSEC-1 {
			if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)}
			}
		}
}
proc setsec(){
	sec=0
	forall {
		if(sec!=0){
		  secRa.set(sec,Ra)
		  secL.set(sec,L)
		  secCm.set(sec,cm)
	}sec=sec+1}
}
proc setmrg(){
	sec=0
	forall {
		if(sec!=0){
		  mrgL.set(sec,L)
		 
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
		  mrgri.set(sec,secri.x(sec))
	}sec=sec+1}
}

proc getnextABP(){local extent, sec
	extent=$1
	secA=0			
	for sec=1,NSEC-1 {if(tosec.x(sec)==-999){			
			    if(order.x(sec)>=order.x(secA)){	
			      if(abs(tocyc.x(sec))==extent){		
				secA=sec
 	}}}}
	secB=0
	for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){	
			    if(sec!=secA){
			      if(tosec.x(sec)==-999){			
				if(abs(tocyc.x(sec))==extent){		
			          secB=sec
	}}}}}
	secP=0
	if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA)
}	
proc shortnms(){
	{AmRa=0 AmL=0 Amdiam=0 Amri=0}
	{BmRa=0 BmL=0 Bmdiam=0 Bmri=0}
	{PmRa=0 PmL=0 Pmdiam=0 Pmri=0}
	if(secA>0) {
		AmRa=secRa.x(secA)
		AmL=mrgL.x(secA)  
		Amdiam=mrgdiam.x(secA) 
		Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
	}
	if(secB>0) {
		BmRa=secRa.x(secB)
		BmL=mrgL.x(secB)  
		Bmdiam=mrgdiam.x(secB) 
		Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
	}
	if(secP>0) {
		PmRa=secRa.x(secP)
		PmL=mrgL.x(secP)  
		Pmdiam=mrgdiam.x(secP) 
		Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
	}
	if(AREABYLONG==1){
	  if(secA>0 && secB>0){
	    if(AmL>BmL){
	      {BmL=AmL  mrgL.set(secB,BmL)}
	      {Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)}
	    }
	    if(BmL>AmL){
	      {AmL=BmL  mrgL.set(secA,AmL)}
	      {Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)}
	    }
	  }
	}
}










proc joinAP(){local newRa, newL, newri, newdiam
		newRa=(AmRa+PmRa)/2
		newL=AmL+PmL
		newri=Amri+Pmri
		newdiam=sqrt(newRa*newL*4/newri/PI/100)

		mrgL.set(secP, newL)
		mrgdiam.set(secP, newdiam)
		mrgri.set(secP, newri)
	tosec.set(secA, secP)
}
proc joinBP(){local newRa, newL, newri, newdiam
		newRa=(BmRa+PmRa)/2
		newL=BmL+PmL
		newri=Bmri+Pmri
		newdiam=sqrt(newRa*newL*4/newri/PI/100)

		mrgL.set(secP, newL)
		mrgdiam.set(secP, newdiam)
		mrgri.set(secP, newri)
	tosec.set(secB, secP)
}
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam

	AmL2BmL=AmL/BmL
	BmL2AmL=BmL/AmL
	if(AmL2BmL>=5) {
		joinAP()
		tosec.set(secB, -997)
	}
	if(BmL2AmL>=5) {
		joinBP()
		tosec.set(secA, -997)
	}
	if(AmL2BmL<5 && BmL2AmL<5>=CUTOFF3 && secpathri.x(sec)<1e10>=CUTOFF2 && secpathri.x(sec)<CUTOFF3>=CUTOFF1 && secpathri.x(sec)<CUTOFF2>=0 && secpathri.x(sec)<CUTOFF1>0) {	
			if (secB==0 && secP!=0) joinAP()
			if (secB!=0 && secP!=0) joinABP()
			if (secB==0 && secP==0) tagAcyc(extent)
			if (secB!=0 && secP==0) tagABcyc(extent)

			getnextABP(extent)
			shortnms()
		}
	}

	if(Lbylong_areabylong==1){
		//Lbylong_areabylong
	  for extent=1,4 {
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, maxL)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbysurf_areabyasis==1){
		//Lbysurf_areabyasis
	  for extent=1,4 {
		totsurf=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  totsurf=totsurf+surf
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbysurf_areabylong==1){
		//Lbysurf_areabylong
	  for extent=1,4 {
		totsurf=0
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  totsurf=totsurf+surf
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbydiam_areabyasis==1){
		//Lbydiam_areabyasis
	  for extent=1,4 {
		totdiam=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  totdiam=totdiam+mrgdiam.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbydiam_areabylong==1){
		//Lbydiam_areabylong
	  for extent=1,4 {
		maxL=0
		totdiam=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  totdiam=totdiam+mrgdiam.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbyavg_areabyasis==1){
		//Lbyavg__areabyasis
	  for extent=1,4 {
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbyavg_areabylong==1){
		//Lbyavg__areabylong
	  for extent=1,4 {
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}

}

arborcharacterize()

cycmake()

print "Section\tLength (um)\tDiameter (um)"

for ii=1,4 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)


mikey
Posts: 37
Joined: Mon Jun 19, 2006 3:06 pm

Post by mikey »

2) NON-WORKING CODE FOR 10 COMPARTMENTS
NON_WORKING - HELP!!!!!!

Code: Select all


/* ---------------------------------------------------------------------

	HOC CODE TO COLLAPSE A DENDRITIC TREE
	=====================================

   Collapse a dendritic tree into 3 compartments "proximal", "middle"
   and "distal".  

   The collapse is made such as the collapsed dendritic structure
   preserves the axial resistance of the original structure. 
   The algorithm works by merging successive pairs of dendritic 
   branches into an equivalent branch (a branch that preserves the 
   axial resistance of the two original branches).
   This merging of branches can be done according to different methods
   selectable in the present code:

   AREABYLONG

	The code rejects the smallest original branch if it is more than
	5 times smaller than the long orginal branch.  If AREABYLONG is
	set to 1, then the small branch is rescaled to the longest one 
	(a new diameter is assigned such that the small branch conserves
	its axial resistance).  This minimizes errors for dendritic 
	structures where branches have a large range of possible lengths.
	AREABYLONG = 0 is the default.


   AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF

	During the merging of two branches, the length of the equivalent 
	branch is calculated by one out of three possible methods.  

	First, a simple average of the two branch length:

		L = (L1+L2)/2

	This is the method originally used by Bush and Sejnowski
	(J. Neurosci. Methods 46: 159-166, 1993).
	It is selected by setting AVGLENGTH=1.

	Second, by weighting this average with the diameters of each branch

		L = (diam1*L1 + diam2*L2)/(diam1+diam2)

	This algorithm produces more fair merging in the case two branches of
	very different diameters are to be merged.  It is selected by setting
	AVGLENGTHWDIAM=1.

	Third, by weighting the average with membrane area:

		L = (area1*L1 + area2*L2)/(area1+area2)

	This algorithm is better in the case two branches of very different
	areas are to be merged.  It is selected by setting AVGLENGTHWSURF=1.


   CUTOFF1, CUTOFF2

	The dendritic branches are assigned to one of the three compartments
	("proximal", "middle" or "distal") based on their path axial 
	resistance (in Meg-Ohms) to the soma.  The cutoff values (CUTOFF1
	and CUTOFF2) determine these cut-offs.  For example, CUTOFF1=10 means
	that all branches laying within 10 Meg-Ohm of axial resistance from 
	the soma will be merged together into the same ("proximal") section.
	

   Written by M. Neubig & A. Destexhe
   Laval University, 1997; CNRS 2000

-----------------------------------------------------------------------*/

CUTOFF1 = 7.3	// define the cutoff between "proximal" and "middle"
CUTOFF2 = 10	// define the cutoff between "middle" and "distal"
CUTOFF3 = 20
CUTOFF4 = 30	
CUTOFF5 = 40	
CUTOFF6 = 50
CUTOFF7 = 60	
CUTOFF8 = 65	
CUTOFF9 = 71

AVGLENGTH = 0		// if 1, merge based on ageraged length
AVGLENGTHWDIAM = 1	// if 1, merge based on avg length, weighted by diam
AVGLENGTHWSURF = 0	// if 1, merge based on avg length, weighted by area

AREABYLONG = 0		// if 1, rescale the length before merging



xopen("tc-cell.geo")		// read geometry file



Lbysurf_areabyasis = 0
Lbysurf_areabylong = 0
Lbydiam_areabyasis = 0
Lbydiam_areabylong = 0
Lbyavg_areabyasis = 0
Lbyavg_areabylong = 0

Lbylong_areabylong = 0

if(AREABYLONG==1) {
  print "Merge by rescaling length"
  if(AVGLENGTH == 1) {
	Lbyavg__areabylong = 1
	print "Merging method based on averaged length"
  } else if(AVGLENGTHWDIAM == 1) {
	Lbydiam_areabylong = 1
	print "Merging method based on averaged length, weigthed by diameters"
  } else if(AVGLENGTHWSURF == 1) {
	Lbysurf_areabylong = 1
  }
} else {
  print "Do not rescale length"
  if(AVGLENGTH == 1) {
	Lbyavg__areabyasis = 1
	print "Merging method based on averaged length"
  } else if(AVGLENGTHWDIAM == 1) {
	Lbydiam_areabyasis = 1
	print "Merging method based on averaged length, weigthed by diameters"
  } else if(AVGLENGTHWSURF == 1) {
	Lbysurf_areabyasis = 1
  }
}

forall nseg=10
forall Ra=260


objectvar secaddress		
objectvar paraddress		
objectvar order			
objectvar parsecndx		
				
objectvar secri			
objectvar secpathri		
objectvar secL			
objectvar secRa			
objectvar secCm			
objectvar mrgL		
objectvar mrgdiam	
objectvar mrgri			
				
objectvar tosec			
objectvar tocyc			
				
objectvar rvp			
objectvar slsoma		
objectvar slroot		

objectvar cycL
objectvar cycdiam
objectvar cycnin





proc initsec(){
	secaddress=new Vector(NSEC,0)		//set/reset
	paraddress=new Vector(NSEC,0)		//set/reset
	order=new Vector(NSEC,0)		//set/reset
	parsecndx=new Vector(NSEC,0)		//set/reset

	secri=new Vector(NSEC,0)		//set/reset
	secpathri=new Vector(NSEC,0.0)		//set/reset
	secL=new Vector(NSEC,0)   		//set/reset
	secRa=new Vector(NSEC,0)   		//set/reset
	secCm=new Vector(NSEC,0)		//set/reset
}
proc initmrg(){
	mrgL=new Vector(NSEC,0)   		//set/reset
	mrgdiam=new Vector(NSEC,0)		//set/reset
	mrgri=new Vector(NSEC,0)		//set/reset

}
proc initcyc(){
	cycL=new Vector(11,0)
	cycdiam=new Vector(11,0)
	cycnin=new Vector(11,0)
}
proc initto(){
	tosec=new Vector(NSEC,-999)		//set/reset
	tocyc=new Vector(NSEC,-999)		//set/reset
}
proc initVectors(){
	initto()
	initsec()
	initmrg()
}
													
proc arborcharacterize(){local sec,sec1,sec2, aaa, x

	{NSEC=0  forall NSEC=NSEC+1}		
	initVectors()
	initcyc()

	sec=0
	forall {
		if(sec==0){
	          secaddress.set(0, this_section())
	          paraddress.set(0, parent_section(0))
	          secRa.set(0,Ra)
	          secL.set(0,L)

	          secCm.set(0,cm)
	          mrgL.set(0,L)
	          {for(x) if(x>0) aaa=aaa+area(x)   mrgdiam.set(0,aaa/PI/L)}
	          mrgri.set(0,secri.x(0))
		}
		if(sec!=0){
		  secaddress.set(sec, this_section())
		  paraddress.set(sec, parent_section(0))
		  {slsoma=new SectionList() rvp=new RangeVarPlot("v")
		   {soma  rvp.begin(.5)}  rvp.end(.5)  rvp.list(slsoma)}
		  {slroot=new SectionList() rvp=new RangeVarPlot("v")
		   {ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)}  ss=ss+1}}
		   rvp.end(.5)  rvp.list(slroot)}
	  	  forsec slroot order.set(sec, order.x(sec)+1)
		  for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x))
		  if(secri.x(sec)>9999) secri.set(sec, -9999)
		  forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x))
	}sec=sec+1}
	setsec()
	setmrg()
	for sec1=1,NSEC-1 {
		for sec2=0,NSEC-1 {
			if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)}
			}
		}
}
proc setsec(){
	sec=0
	forall {
		if(sec!=0){
		  secRa.set(sec,Ra)
		  secL.set(sec,L)
		  secCm.set(sec,cm)
	}sec=sec+1}
}
proc setmrg(){
	sec=0
	forall {
		if(sec!=0){
		  mrgL.set(sec,L)
		 
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
		  mrgri.set(sec,secri.x(sec))
	}sec=sec+1}
}

proc getnextABP(){local extent, sec
	extent=$1
	secA=0			
	for sec=1,NSEC-1 {if(tosec.x(sec)==-999){			
			    if(order.x(sec)>=order.x(secA)){	
			      if(abs(tocyc.x(sec))==extent){		
				secA=sec
 	}}}}
	secB=0
	for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){	
			    if(sec!=secA){
			      if(tosec.x(sec)==-999){			
				if(abs(tocyc.x(sec))==extent){		
			          secB=sec
	}}}}}
	secP=0
	if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA)
}	
proc shortnms(){
	{AmRa=0 AmL=0 Amdiam=0 Amri=0}
	{BmRa=0 BmL=0 Bmdiam=0 Bmri=0}
	{PmRa=0 PmL=0 Pmdiam=0 Pmri=0}
	if(secA>0) {
		AmRa=secRa.x(secA)
		AmL=mrgL.x(secA)  
		Amdiam=mrgdiam.x(secA) 
		Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
	}
	if(secB>0) {
		BmRa=secRa.x(secB)
		BmL=mrgL.x(secB)  
		Bmdiam=mrgdiam.x(secB) 
		Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
	}
	if(secP>0) {
		PmRa=secRa.x(secP)
		PmL=mrgL.x(secP)  
		Pmdiam=mrgdiam.x(secP) 
		Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
	}
	if(AREABYLONG==1){
	  if(secA>0 && secB>0){
	    if(AmL>BmL){
	      {BmL=AmL  mrgL.set(secB,BmL)}
	      {Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)}
	    }
	    if(BmL>AmL){
	      {AmL=BmL  mrgL.set(secA,AmL)}
	      {Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)}
	    }
	  }
	}
}










proc joinAP(){local newRa, newL, newri, newdiam
		newRa=(AmRa+PmRa)/2
		newL=AmL+PmL
		newri=Amri+Pmri
		newdiam=sqrt(newRa*newL*4/newri/PI/100)

		mrgL.set(secP, newL)
		mrgdiam.set(secP, newdiam)
		mrgri.set(secP, newri)
	tosec.set(secA, secP)
}
proc joinBP(){local newRa, newL, newri, newdiam
		newRa=(BmRa+PmRa)/2
		newL=BmL+PmL
		newri=Bmri+Pmri
		newdiam=sqrt(newRa*newL*4/newri/PI/100)

		mrgL.set(secP, newL)
		mrgdiam.set(secP, newdiam)
		mrgri.set(secP, newri)
	tosec.set(secB, secP)
}
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam

	AmL2BmL=AmL/BmL
	BmL2AmL=BmL/AmL
	if(AmL2BmL>=5) {
		joinAP()
		tosec.set(secB, -997)
	}
	if(BmL2AmL>=5) {
		joinBP()
		tosec.set(secA, -997)
	}
	if(AmL2BmL<5 && BmL2AmL<5>=CUTOFF9 && secpathri.x(sec)<1e10>=CUTOFF8 && secpathri.x(sec)<CUTOFF9>=CUTOFF7 && secpathri.x(sec)<CUTOFF8>=CUTOFF6 && secpathri.x(sec)<CUTOFF7>=CUTOFF5 && secpathri.x(sec)<CUTOFF6>=CUTOFF4 && secpathri.x(sec)<CUTOFF5>=CUTOFF3 && secpathri.x(sec)<CUTOFF4>=CUTOFF2 && secpathri.x(sec)<CUTOFF3>=CUTOFF1 && secpathri.x(sec)<CUTOFF2>=0 && secpathri.x(sec)<CUTOFF1>0) {	
			if (secB==0 && secP!=0) joinAP()
			if (secB!=0 && secP!=0) joinABP()
			if (secB==0 && secP==0) tagAcyc(extent)
			if (secB!=0 && secP==0) tagABcyc(extent)

			getnextABP(extent)
			shortnms()
		}
	}

	if(Lbylong_areabylong==1){
		//Lbylong_areabylong
	  for extent=1,10 {
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, maxL)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbysurf_areabyasis==1){
		//Lbysurf_areabyasis
	  for extent=1,10 {
		totsurf=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  totsurf=totsurf+surf
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbysurf_areabylong==1){
		//Lbysurf_areabylong
	  for extent=1,10 {
		totsurf=0
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  totsurf=totsurf+surf
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbydiam_areabyasis==1){
		//Lbydiam_areabyasis
	  for extent=1,10 {
		totdiam=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  totdiam=totdiam+mrgdiam.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbydiam_areabylong==1){
		//Lbydiam_areabylong
	  for extent=1,10 {
		maxL=0
		totdiam=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  totdiam=totdiam+mrgdiam.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbyavg_areabyasis==1){
		//Lbyavg__areabyasis
	  for extent=1,10 {
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbyavg_areabylong==1){
		//Lbyavg__areabylong
	  for extent=1,10 {
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}

}

arborcharacterize()

cycmake()

print "Section\tLength (um)\tDiameter (um)"

for ii=1,10 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)


mikey
Posts: 37
Joined: Mon Jun 19, 2006 3:06 pm

Post by mikey »

The error message returned is thus:

nrniv: syntax error in C:/address near line 651

I checked in the coding and this is line 651:

for ii=1,10 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)
csh
Posts: 52
Joined: Mon Feb 06, 2006 12:45 pm
Location: University College London
Contact:

Post by csh »

mikey wrote:The error message returned is thus:

nrniv: syntax error in C:/address near line 651

I checked in the coding and this is line 651:

for ii=1,10 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)
The error is upstream in line 358:

Code: Select all

  if(AmL2BmL<5 && BmL2AmL<5>=CUTOFF9 && secpathri.x(sec)<1e10>=CUTOFF8 && secpathri.x(sec)<CUTOFF9>=CUTOFF7 && secpathri.x(sec)<CUTOFF8>=CUTOFF6 && secpathri.x(sec)<CUTOFF7>=CUTOFF5 && secpathri.x(sec)<CUTOFF6>=CUTOFF4 && secpathri.x(sec)<CUTOFF5>=CUTOFF3 && secpathri.x(sec)<CUTOFF4>=CUTOFF2 && secpathri.x(sec)<CUTOFF3>=CUTOFF1 && secpathri.x(sec)<CUTOFF2>=0 && secpathri.x(sec)<CUTOFF1>0) {   
         if (secB==0 && secP!=0) joinAP()
         if (secB!=0 && secP!=0) joinABP()
         if (secB==0 && secP==0) tagAcyc(extent)
         if (secB!=0 && secP==0) tagABcyc(extent)

         getnextABP(extent)
         shortnms()
      } /* <-  */
   }
You close a block that has not been opened. I guess the hoc-parser will only complain about that once it has parsed the whole file trying to find matching opening and closing brackets, therefore you get the error message at the very end once it realizes there is a mismatch. You can quite easily pick up such errors using an editor that supports code collapsing - it will save you a lot of time.
When deleting line 647 (a call to a non-existing proc/func cycmake), your code will be read without error messages.
Christoph
mikey
Posts: 37
Joined: Mon Jun 19, 2006 3:06 pm

Post by mikey »

Crisoph, thanks so, so much for taking the time to look at this. I have made the 2 changes you suggested:

1) Take out the line:
cycmake()

2) Take out this bracket on line 358:

if (secB==0 && secP==0) tagAcyc(extent)
if (secB!=0 && secP==0) tagABcyc(extent)
getnextABP(extent)
shortnms()
} I TOOK THIS BRACKET OUT!!!!
}

One thing I should mention about this bracket is that I didn't pick it up as an "orphan" bracket. It did have a corresponding opening bracket according to my Jedit text editor. But I took out this closing bracket regardless.

Upon making these 2 changes - the code did compile without any error messages and returned output. BUT the output was no good - because it gave the length and diameter of all the 10 compartments as 0. Which is obviously no good. So, Im still somewhat stuck. I have put the code with these changes below.


Code: Select all


/* --------------------------------------------------------------------- 

   HOC CODE TO COLLAPSE A DENDRITIC TREE 
   ===================================== 

   Collapse a dendritic tree into 3 compartments "proximal", "middle" 
   and "distal".  

   The collapse is made such as the collapsed dendritic structure 
   preserves the axial resistance of the original structure. 
   The algorithm works by merging successive pairs of dendritic 
   branches into an equivalent branch (a branch that preserves the 
   axial resistance of the two original branches). 
   This merging of branches can be done according to different methods 
   selectable in the present code: 

   AREABYLONG 

   The code rejects the smallest original branch if it is more than 
   5 times smaller than the long orginal branch.  If AREABYLONG is 
   set to 1, then the small branch is rescaled to the longest one 
   (a new diameter is assigned such that the small branch conserves 
   its axial resistance).  This minimizes errors for dendritic 
   structures where branches have a large range of possible lengths. 
   AREABYLONG = 0 is the default. 


   AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF 

   During the merging of two branches, the length of the equivalent 
   branch is calculated by one out of three possible methods.  

   First, a simple average of the two branch length: 

      L = (L1+L2)/2 

   This is the method originally used by Bush and Sejnowski 
   (J. Neurosci. Methods 46: 159-166, 1993). 
   It is selected by setting AVGLENGTH=1. 

   Second, by weighting this average with the diameters of each branch 

      L = (diam1*L1 + diam2*L2)/(diam1+diam2) 

   This algorithm produces more fair merging in the case two branches of 
   very different diameters are to be merged.  It is selected by setting 
   AVGLENGTHWDIAM=1. 

   Third, by weighting the average with membrane area: 

      L = (area1*L1 + area2*L2)/(area1+area2) 

   This algorithm is better in the case two branches of very different 
   areas are to be merged.  It is selected by setting AVGLENGTHWSURF=1. 


   CUTOFF1, CUTOFF2 

   The dendritic branches are assigned to one of the three compartments 
   ("proximal", "middle" or "distal") based on their path axial 
   resistance (in Meg-Ohms) to the soma.  The cutoff values (CUTOFF1 
   and CUTOFF2) determine these cut-offs.  For example, CUTOFF1=10 means 
   that all branches laying within 10 Meg-Ohm of axial resistance from 
   the soma will be merged together into the same ("proximal") section. 
    

   Written by M. Neubig & A. Destexhe 
   Laval University, 1997; CNRS 2000 

-----------------------------------------------------------------------*/ 

CUTOFF1 = 7.3   // define the cutoff between "proximal" and "middle" 
CUTOFF2 = 10   // define the cutoff between "middle" and "distal" 
CUTOFF3 = 20 
CUTOFF4 = 30    
CUTOFF5 = 40    
CUTOFF6 = 50 
CUTOFF7 = 60    
CUTOFF8 = 65    
CUTOFF9 = 71 

AVGLENGTH = 0      // if 1, merge based on ageraged length 
AVGLENGTHWDIAM = 1   // if 1, merge based on avg length, weighted by diam 
AVGLENGTHWSURF = 0   // if 1, merge based on avg length, weighted by area 

AREABYLONG = 0      // if 1, rescale the length before merging 



xopen("tc-cell.geo")      // read geometry file 



Lbysurf_areabyasis = 0 
Lbysurf_areabylong = 0 
Lbydiam_areabyasis = 0 
Lbydiam_areabylong = 0 
Lbyavg_areabyasis = 0 
Lbyavg_areabylong = 0 

Lbylong_areabylong = 0 

if(AREABYLONG==1) { 
  print "Merge by rescaling length" 
  if(AVGLENGTH == 1) { 
   Lbyavg__areabylong = 1 
   print "Merging method based on averaged length" 
  } else if(AVGLENGTHWDIAM == 1) { 
   Lbydiam_areabylong = 1 
   print "Merging method based on averaged length, weigthed by diameters" 
  } else if(AVGLENGTHWSURF == 1) { 
   Lbysurf_areabylong = 1 
  } 
} else { 
  print "Do not rescale length" 
  if(AVGLENGTH == 1) { 
   Lbyavg__areabyasis = 1 
   print "Merging method based on averaged length" 
  } else if(AVGLENGTHWDIAM == 1) { 
   Lbydiam_areabyasis = 1 
   print "Merging method based on averaged length, weigthed by diameters" 
  } else if(AVGLENGTHWSURF == 1) { 
   Lbysurf_areabyasis = 1 
  } 
} 

forall nseg=10 
forall Ra=260 


objectvar secaddress       
objectvar paraddress       
objectvar order          
objectvar parsecndx       
             
objectvar secri          
objectvar secpathri       
objectvar secL          
objectvar secRa          
objectvar secCm          
objectvar mrgL       
objectvar mrgdiam    
objectvar mrgri          
             
objectvar tosec          
objectvar tocyc          
             
objectvar rvp          
objectvar slsoma       
objectvar slroot       

objectvar cycL 
objectvar cycdiam 
objectvar cycnin 





proc initsec(){ 
   secaddress=new Vector(NSEC,0)      //set/reset 
   paraddress=new Vector(NSEC,0)      //set/reset 
   order=new Vector(NSEC,0)      //set/reset 
   parsecndx=new Vector(NSEC,0)      //set/reset 

   secri=new Vector(NSEC,0)      //set/reset 
   secpathri=new Vector(NSEC,0.0)      //set/reset 
   secL=new Vector(NSEC,0)         //set/reset 
   secRa=new Vector(NSEC,0)         //set/reset 
   secCm=new Vector(NSEC,0)      //set/reset 
} 
proc initmrg(){ 
   mrgL=new Vector(NSEC,0)         //set/reset 
   mrgdiam=new Vector(NSEC,0)      //set/reset 
   mrgri=new Vector(NSEC,0)      //set/reset 

} 
proc initcyc(){ 
   cycL=new Vector(11,0) 
   cycdiam=new Vector(11,0) 
   cycnin=new Vector(11,0) 
} 
proc initto(){ 
   tosec=new Vector(NSEC,-999)      //set/reset 
   tocyc=new Vector(NSEC,-999)      //set/reset 
} 
proc initVectors(){ 
   initto() 
   initsec() 
   initmrg() 
} 
                                        
proc arborcharacterize(){local sec,sec1,sec2, aaa, x 

   {NSEC=0  forall NSEC=NSEC+1}       
   initVectors() 
   initcyc() 

   sec=0 
   forall { 
      if(sec==0){ 
             secaddress.set(0, this_section()) 
             paraddress.set(0, parent_section(0)) 
             secRa.set(0,Ra) 
             secL.set(0,L) 

             secCm.set(0,cm) 
             mrgL.set(0,L) 
             {for(x) if(x>0) aaa=aaa+area(x)   mrgdiam.set(0,aaa/PI/L)} 
             mrgri.set(0,secri.x(0)) 
      } 
      if(sec!=0){ 
        secaddress.set(sec, this_section()) 
        paraddress.set(sec, parent_section(0)) 
        {slsoma=new SectionList() rvp=new RangeVarPlot("v") 
         {soma  rvp.begin(.5)}  rvp.end(.5)  rvp.list(slsoma)} 
        {slroot=new SectionList() rvp=new RangeVarPlot("v") 
         {ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)}  ss=ss+1}} 
         rvp.end(.5)  rvp.list(slroot)} 
          forsec slroot order.set(sec, order.x(sec)+1) 
        for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x)) 
        if(secri.x(sec)>9999) secri.set(sec, -9999) 
        forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x)) 
   }sec=sec+1} 
   setsec() 
   setmrg() 
   for sec1=1,NSEC-1 { 
      for sec2=0,NSEC-1 { 
         if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)} 
         } 
      } 
} 
proc setsec(){ 
   sec=0 
   forall { 
      if(sec!=0){ 
        secRa.set(sec,Ra) 
        secL.set(sec,L) 
        secCm.set(sec,cm) 
   }sec=sec+1} 
} 
proc setmrg(){ 
   sec=0 
   forall { 
      if(sec!=0){ 
        mrgL.set(sec,L) 
       
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI)) 
        mrgri.set(sec,secri.x(sec)) 
   }sec=sec+1} 
} 

proc getnextABP(){local extent, sec 
   extent=$1 
   secA=0          
   for sec=1,NSEC-1 {if(tosec.x(sec)==-999){          
             if(order.x(sec)>=order.x(secA)){    
               if(abs(tocyc.x(sec))==extent){       
            secA=sec 
    }}}} 
   secB=0 
   for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){    
             if(sec!=secA){ 
               if(tosec.x(sec)==-999){          
            if(abs(tocyc.x(sec))==extent){       
                   secB=sec 
   }}}}} 
   secP=0 
   if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA) 
}    
proc shortnms(){ 
   {AmRa=0 AmL=0 Amdiam=0 Amri=0} 
   {BmRa=0 BmL=0 Bmdiam=0 Bmri=0} 
   {PmRa=0 PmL=0 Pmdiam=0 Pmri=0} 
   if(secA>0) { 
      AmRa=secRa.x(secA) 
      AmL=mrgL.x(secA)  
      Amdiam=mrgdiam.x(secA) 
      Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100) 
   } 
   if(secB>0) { 
      BmRa=secRa.x(secB) 
      BmL=mrgL.x(secB)  
      Bmdiam=mrgdiam.x(secB) 
      Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100) 
   } 
   if(secP>0) { 
      PmRa=secRa.x(secP) 
      PmL=mrgL.x(secP)  
      Pmdiam=mrgdiam.x(secP) 
      Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100) 
   } 
   if(AREABYLONG==1){ 
     if(secA>0 && secB>0){ 
       if(AmL>BmL){ 
         {BmL=AmL  mrgL.set(secB,BmL)} 
         {Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)} 
       } 
       if(BmL>AmL){ 
         {AmL=BmL  mrgL.set(secA,AmL)} 
         {Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)} 
       } 
     } 
   } 
} 










proc joinAP(){local newRa, newL, newri, newdiam 
      newRa=(AmRa+PmRa)/2 
      newL=AmL+PmL 
      newri=Amri+Pmri 
      newdiam=sqrt(newRa*newL*4/newri/PI/100) 

      mrgL.set(secP, newL) 
      mrgdiam.set(secP, newdiam) 
      mrgri.set(secP, newri) 
   tosec.set(secA, secP) 
} 
proc joinBP(){local newRa, newL, newri, newdiam 
      newRa=(BmRa+PmRa)/2 
      newL=BmL+PmL 
      newri=Bmri+Pmri 
      newdiam=sqrt(newRa*newL*4/newri/PI/100) 

      mrgL.set(secP, newL) 
      mrgdiam.set(secP, newdiam) 
      mrgri.set(secP, newri) 
   tosec.set(secB, secP) 
} 
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam 

   AmL2BmL=AmL/BmL 
   BmL2AmL=BmL/AmL 
   if(AmL2BmL>=5) { 
      joinAP() 
      tosec.set(secB, -997) 
   } 
   if(BmL2AmL>=5) { 
      joinBP() 
      tosec.set(secA, -997) 
   } 
   if(AmL2BmL<5 && BmL2AmL<5>=CUTOFF9 && secpathri.x(sec)<1e10>=CUTOFF8 && secpathri.x(sec)<CUTOFF9>=CUTOFF7 && secpathri.x(sec)<CUTOFF8>=CUTOFF6 && secpathri.x(sec)<CUTOFF7>=CUTOFF5 && secpathri.x(sec)<CUTOFF6>=CUTOFF4 && secpathri.x(sec)<CUTOFF5>=CUTOFF3 && secpathri.x(sec)<CUTOFF4>=CUTOFF2 && secpathri.x(sec)<CUTOFF3>=CUTOFF1 && secpathri.x(sec)<CUTOFF2>=0 && secpathri.x(sec)<CUTOFF1>0) {    
         if (secB==0 && secP!=0) joinAP() 
         if (secB!=0 && secP!=0) joinABP() 
         if (secB==0 && secP==0) tagAcyc(extent) 
         if (secB!=0 && secP==0) tagABcyc(extent) 

         getnextABP(extent) 
         shortnms()  
   } 

   if(Lbylong_areabylong==1){ 
      //Lbylong_areabylong 
     for extent=1,10 { 
      maxL=0 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec) 
           cycnin.set(extent, cycnin.x(extent)+1) 
         } 
      } 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100)) 
           cycdiam.set(extent, cycdiam.x(extent)+newdiam^2) 
         } 
      } 
      if(cycnin.x(extent)!=0) cycL.set(extent, maxL) 
      if(cycnin.x(extent)==0) cycL.set(extent, 0) 
      cycdiam.set(extent, sqrt(cycdiam.x(extent))) 
     } 
   } 
   if(Lbysurf_areabyasis==1){ 
      //Lbysurf_areabyasis 
     for extent=1,10 { 
      totsurf=0 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           surf=mrgL.x(sec)*PI*mrgdiam.x(sec) 
           totsurf=totsurf+surf 
           cycnin.set(extent, cycnin.x(extent)+1) 
         } 
      } 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           surf=mrgL.x(sec)*PI*mrgdiam.x(sec) 
           cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf) 
           cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2) 
         } 
      } 
      if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf) 
      if(cycnin.x(extent)==0) cycL.set(extent, 0) 
      cycdiam.set(extent, sqrt(cycdiam.x(extent))) 
     } 
   } 
   if(Lbysurf_areabylong==1){ 
      //Lbysurf_areabylong 
     for extent=1,10 { 
      totsurf=0 
      maxL=0 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec) 
           surf=mrgL.x(sec)*PI*mrgdiam.x(sec) 
           totsurf=totsurf+surf 
           cycnin.set(extent, cycnin.x(extent)+1) 
         } 
      } 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           surf=mrgL.x(sec)*PI*mrgdiam.x(sec) 
           cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf) 
           newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100)) 
           cycdiam.set(extent, cycdiam.x(extent)+newdiam^2) 
         } 
      } 
      if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf) 
      if(cycnin.x(extent)==0) cycL.set(extent, 0) 
      cycdiam.set(extent, sqrt(cycdiam.x(extent))) 
     } 
   } 
   if(Lbydiam_areabyasis==1){ 
      //Lbydiam_areabyasis 
     for extent=1,10 { 
      totdiam=0 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           totdiam=totdiam+mrgdiam.x(sec) 
           cycnin.set(extent, cycnin.x(extent)+1) 
         } 
      } 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec)) 
           cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2) 
         } 
      } 
      if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam) 
      if(cycnin.x(extent)==0) cycL.set(extent, 0) 
      cycdiam.set(extent, sqrt(cycdiam.x(extent))) 
     } 
   } 
   if(Lbydiam_areabylong==1){ 
      //Lbydiam_areabylong 
     for extent=1,10 { 
      maxL=0 
      totdiam=0 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec) 
           totdiam=totdiam+mrgdiam.x(sec) 
           cycnin.set(extent, cycnin.x(extent)+1) 
         } 
      } 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec)) 
           newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100)) 
           cycdiam.set(extent, cycdiam.x(extent)+newdiam^2) 
         } 
      } 
      if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam) 
      if(cycnin.x(extent)==0) cycL.set(extent, 0) 
      cycdiam.set(extent, sqrt(cycdiam.x(extent))) 
     } 
   } 
   if(Lbyavg_areabyasis==1){ 
      //Lbyavg__areabyasis 
     for extent=1,10 { 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           cycnin.set(extent, cycnin.x(extent)+1) 
         } 
      } 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           cycL.set(extent, cycL.x(extent)+mrgL.x(sec)) 
           cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2) 
         } 
      } 
      if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent)) 
      if(cycnin.x(extent)==0) cycL.set(extent, 0) 
      cycdiam.set(extent, sqrt(cycdiam.x(extent))) 
     } 
   } 
   if(Lbyavg_areabylong==1){ 
      //Lbyavg__areabylong 
     for extent=1,10 { 
      maxL=0 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec) 
           cycnin.set(extent, cycnin.x(extent)+1) 
         } 
      } 
      for sec=1,NSEC-1{ 
         if(tocyc.x(sec)==extent) { 
           cycL.set(extent, cycL.x(extent)+mrgL.x(sec)) 
           newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100)) 
           cycdiam.set(extent, cycdiam.x(extent)+newdiam^2) 
         } 
      } 
      if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent)) 
      if(cycnin.x(extent)==0) cycL.set(extent, 0) 
      cycdiam.set(extent, sqrt(cycdiam.x(extent))) 
     } 
   } 

} 

arborcharacterize() 


print "Section\tLength (um)\tDiameter (um)" 

for ii=1,10 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii) 

mikey
Posts: 37
Joined: Mon Jun 19, 2006 3:06 pm

Post by mikey »

Thanks so much to Dimitris who gave this tip:

"I have the impression that the problem is that you are trying to access memory which is not allocated, e.g. accessing beyond the bounds of an array. If you try a smaller number of compartments (e.g. 5 or 6) do you still get the problem?"

I have tried with 5 compartments - and I do still get this problem - not working.
mikey
Posts: 37
Joined: Mon Jun 19, 2006 3:06 pm

Post by mikey »

I am so, so sorry. But there has been an awful mistake. The coding got corrupted in the cutting and pasting to the forum and the code for the working 4 compartments and the code for the non-working 10 compartments was NOT cut and paste correctly. SO, I now will cut and paste, and check, the code again to see that it gets replicated accurately in the forum. I am so, so sorry. I dont really know how it has happened.

UPDATE - just tried posting again, and the posting corrupts the code again. Will have to contact the administrators about this. Will hopefully get this sorted over next hour or so.
csh
Posts: 52
Joined: Mon Feb 06, 2006 12:45 pm
Location: University College London
Contact:

Post by csh »

When I wrote that the code would be read without errors after deleting the call to cycmake(), I didn't mean this wouldn't alter the behaviour of your code. Actually, cycmake() seems to be the central part of it.
Unfortunately, I had some problems reading your modified code, especially this if-clause:

Code: Select all

   if(AmL2BmL<5 && BmL2AmL<5>=CUTOFF9 && secpathri.x(sec)<1e10>=CUTOFF8 && secpathri.x(sec)<CUTOFF9>=CUTOFF7 && secpathri.x(sec)<CUTOFF8>=CUTOFF6 && secpathri.x(sec)<CUTOFF7>=CUTOFF5 && secpathri.x(sec)<CUTOFF6>=CUTOFF4 && secpathri.x(sec)<CUTOFF5>=CUTOFF3 && secpathri.x(sec)<CUTOFF4>=CUTOFF2 && secpathri.x(sec)<CUTOFF3>=CUTOFF1 && secpathri.x(sec)<CUTOFF2>=0 && secpathri.x(sec)<CUTOFF1>0) {   
/*...*/
}
Did you disable html in your post? Enabling html sometimes leads to strange behaviour when posting code. Some errors might have originated from copying & pasting html-spoiled code.
Anyway, I think it would be a great idea to use the excellent to-do-list Mike Neubig posted earlier. To get you going, I have started refactoring the code simply by introducing a CUTOFF vector instead of distinct variables. For readability puposes, I will only post the parts of the code that I have modified. Line numbers refer to the original code.

Line 72:

Code: Select all

// Modified 09/27/06, CSH------------------------------------------------

n_comp=3 //number of compartments

objectvar CUTOFF
CUTOFF=new Vector(n_comp-1)
CUTOFF.x[0] = 7.3 // define the cutoff between "proximal" and "middle"
CUTOFF.x[1] = 71  // define the cutoff between "middle" and "distal"

// End modified----------------------------------------------------------
Line 415:

Code: Select all

// Modified 09/27/06, CSH------------------------------------------------

 for sec=1,NSEC-1 {
  if(secpathri.x(sec)>=CUTOFF.x[n_comp-2] && secpathri.x(sec)< 1e10) {
   tocyc.set(sec, -1)
  }
 }
 for (i=n_comp-3;i>=0;i-=1) {	
  for sec=1,NSEC-1 {
   if(secpathri.x(sec)>=CUTOFF.x[i] && secpathri.x(sec)< CUTOFF.x[i+1]) {
    tocyc.set(sec, -n_comp+i+1)
   }
  }
 }
 for sec=1,NSEC-1{
  if(secpathri.x(sec)>=0 && secpathri.x(sec)< CUTOFF.x[0]) {
   tocyc.set(sec, -n_comp)
  }
 }
// End modified ----------------------------------------------------------
From there, I think it would be a better idea to use a variable such as n_comp instead of explicit numbers. For instance, instead of using

Code: Select all

for extent=1,10
as you did you could use

Code: Select all

for extent=1,n_comp
That way, you will only need to change the number of compartments once in your code. This can save you hours of debugging.
Christoph
mikey
Posts: 37
Joined: Mon Jun 19, 2006 3:06 pm

Post by mikey »

Thanks so so much csh. Great!! What you have put looks great. What you have posted looks around about the same kind of thing that Michael Hines has just kindly sent me. What do you think? I think I'm right in saying your's and Michael's solutions are equivalent, but am not 100 percent. I have cut and paste the code that Michael Hines has sent me along with the explanatory email. The code should come out right now because I have turned HTML off now (excellent tip csh!!) - and from what I can see that solves the problem.

(NB. the Hines code below has been modified to yield 10 compartments. He originally sent me a version with for 3 compartments. Changed this line:
CUTOFF.append(0, 7.3, 71, 1e10)
to
CUTOFF.append(0, 7.3, 10, 20, 30, 40, 50, 60, 65, 71, 1e10)


---------

Enclosed is a collapse_n.hoc factorized as suggested by Mike Neubig
and based on the diff between the original collapse.hoc and your
putatively working 4 compartment variant. Collapse.hoc
and collapse_n.hoc produce identical results for 3 compartments.
I have no idea if the generalization to n compartments makes sense.
This is merely a code transformation that is identical to the
original for the CUTOFF.size == 4 case (three compartments)
where the first cutoff is 0 and the last cutoff is 1e10. If you wish to make additional compartments then just insert monotonically increasing numbers into

CUTOFF.append(0, 7.3, 71, 1e10)

You have my permission to add all my correspondence with regard to this
matter into the forum. In fact, please do so. Perhaps someone will
comment on the utility of the generalization. Also it is important
for the moderators to be aware of your experience with missing code
in the code fragments you posted.

------------

Code: Select all



/* ---------------------------------------------------------------------

	HOC CODE TO COLLAPSE A DENDRITIC TREE
	=====================================

   Collapse a dendritic tree into 3 compartments "proximal", "middle"
   and "distal".  

   The collapse is made such as the collapsed dendritic structure
   preserves the axial resistance of the original structure. 
   The algorithm works by merging successive pairs of dendritic 
   branches into an equivalent branch (a branch that preserves the 
   axial resistance of the two original branches).
   This merging of branches can be done according to different methods
   selectable in the present code:

   AREABYLONG

	The code rejects the smallest original branch if it is more than
	5 times smaller than the long orginal branch.  If AREABYLONG is
	set to 1, then the small branch is rescaled to the longest one 
	(a new diameter is assigned such that the small branch conserves
	its axial resistance).  This minimizes errors for dendritic 
	structures where branches have a large range of possible lengths.
	AREABYLONG = 0 is the default.


   AVGLENGTH, AVGLENGTHWDIAM, AVGLENGTHWSURF

	During the merging of two branches, the length of the equivalent 
	branch is calculated by one out of three possible methods.  

	First, a simple average of the two branch length:

		L = (L1+L2)/2

	This is the method originally used by Bush and Sejnowski
	(J. Neurosci. Methods 46: 159-166, 1993).
	It is selected by setting AVGLENGTH=1.

	Second, by weighting this average with the diameters of each branch

		L = (diam1*L1 + diam2*L2)/(diam1+diam2)

	This algorithm produces more fair merging in the case two branches of
	very different diameters are to be merged.  It is selected by setting
	AVGLENGTHWDIAM=1.

	Third, by weighting the average with membrane area:

		L = (area1*L1 + area2*L2)/(area1+area2)

	This algorithm is better in the case two branches of very different
	areas are to be merged.  It is selected by setting AVGLENGTHWSURF=1.


   CUTOFF1, CUTOFF2

	The dendritic branches are assigned to one of the three compartments
	("proximal", "middle" or "distal") based on their path axial 
	resistance (in Meg-Ohms) to the soma.  The cutoff values (CUTOFF1
	and CUTOFF2) determine these cut-offs.  For example, CUTOFF1=10 means
	that all branches laying within 10 Meg-Ohm of axial resistance from 
	the soma will be merged together into the same ("proximal") section.
	

   Written by M. Neubig & A. Destexhe
   Laval University, 1997; CNRS 2000

-----------------------------------------------------------------------*/

objref CUTOFF
CUTOFF = new Vector()
CUTOFF.append(0, 7.3, 10, 20, 30, 40, 50, 60, 65, 71, 1e10)
n=CUTOFF.size
//CUTOFF1 = 7.3	// define the cutoff between "proximal" and "middle"
//CUTOFF2 = 71	// define the cutoff between "middle" and "distal"

AVGLENGTH = 0		// if 1, merge based on ageraged length
AVGLENGTHWDIAM = 1	// if 1, merge based on avg length, weighted by diam
AVGLENGTHWSURF = 0	// if 1, merge based on avg length, weighted by area

AREABYLONG = 0		// if 1, rescale the length before merging



xopen("tc-cell.geo")		// read geometry file



Lbysurf_areabyasis = 0
Lbysurf_areabylong = 0
Lbydiam_areabyasis = 0
Lbydiam_areabylong = 0
Lbyavg_areabyasis = 0
Lbyavg_areabylong = 0

Lbylong_areabylong = 0

if(AREABYLONG==1) {
  print "Merge by rescaling length"
  if(AVGLENGTH == 1) {
	Lbyavg__areabylong = 1
	print "Merging method based on averaged length"
  } else if(AVGLENGTHWDIAM == 1) {
	Lbydiam_areabylong = 1
	print "Merging method based on averaged length, weigthed by diameters"
  } else if(AVGLENGTHWSURF == 1) {
	Lbysurf_areabylong = 1
  }
} else {
  print "Do not rescale length"
  if(AVGLENGTH == 1) {
	Lbyavg__areabyasis = 1
	print "Merging method based on averaged length"
  } else if(AVGLENGTHWDIAM == 1) {
	Lbydiam_areabyasis = 1
	print "Merging method based on averaged length, weigthed by diameters"
  } else if(AVGLENGTHWSURF == 1) {
	Lbysurf_areabyasis = 1
  }
}

forall nseg=10
forall Ra=260


objectvar secaddress		
objectvar paraddress		
objectvar order			
objectvar parsecndx		
				
objectvar secri			
objectvar secpathri		
objectvar secL			
objectvar secRa			
objectvar secCm			
objectvar mrgL		
objectvar mrgdiam	
objectvar mrgri			
				
objectvar tosec			
objectvar tocyc			
				
objectvar rvp			
objectvar slsoma		
objectvar slroot		

objectvar cycL
objectvar cycdiam
objectvar cycnin





proc initsec(){
	secaddress=new Vector(NSEC,0)		//set/reset
	paraddress=new Vector(NSEC,0)		//set/reset
	order=new Vector(NSEC,0)		//set/reset
	parsecndx=new Vector(NSEC,0)		//set/reset

	secri=new Vector(NSEC,0)		//set/reset
	secpathri=new Vector(NSEC,0.0)		//set/reset
	secL=new Vector(NSEC,0)   		//set/reset
	secRa=new Vector(NSEC,0)   		//set/reset
	secCm=new Vector(NSEC,0)		//set/reset
}
proc initmrg(){
	mrgL=new Vector(NSEC,0)   		//set/reset
	mrgdiam=new Vector(NSEC,0)		//set/reset
	mrgri=new Vector(NSEC,0)		//set/reset

}
proc initcyc(){
	cycL=new Vector(n,0)
	cycdiam=new Vector(n,0)
	cycnin=new Vector(n,0)
}
proc initto(){
	tosec=new Vector(NSEC,-999)		//set/reset
	tocyc=new Vector(NSEC,-999)		//set/reset
}
proc initVectors(){
	initto()
	initsec()
	initmrg()
}
													
proc arborcharacterize(){local sec,sec1,sec2, aaa, x

	{NSEC=0  forall NSEC=NSEC+1}		
	initVectors()
	initcyc()

	sec=0
	forall {
		if(sec==0){
	          secaddress.set(0, this_section())
	          paraddress.set(0, parent_section(0))
	          secRa.set(0,Ra)
	          secL.set(0,L)

	          secCm.set(0,cm)
	          mrgL.set(0,L)
	          {for(x) if(x>0) aaa=aaa+area(x)   mrgdiam.set(0,aaa/PI/L)}
	          mrgri.set(0,secri.x(0))
		}
		if(sec!=0){
		  secaddress.set(sec, this_section())
		  paraddress.set(sec, parent_section(0))
		  {slsoma=new SectionList() rvp=new RangeVarPlot("v")
		   {soma  rvp.begin(.5)}  rvp.end(.5)  rvp.list(slsoma)}
		  {slroot=new SectionList() rvp=new RangeVarPlot("v")
		   {ss=0 forsec slsoma{ {if(ss==1) rvp.begin(.5)}  ss=ss+1}}
		   rvp.end(.5)  rvp.list(slroot)}
	  	  forsec slroot order.set(sec, order.x(sec)+1)
		  for(x) if(x>0) secri.set(sec, secri.x(sec)+ri(x))
		  if(secri.x(sec)>9999) secri.set(sec, -9999)
		  forsec slroot for(x) if(x>0) secpathri.set(sec,secpathri.x(sec)+ri(x))
	}sec=sec+1}
	setsec()
	setmrg()
	for sec1=1,NSEC-1 {
		for sec2=0,NSEC-1 {
			if(secaddress.x(sec2)==paraddress.x(sec1)) {parsecndx.set(sec1,sec2)}
			}
		}
}
proc setsec(){
	sec=0
	forall {
		if(sec!=0){
		  secRa.set(sec,Ra)
		  secL.set(sec,L)
		  secCm.set(sec,cm)
	}sec=sec+1}
}
proc setmrg(){
	sec=0
	forall {
		if(sec!=0){
		  mrgL.set(sec,L)
		 
mrgdiam.set(sec,sqrt(secRa.x(sec)*secL.x(sec)*4/secri.x(sec)/100/PI))
		  mrgri.set(sec,secri.x(sec))
	}sec=sec+1}
}

proc getnextABP(){local extent, sec
	extent=$1
	secA=0			
	for sec=1,NSEC-1 {if(tosec.x(sec)==-999){			
			    if(order.x(sec)>=order.x(secA)){	
			      if(abs(tocyc.x(sec))==extent){		
				secA=sec
 	}}}}
	secB=0
	for sec=1,NSEC-1 {if(parsecndx.x(sec)==parsecndx.x(secA)){	
			    if(sec!=secA){
			      if(tosec.x(sec)==-999){			
				if(abs(tocyc.x(sec))==extent){		
			          secB=sec
	}}}}}
	secP=0
	if(abs(tocyc.x(parsecndx.x(secA)))==extent) secP=parsecndx.x(secA)
}	
proc shortnms(){
	{AmRa=0 AmL=0 Amdiam=0 Amri=0}
	{BmRa=0 BmL=0 Bmdiam=0 Bmri=0}
	{PmRa=0 PmL=0 Pmdiam=0 Pmri=0}
	if(secA>0) {
		AmRa=secRa.x(secA)
		AmL=mrgL.x(secA)  
		Amdiam=mrgdiam.x(secA) 
		Amri=(AmRa*AmL)/(PI*((Amdiam^2)/4)*100)
	}
	if(secB>0) {
		BmRa=secRa.x(secB)
		BmL=mrgL.x(secB)  
		Bmdiam=mrgdiam.x(secB) 
		Bmri=(BmRa*BmL)/(PI*((Bmdiam^2)/4)*100)
	}
	if(secP>0) {
		PmRa=secRa.x(secP)
		PmL=mrgL.x(secP)  
		Pmdiam=mrgdiam.x(secP) 
		Pmri=(PmRa*PmL)/(PI*((Pmdiam^2)/4)*100)
	}
	if(AREABYLONG==1){
	  if(secA>0 && secB>0){
	    if(AmL>BmL){
	      {BmL=AmL  mrgL.set(secB,BmL)}
	      {Bmdiam=sqrt((BmRa*BmL*4)/(Bmri*PI*100)) mrgdiam.set(secB,Bmdiam)}
	    }
	    if(BmL>AmL){
	      {AmL=BmL  mrgL.set(secA,AmL)}
	      {Amdiam=sqrt((AmRa*AmL*4)/(Amri*PI*100)) mrgdiam.set(secA,Amdiam)}
	    }
	  }
	}
}










proc joinAP(){local newRa, newL, newri, newdiam
		newRa=(AmRa+PmRa)/2
		newL=AmL+PmL
		newri=Amri+Pmri
		newdiam=sqrt(newRa*newL*4/newri/PI/100)

		mrgL.set(secP, newL)
		mrgdiam.set(secP, newdiam)
		mrgri.set(secP, newri)
	tosec.set(secA, secP)
}
proc joinBP(){local newRa, newL, newri, newdiam
		newRa=(BmRa+PmRa)/2
		newL=BmL+PmL
		newri=Bmri+Pmri
		newdiam=sqrt(newRa*newL*4/newri/PI/100)

		mrgL.set(secP, newL)
		mrgdiam.set(secP, newdiam)
		mrgri.set(secP, newri)
	tosec.set(secB, secP)
}
proc joinABP(){local AmBmL,AmBmri,AmBmdiam, newRa,newL,newdiam

	AmL2BmL=AmL/BmL
	BmL2AmL=BmL/AmL
	if(AmL2BmL>=5) {
		joinAP()
		tosec.set(secB, -997)
	}
	if(BmL2AmL>=5) {
		joinBP()
		tosec.set(secA, -997)
	}
	if(AmL2BmL<5 && BmL2AmL<5){
		if(Lbysurf_areabyasis==1 || Lbysurf_areabylong==1) {
			//Lbysurf_areabyasis BSC_Lbysurf_areabylong
			Asurf=AmL*PI*Amdiam
			Bsurf=BmL*PI*Bmdiam
			AmBmL=(AmL*Asurf+BmL*Bsurf)/(Asurf+Bsurf)
			AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)

			AmBmcross=PI*(AmBmdiam^2)/4	
			AmBmRa=(AmRa+BmRa)/2
			AmBmri=AmBmRa*AmBmL/AmBmcross/100	

			newRa=(AmBmRa+PmRa)/2
			newL=AmBmL+PmL
			newri=AmBmri+Pmri
			newdiam=sqrt(newRa*newL*4/newri/PI/100)
		}
		if(Lbydiam_areabyasis==1 || Lbydiam_areabylong==1){
			//Lbydiam_areabyasis  Lbydiam_areabylong
			AmBmL=(AmL*Amdiam+BmL*Bmdiam)/(Amdiam+Bmdiam)
			AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
			AmBmcross=PI*(AmBmdiam^2)/4	

			AmBmRa=(AmRa+BmRa)/2
			AmBmri=AmBmRa*AmBmL/AmBmcross/100	

			newRa=(AmBmRa+PmRa)/2
			newL=AmBmL+PmL
			newri=AmBmri+Pmri
			newdiam=sqrt(newRa*newL*4/newri/PI/100)
		}
		if(Lbyavg__areabyasis==1 || Lbyavg__areabylong==1 || Lbylong_areabylong==1) {
			//Lbyavg__areabyasis Lbyavg__areabylong Lbylong_areabylong
			AmBmL=(AmL+BmL)/2
			AmBmdiam=sqrt(Amdiam^2+Bmdiam^2)
			AmBmcross=PI*(AmBmdiam^2)/4	

			AmBmRa=(AmRa+BmRa)/2
			AmBmri=AmBmRa*AmBmL/AmBmcross/100	

			newRa=(AmBmRa+PmRa)/2
			newL=AmBmL+PmL
			newri=AmBmri+Pmri
			newdiam=sqrt(newRa*newL*4/newri/PI/100)
		}
		mrgL.set(secP, newL)
		mrgri.set(secP, newri)
		mrgdiam.set(secP, newdiam)

		tosec.set(secA, secP)
		tosec.set(secB, secP)
	}
	
}
proc tagAcyc(){
	tosec.set(secA, -$1)
	tocyc.set(secA, $1)
}
proc tagBcyc(){
	tosec.set(secB, -$1)
	tocyc.set(secB, $1)
}
proc tagABcyc(){
	
	tagAcyc($1)
	tagBcyc($1)
}
proc cycmake(){local ic, sec, extent
	initto()
	initcyc()
	setmrg()

	for ic = 2, n {
		for sec = 1, NSEC-1 {
			if (secpathri.x(sec) >= CUTOFF.x[n-ic] \
			    && secpathri.x(sec) < CUTOFF.x[n - ic + 1]) {
				tocyc.set(sec, 1 - ic)
			}
		}
	}
	for extent=1,n-1 {
		getnextABP(extent)
		shortnms()
		while(secA>0) {	
			if (secB==0 && secP!=0) joinAP()
			if (secB!=0 && secP!=0) joinABP()
			if (secB==0 && secP==0) tagAcyc(extent)
			if (secB!=0 && secP==0) tagABcyc(extent)

			getnextABP(extent)
			shortnms()
		}
	}

	if(Lbylong_areabylong==1){
		//Lbylong_areabylong
	  for extent=1,n-1 {
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, maxL)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbysurf_areabyasis==1){
		//Lbysurf_areabyasis
	  for extent=1,n-1 {
		totsurf=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  totsurf=totsurf+surf
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbysurf_areabylong==1){
		//Lbysurf_areabylong
	  for extent=1,n-1 {
		totsurf=0
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  totsurf=totsurf+surf
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*surf)
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totsurf)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbydiam_areabyasis==1){
		//Lbydiam_areabyasis
	  for extent=1,n-1 {
		totdiam=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  totdiam=totdiam+mrgdiam.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbydiam_areabylong==1){
		//Lbydiam_areabylong
	  for extent=1,n-1 {
		maxL=0
		totdiam=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  totdiam=totdiam+mrgdiam.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec)*mrgdiam.x(sec))
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/totdiam)
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbyavg_areabyasis==1){
		//Lbyavg__areabyasis
	  for extent=1,n-1 {
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
			  cycdiam.set(extent, cycdiam.x(extent)+mrgdiam.x(sec)^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}
	if(Lbyavg_areabylong==1){
		//Lbyavg__areabylong
	  for extent=1,n-1 {
		maxL=0
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
			  cycnin.set(extent, cycnin.x(extent)+1)
			}
		}
		for sec=1,NSEC-1{
			if(tocyc.x(sec)==extent) {
			  cycL.set(extent, cycL.x(extent)+mrgL.x(sec))
			  newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
			  cycdiam.set(extent, cycdiam.x(extent)+newdiam^2)
			}
		}
		if(cycnin.x(extent)!=0) cycL.set(extent, cycL.x(extent)/cycnin.x(extent))
		if(cycnin.x(extent)==0) cycL.set(extent, 0)
		cycdiam.set(extent, sqrt(cycdiam.x(extent)))
	  }
	}

}

arborcharacterize()

cycmake()

print "Section\tLength (um)\tDiameter (um)"

for ii=1,n-1 print ii, "\t", cycL.x(ii),"\t", cycdiam.x(ii)

Post Reply