e_increase++;
		}
		quality.energy[5] = e_increase;
		quality.quality[5] = (next_Kfree + 0.5)*pow(parameter.Exp, e_increase);
	}

	for(int i=0; i<6; i++)
		quality.totalQuality += quality.quality[i];
	quality.averageQuality = quality.totalQuality / quality.k_free;
}


/***when restart simulation we must do this method***/ 
void Simulation::putFirstTwoPoint(){
	Point tempPoint(-1,-1,-1);
	int temp = sequence.getLength()-1;
	tempPoint.set(temp, temp, temp);
	tempResult.push_back(tempPoint);
	tempPoint.set(temp+1,temp,temp);
	tempResult.push_back(tempPoint);
	tempWeight.push_back(1);
	tempWeight.push_back(1);
	tempEnergy.push_back(0);
	tempEnergy.push_back(0);
	setSpaceNode(temp,temp,temp,sequence.getSequence().GetAt(0));
	setSpaceNode(temp+1,temp,temp,sequence.getSequence().GetAt(1));
}


/***choose a direction to put  node in***/
int Simulation::selectOneDirection(Quality &tempQuality){
	long double operation;
	int i;
	long double partSum[7];
	partSum[0] = 0;
	for(i=1; i<7; i++)
		partSum[i] = partSum[i-1] + tempQuality.quality[i-1]/tempQuality.totalQuality;
	while((operation = (long double)(rand()%1000)/1000) == 0){
	}
	for(i=1; i<7; i++){
		if(operation <= partSum[i] && operation > partSum[i-1]){
			break;			
		}			
	}	
	return (i-1);
}

/***put the  next point in to the space and set some parameter***/
void Simulation::putOneNode(Quality &tempQuality, int selectedDire){
	long double weight;
	Point point(-1,-1,-1);
	int energy =  tempQuality.energy[selectedDire];
	weight = tempWeight.back()*pow(parameter.Exp, energy)*tempQuality.k_free;
	tempWeight.push_back(weight);
	tempEnergy.push_back(energy);
	tempTotalEnergy += energy;
	//put the point into the site;
	point.set(tempResult.back().x+parameter.directions[selectedDire][0],
			  tempResult.back().y+parameter.directions[selectedDire][1],
			  tempResult.back().z+parameter.directions[selectedDire][2]);
	tempResult.push_back(point);
	setSpaceNode(point.x,point.y,point.z,sequence.getSequence().GetAt(tempResult.size()-1));
}

/***initialize a structure***/
void Simulation::initialization(){
	Quality quality;
	Point tempPoint(-1,-1,-1); //we use the point to insert the vector
	int select;
	//begin the construction
	while(tempResult.size() != sequence.getLength()){
		resetTemp();
		resetMatrix();
		putFirstTwoPoint();
		//find the sturction and give up the temp structure when a collide happen
		while(tempResult.size() < sequence.getLength()){
			computingQuality(quality);
			if(quality.k_free == 0)
				break;
			else{
				select = selectOneDirection(quality);
				putOneNode(quality, select);
			}
		}
	}
	bestRecoder = tempTotalEnergy;
	bestResult = tempResult;
	parameter.initial(tempWeight);
	
#ifdef MYDEBUG
/*
		std::vector<Point>::iterator pos;
		std::vector<long double>::iterator weight = parameter.averageOfWei.begin();
		for(pos = tempResult.begin(); pos<tempResult.end();pos++){
			cout<<pos->x<<"  "<<pos->y<<"  "<<pos->z<<" "<<"Weight:"<<*weight<<endl;
			weight++;
		}
		cout<<"energy of the structure:"<<tempTotalEnergy<<endl;
		cout<<"<The initialization of the simulaiton is stopping>"<<endl;*/
	
#endif
}

/***back to the last state**/
void Simulation::releaseOneNode(){
	tempWeight.pop_back();//delete
	int energy = tempEnergy.back();
	tempTotalEnergy -= energy;//delete
	tempEnergy.pop_back();
	Point point = tempResult.back();
	tempResult.pop_back();//delete
	setSpaceNode(point.x, point.y, point.z, '0');//delete
}

	
long double Simulation::selectMultiDirection(Quality &tempQuality, std::vector<int> &direction, const int &number){
	int multiple = linage(tempQuality.k_free, number)*number / tempQuality.k_free;
	long double allQuality = tempQuality.totalQuality * multiple;
	long double partQuality = 0;
	int tempIndex;
	int flag[6] = {0,0,0,0,0,0};
	while(direction.size() < number){
		tempIndex = rand()%6;
		if(tempQuality.quality[tempIndex] != 0 && flag[tempIndex] ==0){
			flag[tempIndex] = 1;
			direction.push_back(tempIndex);
			partQuality += tempQuality.quality[tempIndex];
		}
	}
	return (allQuality/partQuality);
}

/***conputing c(N,m)***/
int Simulation::linage(int big, int little){
	static const int mul[6] = {1, 1, 2, 6 ,24, 120};
	return (mul[big] / (mul[little]*mul[big-little]));
}

/***start the recursion procedure***/	
void Simulation::recursion(){
#ifdef MYDEBUG
		//cout<<"The program begin???"<<1<<endl;
#endif
	int finished = tempResult.size();
	int select;
	//find the better result, the reset the best recorder
	if(finished == sequence.getLength()){
#ifdef MYDEBUG
		//cout<<"((((((There is a new structure)))))))))"<<endl;
#endif		
		if(tempTotalEnergy >= bestRecoder){
			bestRecoder = tempTotalEnergy;
			bestResult = tempResult;
		}
		return;
	}

	Quality tempQuality;
	computingQuality(tempQuality);
	//the structure have collided
	if(tempQuality.k_free == 0)
		return;

	long double weightPred  = tempWeight.back()*tempQuality.totalQuality;
#ifdef MYDEBUG
		//cout<<parameter.averageOfWei.size()<<endl;
		//cout<<parameter.c<<endl;
		//cout<<parameter.averageOfWei[finished]<<endl;
		//cout<<parameter.numberOfLen.size()<<endl;
#endif
	long double uperThreshold 
		= parameter.c *(parameter.averageOfWei[finished] / parameter.z0)
					  *(parameter.numberOfLen[finished] / parameter.c0)
					  *(parameter.numberOfLen[finished] / parameter.c0);
	long double downThreshold = uperThreshold * 0.2;

	if(weightPred >= downThreshold && weightPred < uperThreshold){
		select = selectOneDirection(tempQuality);
		putOneNode(tempQuality, select);

#ifdef MYDEBUG
	//cout<<"<*********Enter the middle area***********>"<<endl;
	/*cout<<"Direciton selected is "<<select<<endl;
	cout<<"TempTotalEnergy is "<<tempTotalEnergy<<endl;*/		
#endif

		parameter.updatePara(tempWeight.back(), finished);
		recursion();
		releaseOneNode();
	}
	else if(weightPred < downThreshold){
#ifdef MYDEBUG

	//cout<<"<**********Enter the low area*********>"<<endl;
	/*cout<<"Direciton selected is "<<select<<endl;
	cout<<"TempTotalEnergy is "<<tempTotalEnergy<<endl;*/
		
#endif
		if((double)(rand()%1000) / 1000 <0.5)
			return;
		else{
			select = selectOneDirection(tempQuality);
			putOneNode(tempQuality, select);
			int tail = tempWeight.size() -1;
			tempWeight[tail] = 2*tempWeight[tail]; 
			parameter.updatePara(tempWeight.back(), finished);
			recursion();
			releaseOneNode();
		}
	}
	else if(weightPred > uperThreshold){
		unsigned int numberOfCopy = (int)(weightPred/uperThreshold)+1;
		if(numberOfCopy > tempQuality.k_free)
			numberOfCopy = tempQuality.k_free;
		std::vector<int> direction;
		//we need selectMultiDIrection return a random sequence direcitons !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		double PA = selectMultiDirection(tempQuality, direction, numberOfCopy);
#ifdef MYDEBUG

		//cout<<"<*********Enter the up area**********>"<<endl;
/*		cout<<"Number of the direciton selected are "<<direction.size()<<endl;
		for(int j=0; j<direction.size(); j++)
			cout<<"The direction of selected is: "<<direction[j]<<endl;
		cout<<"TempTotalEnergy is "<<tempTotalEnergy<<endl;*/
		
#endif
		//start the loop
		for(int i=0; i<numberOfCopy; i++){
			putOneNode(tempQuality, direction[i]);
			int tail = tempWeight.size() -1;
			tempWeight[tail] = tempWeight[tail] / (numberOfCopy*PA*linage(tempQuality.k_free, numberOfCopy)); 			
			parameter.updatePara(tempWeight.back(), finished);
			recursion();
			releaseOneNode();
		}
	}
}