#include "SphierStruct.h"
#include <cmath>

using namespace std;

typedef pair<Node*,int> Neighbor;

/*Sphier struct*/
SphierStruct::SphierStruct():
   staticLayer(0),
   ID0(0),
   ID1(0),
   node0Array(),
   node1Array(),
   improvement(false)
   {
    
   }

/*Another constructir for SphierStruct*/
SphierStruct::SphierStruct(int rows, int cols, double threshold,int rowt, int colt, char* fname):
	staticLayer(0),
	ID0(0),
	ID1(0),
	rows(rows),
	cols(cols),
	rowt(rowt),
	colt(colt),
	threshold(threshold),
	node0Array(rows,0),
    node1Array(cols,0),
    fname(fname),
    improvement(false)
	{
	}
/*Destructor*/
SphierStruct::~SphierStruct()
{
//	delete node0Array;
//	delete node1Array;
//	delete table;
}

/*Add a node*/
void SphierStruct::addNode(Node& node)
{
       //.insert(pair<int,int>(node.ID,node.position));
	if(node.layerNum == 0){
		node0Array.add(node);
		//ID0 = ID0+1;
	}
	if(node.layerNum == 1){
		node1Array.add(node);
		//ID1 = ID1+1;
	}
}
Node* SphierStruct::findNode(Node* node)
{
	if(node->layerNum == 0)
		return &node0Array[node->ID];
	if(node->layerNum == 1)
		return &node1Array[node->ID];
	
}

Node* SphierStruct::findNode(int ID, int layer)
{
	Node* result;
//	Node* node = new Node(ID, layer);
//        Node* temp = findNode(node);
//        if(temp->ID != ID) cout<<"Warning: ID of the extracted node dont match with the given ID"<<endl;
    if(layer == 0)
	    result = &node0Array[ID];
	else
	    result = &node1Array[ID];
	return result;
}

void SphierStruct::addNode(std::string name, int Layer){
	
	//cout<<"Inside node add"<<endl;
	//Node temp(ID, Layer, name);
	if (Layer == 0){
		Node temp(ID0, Layer, name);
		
		//addNode(temp);
		 node0Array[ID0] = temp;
		 //cout<<"node added in array"<<endl;
		 table[0][ID0]= temp.position;
		 ID0 = ID0+1;
		//cout<<"added to table"<<endl;
	}
	if(Layer == 1){
		Node temp(ID1, Layer, name);
		//addNode(temp);
		node1Array[ID1] = temp;
		table[1][ID1]= temp.position;
		ID1 = ID1+1;
		
	}
	
}
/*set the static layer*/
int SphierStruct::invertStaticLayer()
{
	int result = staticLayer;
    staticLayer = 1-staticLayer;
	return result;
}
void SphierStruct::setPositionBary(int layer)
{
    improvement = false;
	if(layer == 0){
		int length = node0Array.Length();
		for(int i=0; i< length; i++){
			node0Array[i].setPositionBary();
			//cerr<<"node has position changed"<<endl;
            if(node0Array[i].positionChanged == true) improvement = true;
            else improvement = false;
                //cerr<<"SS::SPB-> Node "<<node0Array[i].name<<" has now position ="<<node0Array[i].position<<endl;
			table[0][node0Array[i].ID]= node0Array[i].position;
                //cerr<<"SS::SPB-> Table Node0"<<node0Array[i].ID<<" has now position ="<<node0Array[i].position<<endl;
		}
	}
	if(layer == 1){
		int length = node1Array.Length();
		for(int i=0; i< length; i++){
			node1Array[i].setPositionBary();
            if(node1Array[i].positionChanged == true) improvement = true;
            else improvement = false;
            	//cerr<<"SS::SPB-> Node "<<node1Array[i].name<<" has now position ="<<node1Array[i].position<<endl;
			table[1][node1Array[i].ID]= node1Array[i].position;
		}
	}
        adjustNodePositions(layer);
}

void SphierStruct::adjustNodePositions(int layer)
{
    multimap <int,int> rev;
    map <int,int>::iterator iter;
    Array <Node>* ptr;
    if(layer == 0) ptr = &node0Array;
    if(layer == 1) ptr = &node1Array;
    for(iter = table[layer].begin(); iter != table[layer].end(); iter++)
    {
        //cout << (*iter).first << " is the first character in the word " << iter->second << endl;    
        //rev[iter->second] = iter->first;
          rev.insert(pair<int,int>(iter->second,iter->first));
    }
    int index = 0;
    multimap <int,int>::iterator iter2;
    for(iter2=rev.begin(); iter2!= rev.end(); iter2++){
       Node* temp = &ptr->elementAt(iter2->second);
       temp->setPosition(index);
       table[layer][temp->ID]= temp->position;
      // cerr<<"SS::ANP-> Node "<<temp->name<<" with ID="<<temp->ID<<" has adjusted position ="<<index<<endl;
       index++;
	   //cerr<<"adjust completed"<<endl;
    }
    //test
//    for(int i=0; i < node0Array.Length(); i++){
//    	cerr<<"Row "<<node0Array[i].name<<" Position= "<<node0Array[i].position<<endl;
//    }
//    for(int i=0; i < node1Array.Length(); i++){
//    	cerr<<"Column "<<node1Array[i].name<<" Position= "<<node1Array[i].position<<endl;
//    }
}

/*Barycenter heuristic based crossing minimization*/
void SphierStruct::setPositionBary()
{
    setPositionBary(invertStaticLayer());	
}

/* print the structure*/
void SphierStruct::printStruct()
{
	
	for(int layer =0; layer < 2; layer++){
		int length;
		if(layer == 0) length=node0Array.Length();
		else length = node1Array.Length();
		cout<<" Length = "<<length;
    	cout<<" Layer = "<<layer;
		cout<<endl;
		for(int i=0; i < length; i++){
			Node* ptr = findNode(i,layer);
			cout<<"Node: = "<<ptr->name;
			cout <<"   Neighbors: ";
			LinkedList<Neighbor> neighbors = ptr->neighbors;
			const ListElement<Neighbor>* finger = neighbors.Head();
			for(int j=0; j<neighbors.size; j++){
				Neighbor nb = finger->Datum();
				Node* n = nb.first;
				cout<<n->name<<" "; 
				finger = finger->Next();
			}
			cout<<endl;
		}
		cout<<endl;
	}

}

double SphierStruct::evaluate(double num, int algo){
	double retval;
	if (algo == 0){ // KL Divergence
//		if(num ==0 ) retval = 1;
//		else
			retval = (num)*log(abs(num));
	}
	if(algo == 1){ // Euclidean Distance
		retval = num*num;
	}
	if(algo == 2){ //Logistic Loss
//		if(num == 0 || (1-num)==0)
//			retval = 1;
//		else
			retval = num*log(abs(num))+(1-num)*log(abs(1-num));
	}
	if(algo == 3){ // Itakura-Saito Distance
//		if (num == 0) retval =1;
//		else
			retval = -1*log(abs(num));
	}
	if(algo ==4){ // Hinge Loss
		retval = abs(num);
	}
	if(algo == 5){// mahalonobis distance
		retval = num*num;
	}
	return retval;
}

double SphierStruct::evalGrad(double fq, int algo){
	double result;
	if(algo == 0){// KL divergenc
//		if(fq == 0) return 1;
		result = 1+log(abs(fq));
	}
	if(algo == 1){//Euclidean Distance
//		if(fq == 0) return 1;
//		else
			result = fq/abs(fq)*2*fq;
	}
	if(algo == 2){// Logistic Loss
//		if(fq ==0 || (1-fq) == 0) return 1;
		result = 1+log(fq)-log(1-fq)+1;
	}
	if(algo == 3){// Itakura-Saito Distance
//		if(fq == 0) return 1;
		result = -1/fq;
	}
	if(algo == 4){// HInge loss
		result = 1;
	}
	if(algo == 5){// Mahalonobis Distance
		result = 2*fq;
	}
	return result;
}

double SphierStruct::calculateBD(double fp, double fq, double p, 
									double q, double fqd){
	double retval;
	retval = fp-fq-(fqd*(p-q));
	return retval;
}


void SphierStruct::findBiclusters(int algo)
{
//	for(int i=0; i < node0Array.Length(); i++){
//    	cerr<<"Row "<<node0Array[i].name<<" Position= "<<node0Array[i].position<<endl;
//    }
//    for(int i=0; i < node1Array.Length(); i++){
//    	cerr<<"Column "<<node1Array[i].name<<" Position= "<<node1Array[i].position<<endl;
//    }
	//delete [] table;
	int rows = node0Array.Length();
    int cols = node1Array.Length();
    //cerr<<"rows "<<rows<<" cols "<<cols<<endl;
    int matrix[rows][cols];
    Node* matrix2[rows][cols];
    Node* matrix3[rows];
   // cerr<<"going into loop"<<endl;
    for(int i=0; i<rows ;i++){
    	//if(&node0Array[i] != 0){
	    	//int x = table[0][i];
	    	int x = node0Array[i].position;
    		matrix3[x] = &node0Array[i];
    		//cerr<<"matrix3[x] for x";
    		//cerr<<" = "<< matrix3[x]->name<<endl;
        	for(int j=0; j<cols; j++){
            	//int y = table[1][j];
            	int y = node1Array[j].position;
            	//cerr<<"SSP::PM-> x ="<<x<<"  and y="<<y<<endl;
            	pair<Node*, int> res = node0Array[i].findNeighbor(&node1Array[j]);
            	if(res.first != NULL){
            		//cerr<<" Neighbor was "<<endl;
            		//cerr<<res.first->name<<endl;
                	matrix[x][y] = res.second;
                	matrix2[x][y] = res.first;
                	//cerr<<"leaving if"<<endl;
                	//cerr<<"matrix2[x][y] with x = "<< x <<"and y= "<< y <<"="<< matrix2[x][y]->name<<endl;
                   
            	}
           		else{
           			//cerr<<"SSP:: in else"<<endl;
                	matrix[x][y]=0;
                	matrix2[x][y]=0;
                	//cerr <<"out of else"<<endl;
            	}
        	}
    	//} 
    }
    
	LinkedList<string> rowcluster;
	LinkedList<string> prevrc;
	LinkedList<string> colcluster;
	LinkedList<string> prevcc;
	typedef pair<LinkedList<string>, LinkedList<string> > Bicluster;
	map<int, Bicluster> biclusters;
	map<int,double> scores;
	
	//cerr<<"now finding bixes"<<endl;
	ofstream mfile;
    mfile.open("output-matrix.cHawk");
    mfile <<"Threshold ="<<threshold<<endl;
    mfile<<"colt"<<colt<<" rowt ="<<rowt<<endl;
	
    	for(int i=0; i < rows; i++){
    		for(int j =0; j< cols; j++){
    			mfile<<matrix[i][j]<<" ";
    		}
    		mfile<<endl;
    	}
    
	int colcount=0;
	int rowcount=0;
	int startrow =0;
	int startcol =0;
	int index =0;
	int prevrow=startrow;
	int prevcol =startcol;
	int stopcol;
	double fp=0;
	double fq=0;
	double p=0;
	double q=0;
	double dist=0;
	double blockmean;
	int ovlprow;
	int ovlpcol;
	bool set=false;
	int coldist=0;
	int coldist2=0;
	double prevp,prevq,prevfp, prevfq,prevpd,prevqd;
	for(int i=startrow; i < rows-1; i++){
		for(int j=startcol; j< cols;j++){
			//if(matrix[i][j] == 0) continue;
			//cerr<<"i ="<<i<<" j ="<<j<<endl;
			int den = (i-startrow)*(j-startcol)+1;
			fp = evaluate(matrix[i][j],algo);
			
			//fp = fp/den;
			fq = evaluate(matrix[i+1][j],algo);
			//fq = fq/den;
			
			
			if(j != startcol){
				prevfp += evaluate(matrix[i][j],algo);
				prevfq += evaluate(matrix[i+1][j],algo);
				//mfile<<"j > startcol"<<endl;
			}
			else{
				prevfp = fp; 
				prevfq = fq;
			}
			
			p = (double)matrix[i][j];
			//p = p/den;
			q = (double)matrix[i+1][j];
			double newmean = p+q/2;
			//if(abs(blockmean-newmean) > newmean ) continue;
			blockmean += newmean;
			//double prevp, prevq;
			if(j != startcol){
				prevp = abs(prevp-(matrix[i][j]));
				prevq = abs(prevq-(matrix[i+1][j]));
				blockmean=(p+q)/2;
				//mfile<<"j > startcol"<<endl;
			}
			else{
				prevp = p; 
				prevq = q;
			}
			//q = q/den;
			double fqd= evalGrad(q,0);
			prevpd = evalGrad(prevp,0);
			prevqd = evalGrad(prevq,0);
			//mfile<<"p ="<<p<<" q= "<< q <<" fqd="<<fqd<<endl;
			//mfile<<"prevp ="<<prevp<<" prevq= "<< prevq <<" prevpd="<<prevpd<<endl;
			//mfile<<"fp "<<fp<<" fq="<<fq<<endl;
			dist = calculateBD(fp,fq,p,q,fqd);
			double distp = calculateBD(prevfp,fp,prevp,p,prevpd);			
			double distq = calculateBD(prevfq,fq,prevq,q,prevqd);
			//double distp = abs(p-prevp);
			//double distq = abs(q-prevq);
			//mfile<<"distp ="<<distp<<" distq ="<<distq<<endl;
			//dist = dist/den;
			//double dist2 = sqrt(distp*distp+distq*distq);
			double dist2 = sqrt((distp-distq)*(distp-distq));
			if((j-startcol) != 0){
				 coldist = matrix[i][j]-matrix[i][j-1];
				 coldist2 = matrix[i+1][j]-matrix[i+1][j-1];
			}
			//mfile<<"dist2 ="<<dist2<<" distbd ="<<dist<<endl;
			//dist = (dist*dist2)/(dist+dist2);
			//dist = dist < dist2 ? dist:dist2;
			//mfile<<"dist"<<abs(dist)<<endl;
			//if( abs(dist) <= threshold || dist2 <= threshold){
//			dist = (dist+0.5*dist2)/2;
			if( (abs(dist) <= threshold) && (abs(coldist2-coldist) == 0)){
				if(rowcount == 0 ){
					if(!colcluster.contains(matrix2[i][j]->name))
						colcluster.Append(matrix2[i][j]->name);
				}
				else if(!colcluster.contains(matrix2[i][j]->name)){
					//cerr<<"ovlprow ="<<ovlprow<<" i="<<i<<endl;
					if(i != ovlprow && set == false){
						ovlprow=i+1;
						set = true;
					}
				}
				colcount += 1;
				//mfile<<"column "<<matrix2[i][j]->name<<" added to colclus"<<endl;
				
			}
			//else break;
			if((i == rows-2)&& (rowcount >= rowt && colcount > colt)){
				//mfile<<"i == rows-2 "<<i<<endl;
				rowcluster.Append(matrix3[i+1]->name);
				Bicluster bix(rowcluster, colcluster);
				biclusters[index] = bix;
				//mfile<<"Bix added"<<endl;
			    rowcluster.Purge();		
			    colcluster.Purge();
				index += 1;
				//cerr<<mfile<<"total bix"<<index<<endl;
				rowcount=0;
				colcount =0;
				fp=0;fq=0;p=0;q=0;dist=0;
			}
					
			if((dist > threshold)||(j == (cols-1))){
			//if(j == (cols -1)){
				//prevcol = startcol;
				fp=0;fq=0;p=0;q=0;dist=0;
				prevp=0;prevq=0;prevfp=0;prevfq=0;prevpd=0;prevqd=0;
				//mfile<<"dist > threshold or j= cols-1"<<endl;
				stopcol = j;
				if(colcount >= colt){
					//mfile<<"colcount > colt"<<colcount<<endl;
					if(!rowcluster.contains(matrix3[i]->name)){
						rowcluster.Append(matrix3[i]->name);
						rowcount += 1;
						//mfile<<" row "<<matrix3[i]->name<< " added"<<endl;
					}
					if(!rowcluster.contains(matrix3[i+1]->name)){
						rowcluster.Append(matrix3[i+1]->name);
						rowcount += 1;
						//mfile<<" row "<<matrix3[i+1]->name<< " added"<<endl;
					}
					//mfile<<"Rowcount = "<<rowcount<<"& i ="<<i<<endl;
					if(i == (rows-2)){
						Bicluster bix(rowcluster, colcluster);
						biclusters[index] = bix;
					    rowcluster.Purge();		
					    colcluster.Purge();
						index += 1;
						rowcount =0;
						colcount =0;
					}
					//stopcol = j;
					colcount =0;
					
				}
				else if((rowcount >= rowt)){
				//if((rowcount >= rowt)){
					//mfile<<"rowcount > rowt"<< rowcount<<" "<<i<<endl;
					Bicluster bix(rowcluster, colcluster);
					biclusters[index] = bix;
				    rowcluster.Purge();		
				    colcluster.Purge();
					index += 1;
					rowcount =0;
					colcount =0;
//					if(i < prevrow) startrow=i;
//					else startrow = prevrow;
					prevrow = startrow;
				}
				if(j >= (cols-1)){
					//startrow = i;
					startcol = 0;
					colcount =0;
				}
				else{
					//startcol = stopcol;
					colcount =0;
					prevrow = i;
					if(set == true){
						i = ovlprow;
						set = false;
					}
					//i = i-1;
				}
			}
		}
	}
	// Refine the biclusters
    //cerr<<"Refining biclusters"<<endl;
//	map<int, Bicluster>::iterator iter;
//	for( iter = biclusters.begin(); iter != biclusters.end(); iter++ ) {
//   		 //cout << "word: " << iter->first << ", count: " << iter->second << endl;
//		map<int, Bicluster>::iterator iter2;   		 
//		for( iter2 = biclusters.begin(); iter2 != biclusters.end(); iter2++ ) {
//			if(iter != iter2 && iter->second.first == iter2->second.first){
//				LinkedList<string> tmp = iter2->second.second;
//				tmp += iter->second.second;
//				Bicluster bic(iter->second.first,tmp);
//				map<int, Bicluster>::iterator iter3 = biclusters.find(iter->first);
//				biclusters.erase(iter3);
//				biclusters.insert(make_pair(iter->first,bic));
//				map<int, Bicluster>::iterator iter4 = biclusters.find(iter2->first);
//				biclusters.erase(iter4);
//			}
//		}
//
//    }
//	for(int i=0; i <index; i++){
//		if(biclusters[i] == NULL) skip;
//		for(int j=0; j< index; j++){
//			if(biclusters[j] == NULL) skip;
//			if(j != i && biclusters[i].first == biclusters[j].first){
//				LinkedList<string> tmp = biclusters[j].second;
//				tmp += biclusters[i].first;
//				biclusters[i].first = tmp;
//				map<int, Bicluster>::iterator iter = biclusters.find(j);
//				biclusters.erase(j);
//			}
//		}
//	}
	
	int count =0;
	mfile<<"total bixes ="<<biclusters.size()<<endl;
    ofstream cfile;
	//cout<<"bix file "<<bixfile<<endl;
	cfile.open(fname);
    map<int, Bicluster>::iterator it;
    for(it=biclusters.begin(); it != biclusters.end(); it++){
		     mfile<<"---------------------------"<<endl;
    		 mfile<<"Bicluster #  "<< count++;
    	     pair<LinkedList<string>, LinkedList<string> > bix = it->second;
    	     LinkedList<string> rownodes = bix.first;
    	     LinkedList<string> colnodes = bix.second;
    	     mfile<<" Rows = "<<rownodes.size<<" Cols = "<<colnodes.size<<endl;
    	     //mfile<<" MSR = "<<scores[i]<<endl;
    	     mfile<<"Rows are"<<endl;
    	     for(int rows=0; rows < node0Array.Length(); rows++){
    	     	if(rownodes.contains(node0Array[rows].name)){
    	     		cfile<<"1";
    	     		mfile<<node0Array[rows].name<<" "<<endl;
    	     	}
    	     	else
    	     		cfile<<"0";
    	     	if(rows != (node0Array.Length()-1))
    	     	    cfile << " ";	
    	     }
    	     cfile<<endl;
    	     mfile <<"Columns are::"<<endl;
    		 for(int cols=0; cols < node1Array.Length(); cols++){
    	     	if(colnodes.contains(node1Array[cols].name)){
    	     		cfile<<"1";
    	     		mfile<<node1Array[cols].name<<" "<<endl; 
    	     	}
    	     	else
    	     		cfile<<"0";
    	     	if(cols != (node1Array.Length()-1))
    	     	    cfile << " ";	
    	     }
    	     cfile<<endl;
		//cerr<<"Done"<<endl;
	   
   }
   cfile.close();
   mfile.close();
   
}

void SphierStruct::printMatrix()
{
    int rows = node0Array.Length();
    int cols = node1Array.Length();
    int matrix[rows][cols];
    Node* matrix2[rows][cols];
    Node* matrix3[rows];
    for(int i=0; i<rows;i++){
    	int x = table[0][i];
    	matrix3[x] = &node0Array[i];
    	//cout<<"matrix3[x] for x" << x<<" = "<< matrix3[x]->name<<endl;
        for(int j=0; j<cols; j++){
            int y = table[1][j];
            //cout<<"SS::PM-> x ="<<x<<"  and y="<<y<<endl;
            pair<Node*, int> res = node0Array[i].findNeighbor(&node1Array[j]);
            //cout<<" Neighbor was "<<res.second<<endl;
            if(res.first != NULL){
                matrix[x][y] = res.second;
                matrix2[x][y] = res.first;
               //cout<<"matrix2[x][y] with x = "<< x <<"and y= "<< y <<"="<< matrix2[x][y]->name<<endl;
                   
            }
            else{
                matrix[x][y]=0;
                matrix2[x][y]=0;
            }
        } 
    }
    ofstream mfile;
    mfile.open("output-matrix.cHawk");
    mfile <<"Threshold ="<<threshold<<endl;
   // if(stage == 1){
//    	for(int i=0; i < rows; i++){
//    		for(int j =0; j< cols; j++){
//    			mfile<<matrix[i][j]<<" ";
//    		}
//    		mfile<<endl;
//    	}
    //}
    mfile<<endl;
  
    	int startrow = 0;
		int startcol = 0;
		int stopcol = 0;
		bool add = false;
		bool blocked;
		double residue[cols];
		LinkedList<string> rowcluster;
		LinkedList<string> prevrc;
		LinkedList<string> colcluster;
		LinkedList<string> prevcc;
		typedef pair<LinkedList<string>, LinkedList<string> > Bicluster;
		map<int, Bicluster> biclusters;
		map<int,double> scores;
		double mean_squared_res = 0;
		double prevmsr;
		int colsum[cols];
		int blocksum[cols];
		int rowsum=0;
		int index =0;
		for(int i=0; i < cols; i++){
			colsum[i] =0;
			residue[i]=0;
		}
		for(int i=startrow; i < rows; i++){
			rowsum = 0;
			prevcc = colcluster;
			colcluster.Purge();
			add = false;
			blocked = false;
			if(i == startrow){
				 mfile<<"inside i = rowcluster"<<endl;
				 rowcluster.Append(matrix3[i]->name);
				 mfile<<matrix3[i]->name<<" added to the row cluster"<<endl;
				 //continue;
			}
			for(int j=startcol; j< cols; j++){
				mfile<<"matrix[i][j] = "<< matrix[i][j]<<endl;
				if(matrix[i][j] == 0 || matrix[i][j] == -1)
				    continue;
				mfile<<"i = "<<i<<" and startrow ="<<startrow<<endl;
				if(i == startrow){
					residue[j] =0;
			   		colsum[j] = matrix[i][j];
			   		mfile<<"colsum[j] = "<<colsum[j]<<endl;
			   		blocksum[j] = matrix[i][j];
			   		mfile<<"blocksum[j] = "<<blocksum[j]<<endl;
			   		if(matrix2[i][j]){
			   		   colcluster.Append(matrix2[i][j]->name);
			   		  mfile<<matrix2[i][j]->name<<" column added to col cluster"<<endl;
			   		}
				}
				else{
			   		colsum[j] = colsum[j]+matrix[i][j];
			   		mfile<<"colsum["<<j<<"] = "<<colsum[j]<<endl;
			   		rowsum = rowsum+matrix[i][j];
			   		mfile<<"rowsum "<<rowsum<<endl;
			   		if(j == startrow){
			   	  		blocksum[j] = colsum[j];
			   		}
			  		else{
			   	  		blocksum[j] = blocksum[j-1]+colsum[j];
			   		}
			   		mfile<<"blocksum "<<blocksum[j]<<endl;
			  		int blockcount = (i-startrow+1)*(j-startcol+1);
			  		mfile<<"block count "<<blockcount<<endl;
			   		double rowmean = rowsum/(double)(j-startcol+1);
			   		mfile<<"row mean "<<rowmean<<endl;
			  		double colmean = colsum[j]/double(i-startrow+1);
			  		mfile<<"col mean "<<colmean<<endl;
			   		double blockmean = blocksum[j]/(double)blockcount;
			   		mfile<<"block mean "<<blockmean<<endl;
			   		double res = matrix[i][j]-rowmean-colmean+blockmean;
			   		mfile<<"res = "<<res<<endl;
			   		
			   		if(j == startrow)
			   		    residue[j] = res*res;
			   		else
			   		    residue[j] = residue[j-1]+(res*res);
			   		    
			   		mfile<<"residue["<<j<<"] = "<<residue[j]<<endl;
			   		double msr = residue[j]/(double)blockcount;
			   			
			   		mfile<<"msr = "<<msr<<endl;
	   			    mfile<<"threshold = "<<threshold<<endl;
			   		if(msr <= threshold){
			   		    mfile<<"msr < threshold"<<endl;
			   	  		if(!colcluster.contains(matrix2[i][j]->name)){
				   	  		colcluster.Append(matrix2[i][j]->name);
				    			mfile<<matrix2[i][j]->name<<" is added into col cluster"<<endl;
				    	  	if(msr != 0){
				   	    	    mean_squared_res = msr;
				   	             mfile<<"mean squared res = "<<mean_squared_res<<endl;
				   	        }
			   	  	    }
			   	  		add = true;
			   	  		mfile<<"add = true"<<endl;
			   		} 
			   		else{
			   				//if(add == true){
			   		    stopcol = j;
			   		    blocked = true;
			   		    mfile<<"blocked = true"<<j<<endl;
			   		    break;
			   				//}
			   		}
			   
				}			    
		}//end for j
		if(i != startrow){
			if(add ==true ){
				rowcluster.Append(matrix3[i]->name);
				mfile<<matrix3[i]->name<<" added to rowcluster"<<endl;
				add = false;
				
			}
		}
	    if(blocked == true){
		   	mfile<<"Bicluster# "<<index<<" created with rows ="<<rowcluster.size;
		   	mfile<<"and columns = "<<colcluster.size<<endl;
		   	if( rowcluster.size >= rowt && prevcc.size >= colt ){
		    	Bicluster tmp(rowcluster, prevcc);
		    	biclusters[index] = tmp;
		    	scores[index] = mean_squared_res;
		    	//cout<<"in bix msr ="<<scores[index];
		    	index = index+1;
	         }
	         prevrc = rowcluster;
		     rowcluster.Purge();
		     colcluster.Purge();
		     startrow = i;
		  	 i = i-1;
		    	if(stopcol > 0)
		      		startcol = stopcol;
		   		else
		      		startcol = 0;
		        for(int c=0; c< cols; c++){
		           colsum[c] = 0;
		           blocksum[c]=0;
		           residue[c] = 0;
		        }
		        rowsum =0;
		    	//cout<<"startcol = "<<startcol<<endl;
		    	//cout<<"startrow = "<<startrow<<endl;
		    
		}
		if((rowcluster.size >= rowt)&& (i == rows-2)){
			    
		    	Bicluster tmp(rowcluster, colcluster);
		    	//rowcluster.Purge();
		    	//colcluster.Purge();
		    	biclusters[index] = tmp;
		    	scores[index] = mean_squared_res;
			    mfile<<"in bix msr ="<<scores[index];
		    	index = index+1;
		}	
		prevrc = rowcluster;
	}// end for i
	   
	ofstream cfile;
	//cout<<"bix file "<<bixfile<<endl;
	cfile.open(fname);
	for(int i=0; i < index ; i++){
		     mfile<<"---------------------------"<<endl;
    		 mfile<<"Bicluster #  "<< i;
    	     pair<LinkedList<string>, LinkedList<string> > bix = biclusters[i];
    	     LinkedList<string> rownodes = bix.first;
    	     LinkedList<string> colnodes = bix.second;
    	     mfile<<" Rows = "<<rownodes.size<<" Cols = "<<colnodes.size;
    	     mfile<<" MSR = "<<scores[i]<<endl;
    	     mfile<<"Rows are"<<endl;
    	     for(int rows=0; rows < node0Array.Length(); rows++){
    	     	if(rownodes.contains(node0Array[rows].name)){
    	     		cfile<<"1";
    	     		mfile<<node0Array[rows].name<<" "<<endl;
    	     	}
    	     	else
    	     		cfile<<"0";
    	     	if(rows != node0Array.Length()-1)
    	     	    cfile << " ";	
    	     }
    	     cfile<<endl;
    	     mfile <<"Columns are::"<<endl;
    		 for(int cols=0; cols < node1Array.Length(); cols++){
    	     	if(colnodes.contains(node1Array[cols].name)){
    	     		cfile<<"1";
    	     		mfile<<node1Array[cols].name<<" "<<endl; 
    	     	}
    	     	else
    	     		cfile<<"0";
    	     	if(cols != node1Array.Length()-1)
    	     	    cfile << " ";	
    	     }
    	     cfile<<endl;
	   
   }
   cfile.close();
   mfile.close();
   
}

char*
SphierStruct::itoa( int value, char* result, int base ) {
	
	// check that the base if valid
	
	if (base < 2 || base > 16) { *result = 0; return result; }
	char* out = result;
	int quotient = value;
	do {
    	*out = "0123456789abcdef"[ std::abs( quotient % base ) ];
    	++out;
    	quotient /= base;
	
	} while ( quotient );
	// Only apply negative sign for base 10
	if ( value < 0 && base == 10) *out++ = '-';
	std::reverse( result, out );
	*out = 0;
	return result;
}

/*int main(){
   SphierStruct spt(5,5);
   for (int i=0; i< 5; i++){
        string str1(i,'n');    
         string str2(i,'m');
       spt.addNode(str1,0);
       spt.addNode(str2,1);
   } 
   for(int i=0; i< 5; i++){
      Node* temp1 = spt.findNode(i,0);
      Node* temp2 = spt.findNode(i,1);
      Node* temp3 = spt.findNode((i+1)%5,1);
      temp1->addNeighbor(*temp2);
      temp1->addNeighbor(*temp3);
      temp2->addNeighbor(*temp1);
      temp3->addNeighbor(*temp1);
   }
   spt.printMatrix();
   spt.printStruct();
   for(int i=0; i<3; i++){
      spt.setPositionBary();
   }
   spt.printMatrix();
}*/
