Prune.C

#include"Prune.h"

/*****************************************************************/

void Prune :: initialize()
{
	int i;
	isRead();
	ShollAnalysis();

	npc=ncomp;
	PruCom=new SWCData[npc];
	for(i=0; i<npc; i++)
		PruCom[i]=Comptmt[i];


	delete(Comptmt);

	getSoma();
	getStems();
	getBP();
	getTP();



	Children=new Subset[npc] ;
	for(i=0; i<npc; i++){
		Children[i].npts=0;
		Children[i].Pts=new int [3];	// Maximum is trifurcation
	}	

	findChildren();

 	//generateSpecs(); exit(1);

	abpprune=new double[sdata.na]; bbpprune=new double[sdata.nb];
	adlprune=new double[sdata.na]; bdlprune=new double[sdata.nb];
	abpP    =new double[sdata.na]; bbpP    =new double[sdata.nb];
	adlP    =new double[sdata.na]; bdlP    =new double[sdata.nb];

	for(i=0; i<sdata.na; i++)		// Amount of initial pruning is zero
		adlprune[i]=abpprune[i]=0;
	for(i=0; i<sdata.nb; i++)
		bdlprune[i]=bbpprune[i]=0;

	cmr=0;	// Reduction in number of compartments
	dlr=0;	// Reduction in dl.

	setZeta(basefilename);
	//setBasalMaxZeta(basefilename);
	setProbabilities();

/* findTPBP requires probabilities to be initialized and hence is called
 * after setProbabilities() */

	TPBP=new int[TP.npts];
	bpPruneP=new double[TP.npts];
	findTPBP();
}

/*****************************************************************/

void Prune :: PruneTree()
{
	int i;
	initialize();
	prune();
//	cerr << "\nNo. of Soma Pts: " << Soma.npts ;
//	cerr << "\nNo. of Stems: " << Stems.npts ;
//	cerr << "\nNo. of BP's: " << BP.npts ;
//	cerr << "\nNo. of TP's: " << TP.npts << endl ;
	char  * filename=new char[35];
	sprintf(filename,"%s_final.swc",outfilename);
	cerr << "\nFinal save in " << filename << " after pruning " 
		<< dlr << " micron of length" << endl ;
	saveSWC(filename);
	
	//printStatus();
}

/*****************************************************************/

void Prune :: prune()
{
	int i;
	struct timeval tv;
	struct timezone tz;
	gettimeofday(&tv,&tz);
	int SEED=tv.tv_usec;
	srand48(SEED);

	int count =0;
	int scount = 1000; // Save for every 1000 um reduction in DL.
	char  * filename=new char[35];

	//while(convmeasure()>0.05){
	while(dlr<dlred){
	//while(dlr<dlsum){
	//while (count < 5000000){
		if(TP.npts==0){
			cerr << "\nNo more to Prune...\n";
			break;
		}	
		for(i=0; i<TP.npts; i++){
			bpPrune(i);
			tpPrune(i);
		}	
		//cerr << count << "\t" << cmr << "\t" << dlr << "\t" 
		//	 << convmeasure() << endl ;
		count ++;
		if(dlr>scount){
			sprintf(filename,"%s_%d.swc",outfilename,scount);
			cerr << "\nSaving file " << filename << " after pruning " 
				<< dlr << " micron of length" << endl ;
			saveSWC(filename);
			scount += 1000;
		}
	}	

	//cerr << "\nReduction in DL: " << dlr ;
	//cerr << "\nReduction in Cmptmt: " << cmr << endl ;
}

/*****************************************************************/
// Prunes TP[i]

void Prune :: tpPrune(int i)
{
	double dist;
	int l;
	int sn;	// Segment number
	int cn;		// Compartment number
	int PRUNE;

	/* BP has been pruned and the number of TP's has been reduced.
	 * So, the value of i will be greater than what is allowed. We
	 * just return back without bothering much as it does not mean
	 * anything */ 

	if(i>(TP.npts-1)){
		return;
	}	

	cn    = TP.Pts[i];
	dist  = 0.0;
	dist += (somax-PruCom[cn].x)*(somax-PruCom[cn].x);
	dist += (somay-PruCom[cn].y)*(somay-PruCom[cn].y);
	dist += (somaz-PruCom[cn].z)*(somaz-PruCom[cn].z);
	dist  = sqrt(dist);
	sn   = (int)(dist/(double)STLENGTH); // Segment where TP[i] belongs
	PRUNE=0;
	if (PruCom[cn].tc==4){	// Apical
		if(drand48()*maxdlP<adlP[sn]) PRUNE=1;
	}
	else if (PruCom[cn].tc==3){	// Basal
		if(drand48()*maxdlP<bdlP[sn]) PRUNE=1;
	}
	else{
		cerr << "\ntpP :: Terminal point is not dendrite!!!\n\n";
		checkTPs();
		cerr <<  i << endl << PruCom[cn] ;
		cerr <<  Children[cn].npts << endl  ;
		for(l=0; l<Children[cn].npts; l++){
			cerr << PruCom[Children[cn].Pts[l]];
		}
		exit(1);
	}	

	if(PRUNE){
		pruneTP(cn,sn);
		//setProbabilities();		// Update probability
		findChildren();			// Update children list
		findTPBP();				// Update TPBP list
	}	
}	

/*****************************************************************/
// Prunes TPBP[i]

void Prune :: bpPrune(int i)
{
	double dist;
	int sn;	// Segment number
	int cn;		// Compartment number
	int PRUNE;
	cn    = TPBP[i];
	
	if(cn==-2) return; // Case where the TP does not end on a BP.
	
	dist  = 0.0;
	dist += (somax-PruCom[cn].x)*(somax-PruCom[cn].x);
	dist += (somay-PruCom[cn].y)*(somay-PruCom[cn].y);
	dist += (somaz-PruCom[cn].z)*(somaz-PruCom[cn].z);
	dist  = sqrt(dist);
	sn   = (int)(dist/(double)STLENGTH); // Segment where TP[i] belongs
	PRUNE=0;
	if (PruCom[cn].tc==4){	// Apical
		//if(drand48()*maxdlP<(bpPruneP[i]*abpP[sn])) PRUNE=1;
		if(drand48()*maxdlP<bpPruneP[i]) PRUNE=1;
	}
	else if (PruCom[cn].tc==3){	// Basal
		//if(drand48()*maxdlP<(bpPruneP[i]*bbpP[sn])) PRUNE=1;
		if(drand48()*maxdlP<bpPruneP[i]) PRUNE=1;
	}
	else{
		cerr << "\nBranching point is not dendrite!!!\n\n";
		//cerr << j << " " << i << " "<< PruCom[cn] ;
		exit(1);
	}	
	if(PRUNE){
		pruneBP(i);
		//setProbabilities();		// Update probability
		findChildren();
		findTPBP();
	}	
}

/*****************************************************************/
// Convergence measure.

double Prune :: convmeasure()
{
	double cm=0.0;
	int i;
	
/* If the percentage reduction is negative it is not considered. Can be
 * reviewed at a later stage. TEST */

	for(i=0; i<sdata.na; i++){
		if((adlzeta[i]-adlprune[i])/adlzeta[i]>cm)
			cm=fabs((adlzeta[i]-adlprune[i])/adlzeta[i]);
		if((abpzeta[i]-abpprune[i])/abpzeta[i]>cm)
			cm=fabs((abpzeta[i]-abpprune[i])/abpzeta[i]);
	}		

	for(i=0; i<sdata.nb; i++){
		if((bdlzeta[i]-bdlprune[i])/bdlzeta[i]>cm)
			cm=fabs((bdlzeta[i]-bdlprune[i])/bdlzeta[i]);
		if((bbpzeta[i]-bbpprune[i])/bbpzeta[i]>cm)
			cm=fabs((bbpzeta[i]-bbpprune[i])/bbpzeta[i]);
	}		

	return cm;
}

/*****************************************************************/

void Prune :: saveSWC(char * filename)
{
	char fname[35];
	if(!filename) strcpy (fname,outfilename);
	else strcpy(fname,filename);

   	ofstream outfile (fname);
   	for(int i=0; i<npc; i++)
       	outfile << PruCom[i] ;
   	outfile.close();
}

/*****************************************************************/

void Prune :: checkTPs()
{
	int cn;
	for (int i=0; i<TP.npts; i++){
		cn=TP.Pts[i];
		if((PruCom[cn].tc != 3) && (PruCom[cn].tc != 4)){
    		cerr << "\ncheckTPs :: TP not dendrite\n";
			exit(1);
		}
	}
}

/*****************************************************************/