costToGoMatrix[0][Q]= expectedCost; //costToGoMatrix

  return costToGo_j[Q]; 
}



//Given the cost-to-go vector from succ_j th customer, compute cost-to-go from j,
//in the hypothesis of going first to the depot for restocking
//and then to succ_jth customer.  
double 
Solution::computeCostIfPreventiveRestock(int j,int succ_j,const vector<double>& cost_succ_j){
  double temp = 0.0;
  vector<PossibleDemand>::iterator d;
  for(d= problem->customerDemand[(*this)[succ_j]].begin();
      d< problem->customerDemand[(*this)[succ_j]].end();
      d++ ){
    temp += cost_succ_j[problem->capacity - d->demand_ ] * d->probability_ ;
  }

  temp +=  problem->distanceMatrix[(*this)[j]][0] +
    problem->distanceMatrix[0][(*this)[succ_j]];
  //  return  temp + problem->distanceMatrix[(*this)[j]][0] +
  //  problem->distanceMatrix[0][(*this)[succ_j]];
  //  cout <<"computeCostIfPreventiveRestock()= " << temp << endl;//////////////////

  return temp;
}



//Given the cost-to-go vector from succ_j th, compute cost-to-go from jth,
//in the hypothesis of proceeding directly from the jth customer
//to the succ_jth customer.  
double 
Solution::computeCostIfProceed(int q, int j,int succ_j,const vector<double>& cost_succ_j){

  double temp = 0.0;
  vector<PossibleDemand>::iterator d;
  for(d= problem->customerDemand[(*this)[succ_j]].begin();
      d< problem->customerDemand[(*this)[succ_j]].end();
      d++ ){
    if(d->demand_ <= q)
      temp += cost_succ_j[q - d->demand_ ] * d->probability_ ;
    else 
      temp += (problem->distanceMatrix[(*this)[succ_j]][0] + 
	       problem->distanceMatrix[0][(*this)[succ_j]] + 
	       cost_succ_j[q + problem->capacity - d->demand_ ] ) * d->probability_;
  }

  temp += problem->distanceMatrix[(*this)[j]][(*this)[succ_j]];
  // cout <<"computeCostIfProceed(): q "<< q << " temp " << temp << endl ;/////////////
  return  temp;
}



//Compute the approximated variation in expected cost 
//due to the application of the current Oropt move to this solution
double 
Solution::computeProxyDelta(){
  int i = move[0];
  int k = move[1];
  int j = move[2];
  
  double PDS = computeProxyDeletionSaving(costToGoMatrix[i], costToGoMatrix[i+k+1] );
  double PIC =computeProxyInsertionCost(costToGoMatrix[j], 
					costToGoMatrix[(j+1)%(problem->numberOfCustomers)] );
  //cout <<"computeProxyDelta(): PIC= " << PIC << " PDS= " << PDS << endl  ;////////////
  //cout <<"computeProxyDelta()= " << PIC-PDS << endl;///////////
  return PIC - PDS;
}



//Compute the approximated saving  for deleting a string
//of consecutive customers, as given by the current Oropt
//move of this solution.
//Input parameter is the cost-to-go vector from the first node after 
//the string in the current move.
double
Solution::computeProxyDeletionSaving(const vector<double>& cost_prev ,
				     const vector<double>& cost_next ){
  int  Q = problem->capacity;
  double average=0.0;
  int i = move[0];
  int k = move[1];
  double cost_proceed;
  double cost_pruned;
  double cost_restock =  computeCostIfPreventiveRestock( i, i+k+1, cost_next );

  cost_pruned = cost_restock;
  if( i==0 ) average = cost_prev[Q]-cost_pruned; 
  
  else{
    for(int q=Q; q>=0; q--){
      cost_pruned = cost_restock;
      cost_proceed = computeCostIfProceed( q, i, i+k+1, cost_next ) ;
      if( !(cost_proceed > cost_restock) ) cost_pruned = cost_proceed;
      
      average += cost_prev[q] - cost_pruned; //update average
    }

    average = average/(double)(Q+1);
  }
  return average;
}



//Compute the approximated cost for inserting a string
//of consecutive customers, as given by the current Oropt
//move of this solution.
//Input parameter is the cost-to-go vector from the customer 
//after which the string would be inserted, according to this solution 
//Oropt move.
double 
Solution::computeProxyInsertionCost(const vector<double>& cost_prev,
				    const vector<double>& cost_next){
  double average=0.0;
  int n = problem->numberOfCustomers;
  int  Q = problem->capacity;
  int j = move[2];
  int i = move[0];
  int k = move[1]; //String length.

  //Initialization of cost-to-go vector as initial data of the Dynamic programming
  //iteration.
  vector<double> costToGo_jplus1(Q+1);
  for(int q=0; q<=Q; q++){
      if( j < n -1 ) costToGo_jplus1[q] = cost_next[q];
      else
	costToGo_jplus1[q]=problem->distanceMatrix[ (*this)[i+k] ][0];
    }
  //cout <<"init costToGo_jplus1[0]"<<  costToGo_jplus1[0] << endl;//////////////

  //Number of Dynamic programming iterations.
  int iter;
  if(j < n-1) iter = k;
  else iter = k-1;

  //Dynamic programming recursion.
  double costRestock;
  double costNoRestock;
  vector<double> costToGo_j(Q+1);

  for(int l=iter; l>=1; l--){
    int next = -1;
    if(l == k) next = (j+1)%n;    //Next is first customer after the inserted string.
    else next = i+l+1;            //Next is inside the string.
    costRestock = computeCostIfPreventiveRestock(i+l,next,costToGo_jplus1);
    //cout << i+l <<" h' " << costRestock << endl;//////////////
    
    for(int q=Q; q>=0; q--){
      //compute costToGo in case of proceeding to the next customer
      costNoRestock = computeCostIfProceed(q,i+l,next,costToGo_jplus1);

      //cout << j <<" h(" << q << ")= "<< costNoRestock ;//////////////
      
      //choose smaller cost btw preventiveRestock and Proceed, and set threshold
      if(costRestock < costNoRestock)
	costToGo_j[q]= costRestock;
      else 
	costToGo_j[q] = costNoRestock;
      // cout << i+l <<" f(" << q << ")= "<<  costToGo_j[q] << endl;//////////////
    }
    
    //store chosen cost vector in costToGo_jplus1   
    for(int q=Q; q>=0; q--)
      costToGo_jplus1[q] = costToGo_j[q];
  }  
  
  //customer = j: costToGo from the customer after which the string is inserted
  costRestock = computeCostIfPreventiveRestock(j,i+1,costToGo_jplus1);
  for(int q=Q; q>=0; q--){
    costNoRestock = computeCostIfProceed(q,j,i+1,costToGo_jplus1);
    if(costRestock < costNoRestock)
      costToGo_j[q]= costRestock;
    else 
      costToGo_j[q] = costNoRestock;
    //cout << j <<" f(" << q << ")= " << costToGo_j[q] << endl;/////////////
    //    cout  <<" cost_prev(" << q << ")= " << cost_prev[q] << endl;/////////////
    average+= costToGo_j[q]-cost_prev[q];
  }

  average = average/(Q+1);

  return average;
}



double 
Solution::computeExactDelta(){
  int n = problem->numberOfCustomers;
  int  Q = problem->capacity;
  int j = move[2];
  int i = move[0];
  int k = move[1]; //String length.

  // cout << "compExactDelta(): " << i << k << j << endl;////////////

  //Initialization of cost-to-go vector, either from j+1 or from the
  //last customer of the string.
  vector<double> costToGo_jplus1(Q+1);
  for(int q=0; q<=Q; q++){
    if( j < n -1 ){
      costToGo_jplus1[q] = costToGoMatrix[j+1][q];
      //    cout << j+1 << " " << (j+2)%n << endl;//////////////
    }
    else{
      costToGo_jplus1[q]=problem->distanceMatrix[ (*this)[i+k] ][0];
      //cout << i+k << " " << 0 << endl;//////////////
    }
  }
  

  //Number of Dynamic programming iterations.
  int iter;
  if(j < n-1) iter = k;
  else iter = k-1;

  //Recursion backwards in the string.
  double costRestock;
  double costNoRestock;
  vector<double> costToGo_j(Q+1);

  for(int l=iter; l>=1; l--){
    int next = -1;
    if(l == k) next = (j+1)%n;    //Next is first customer after the inserted string.
    else next = i+l+1;            //Next is inside the string.
    costRestock = computeCostIfPreventiveRestock(i+l,next,costToGo_jplus1);
    //  cout << i+l << " " << next << endl;//////////////
    
    for(int q=Q; q>=0; q--){
      //compute costToGo in case of proceeding to the next customer
      costNoRestock = computeCostIfProceed(q,i+l,next,costToGo_jplus1);

      //cout << j <<" h(" << q << ")= "<< costNoRestock ;//////////////
      
      //choose smaller cost btw preventiveRestock and Proceed, and set threshold
      if(costRestock < costNoRestock)
	costToGo_j[q]= costRestock;
      else 
	costToGo_j[q] = costNoRestock;
      // cout << i+l <<" f(" << q << ")= "<<  costToGo_j[q] << endl;//////////////
    }
    
    //store chosen cost vector in costToGo_jplus1   
    for(int q=Q; q>=0; q--)
      costToGo_jplus1[q] = costToGo_j[q];
  }  
  
  //Recursion from the beginning of the string back to j.
  costRestock = computeCostIfPreventiveRestock(j,i+1,costToGo_jplus1);
  for(int q=Q; q>=0; q--){
    costNoRestock = computeCostIfProceed(q,j,i+1,costToGo_jplus1);
    if(costRestock < costNoRestock)
      costToGo_j[q]= costRestock;
    else 
      costToGo_j[q] = costNoRestock;
  }
  //store chosen cost vector in costToGo_jplus1   
  for(int q=Q; q>=0; q--)
    costToGo_jplus1[q] = costToGo_j[q];

  //cout << j << " " << i+1 << endl;//////////////

  //Recursion from j back to  i+k+1.
  for(int l=j; l>i+k+1; l--){
    costRestock = computeCostIfPreventiveRestock(l-1,l,costToGo_jplus1);
    //cout << l-1 << " " << l << endl;//////////////
    
    for(int q=Q; q>=0; q--){
      //compute costToGo in case of proceeding to the next customer
      costNoRestock = computeCostIfProceed(q,l-1,l,costToGo_jplus1);

      //cout << j <<" h(" << q << ")= "<< costNoRestock ;//////////////
      
      //choose smaller cost btw preventiveRestock and Proceed, and set threshold
      if(costRestock < costNoRestock)
	costToGo_j[q]= costRestock;
      else 
	costToGo_j[q] = costNoRestock;
      // cout << i+l <<" f(" << q << ")= "<<  costToGo_j[q] << endl;//////////////
    }
    
    //store chosen cost vector in costToGo_jplus1   
    for(int q=Q; q>=0; q--)
      costToGo_jplus1[q] = costToGo_j[q];
  }  

  //Last itaration from i+k+1 back to i.
  if(i==0) {//Last iteration.
    costToGo_j[Q] = computeCostIfProceed(Q,i,i+k+1,costToGo_jplus1);
    //cout << i << " " << i+k+1 << endl;//////////////
    return costToGo_j[Q]-expectedCost;
  }

  else{//i>0.
    //Single itaration from i+k+1 back to i.
    costRestock = computeCostIfPreventiveRestock(i,i+k+1,costToGo_jplus1);
    // cout << i << " " << j << endl;//////////////
    for(int q=Q; q>=0; q--){
      costNoRestock = computeCostIfProceed(q,i,i+k+1,costToGo_jplus1);
      if(costRestock < costNoRestock)
      costToGo_j[q]= costRestock;
      else 
	costToGo_j[q] = costNoRestock;
    }
    //store chosen cost vector in costToGo_jplus1   
    for(int q=Q; q>=0; q--)
      costToGo_jplus1[q] = costToGo_j[q];

    //Iterations from i back to the depot.
    for(int cr=i-1; cr>=1; cr-- ){
      costRestock = computeCostIfPreventiveRestock(cr,cr+1,costToGo_jplus1);
      //cout << j <<" h' " << costRestock << endl;//////////////
      //cout << cr << " " << cr+1 << endl;//////////////	
      
      for(int q=Q; q>=0; q--){
	//compute costToGo in case of proceeding to the next customer
	costNoRestock = computeCostIfProceed(q,cr,cr+1,costToGo_jplus1);
	//cout << j <<" h(" << q << ")= "<< costNoRestock ;//////////////
	
	//choose smaller cost btw preventiveRestock and Proceed, and set threshold
	if(costRestock < costNoRestock)
	  costToGo_j[q]= costRestock;
	else 
	  costToGo_j[q] = costNoRestock;
	//cout << " f(" << q << ")= "<<  costToGo_j[q] << endl;//////////////
      }
      
      //store chosen cost vector in costToGo_jplus1   
      for(int q=Q; q>=0; q--)
	costToGo_jplus1[q] = costToGo_j[q];
    }  
  
    //cr==0: costToGo from the depot on (last iteration)
    costToGo_j[Q] = computeCostIfProceed(Q,0,1,costToGo_jplus1);
    //cout << 0 << " " << 1 << endl;//////////////
    //cout << " costToGo_j[Q]= " << costToGo_j[Q] << endl;/////////////
    
    return costToGo_j[Q]-expectedCost;
    
  }

}



//Print on output stream the solution 
void
Solution::printOn( ostream &os ) {
	os << "#sequence of "<< size()<< " customers:" << endl;
	for( iterator i = begin(); i != end(); i++ )
	  //os << (*this)[*i] << " ";
	  os << *i << " ";
	os << endl;
	
}



//Print on output stream the solution and the thresholds
void
Solution::printOn( ostream &os,const vector<int>& thresholds ) {
	os << "#sequence of "<< thresholds.size()<< " couples (customer,threshold):"
	   << endl;
	for( iterator i = begin(); i != end(); i++ )
		os << *i << "," << thresholds[*i] << "   ";
	os << endl;
	
}


//Print on output stream the costToGo matrix given in input
void
Solution::printOn( ostream &os,const vector<vector<double> >& costToGoMatrix ) {
  if (!computedExpectedCost) cout <<"PrintOn():Warning! printing not up-to-date data! " << endl;
	os << "#costToGo matrix" << endl;

	for( int i =0 ; i <= problem->numberOfCustomers-1; i++ ){
	  cout << i <<"th customer: " ;
	  for(int q = 0; q <= problem->capacity ; q++)
	    cout << costToGoMatrix[i][q] << " ";
	  cout << endl;
	}
}



//Copy a Solution into this one.
void
Solution::copySolution(const Solution& orig ) {
  //(*this) = orig;
  for(int i=0; i< problem->numberOfCustomers; i++) (*this)[i] = orig[i];
  for(int s=0; s<3; s++) move[s] = orig.move[s];
  expectedCost = orig.expectedCost;
  computedExpectedCost = orig.computedExpectedCost;
  rg = orig.rg;
  control_ = orig.control_;
}



//Get the pointer to the problem.
Problem*
Solution::getProblem( ) {
  return problem;
}



//Get the current move of the solution
vector<int>
Solution::getMove( ) {
  return move;
}



void
Solution::allocateCostToGoMatrix(){
  int n = problem->numberOfCustomers;
  int Q = problem->capacity;

  costToGoMatrix = vector< vector<double> > ( n );
  for( int i=0; i< n ; i++ )
    costToGoMatrix[i] = vector<double>( Q+1, -1.0 );
}