/**
 * Solves (Solomon's benchmark instance to the) Vehicle Routing Problem
 * 
 * @author (Stefan Edelkamp) 
 * @version (2013)
 */

import java.util.Random;  

public class VRP
{
  public final static int N=101;
  public final static int V=19; 
  public final static int ITERATIONS=40;  
  public final static int LEVEL=6; 
  public final static int SHOW=(LEVEL-3);  
  public final static double MAXVALUE=100000000.0;  

  double max_capacity = 200.0;        // Benchmark instance r101 
  double [] coordx = {35.0,
     41.0,35.0,55.0,55.0,15.0,25.0,20.0,10.0,55.0,30.0,
     20.0,50.0,30.0,15.0,30.0,10.0, 5.0,20.0,15.0,45.0,
     45.0,45.0,55.0,65.0,65.0,45.0,35.0,41.0,64.0,40.0,
     31.0,35.0,53.0,65.0,63.0, 2.0,20.0, 5.0,60.0,40.0,
     42.0,24.0,23.0,11.0, 6.0, 2.0, 8.0,13.0, 6.0,47.0,
     49.0,27.0,37.0,57.0,63.0,53.0,32.0,36.0,21.0,17.0,
     12.0,24.0,27.0,15.0,62.0,49.0,67.0,56.0,37.0,37.0,
     57.0,47.0,44.0,46.0,49.0,49.0,53.0,61.0,57.0,56.0,
     55.0,15.0,14.0,11.0,16.0, 4.0,28.0,26.0,26.0,31.0,
     15.0,22.0,18.0,26.0,25.0,22.0,25.0,19.0,20.0,18.0};
  double [] coordy = {35.0,
     49.0,17.0,45.0,20.0,30.0,30.0,50.0,43.0,60.0,60.0,
     65.0,35.0,25.0,10.0, 5.0,20.0,30.0,40.0,60.0,65.0,
     20.0,10.0, 5.0,35.0,20.0,30.0,40.0,37.0,42.0,60.0,
     52.0,69.0,52.0,55.0,65.0,60.0,20.0, 5.0,12.0,25.0,
      7.0,12.0, 3.0,14.0,38.0,48.0,56.0,52.0,68.0,47.0,
     58.0,43.0,31.0,29.0,23.0,12.0,12.0,26.0,24.0,34.0,
     24.0,58.0,69.0,77.0,77.0,73.0, 5.0,39.0,47.0,56.0,
     68.0,16.0,17.0,13.0,11.0,42.0,43.0,52.0,48.0,37.0,
     54.0,47.0,37.0,31.0,22.0,18.0,18.0,52.0,35.0,67.0,
     19.0,22.0,24.0,27.0,24.0,27.0,21.0,21.0,26.0,18.0};
  double [] l = {0.0,
    161.0, 50.0,116.0,149.0, 34.0, 99.0, 81.0, 95.0, 97.0,124.0,
     67.0, 63.0,159.0, 32.0, 61.0, 75.0,157.0, 87.0, 76.0,126.0, 
     62.0, 97.0, 68.0,153.0,172.0,132.0, 37.0, 39.0, 63.0, 71.0, 
     50.0,141.0, 37.0,117.0,143.0, 41.0,134.0, 83.0, 44.0, 85.0,
     97.0, 31.0,132.0, 69.0, 32.0,117.0, 51.0,165.0,108.0,124.0, 
     88.0, 52.0, 95.0,140.0,136.0,130.0,101.0,200.0, 18.0,162.0,
     76.0, 58.0, 34.0, 73.0, 51.0,127.0, 83.0,142.0, 50.0,182.0,
     77.0, 35.0, 78.0,149.0, 69.0, 73.0,179.0, 96.0, 92.0,182.0,
     94.0, 55.0, 44.0,101.0, 91.0, 94.0, 93.0, 74.0,176.0, 95.0,
    160.0, 18.0,188.0,100.0, 39.0,135.0,133.0, 58.0, 83.0,185.0};        
   double [] r = {230,
    171.0, 60.0,126.0,159.0, 44.0,109.0, 91.0,105.0,107.0,134.0,
     77.0, 73.0,169.0, 42.0, 71.0, 85.0,167.0, 97.0, 86.0,136.0, 
     72.0,107.0, 78.0,163.0,182.0,142.0, 47.0, 49.0, 73.0, 81.0, 
     60.0,151.0, 47.0,127.0,153.0, 51.0,144.0, 93.0, 54.0, 95.0,
    107.0, 41.0,142.0, 79.0, 42.0,127.0, 61.0,175.0,118.0,134.0, 
     98.0, 62.0,105.0,150.0,146.0,140.0,111.0,210.0, 28.0,172.0,
     86.0, 68.0, 44.0, 83.0, 61.0,137.0, 93.0,152.0, 60.0,192.0,
     87.0, 45.0, 88.0,159.0, 79.0, 83.0,189.0,106.0,102.0,192.0,
    104.0, 65.0, 54.0,111.0,101.0,104.0,103.0, 84.0,186.0,105.0,
    170.0, 28.0,198.0,110.0, 49.0,145.0,143.0, 68.0, 93.0,195.0};        
  double [] w = {0.0,
    10.0, 7.0,13.0,19.0,26.0, 3.0, 5.0, 9.0,16.0,16.0,
    12.0,19.0,23.0,20.0, 8.0,19.0, 2.0,12.0,17.0, 9.0,
    11.0,18.0,29.0, 3.0, 6.0,17.0,16.0,16.0, 9.0,21.0,
    27.0,23.0,11.0,14.0, 8.0, 5.0, 8.0,16.0,31.0, 9.0,
     5.0, 5.0, 7.0,18.0,16.0, 1.0,27.0,36.0,30.0,13.0,
    10.0, 9.0,14.0,18.0, 2.0, 6.0, 7.0,18.0,28.0, 3.0,
    13.0,19.0,10.0, 9.0,20.0,25.0,25.0,36.0, 6.0, 5.0,
    15.0,25.0, 9.0, 8.0,18.0,13.0,14.0, 3.0,23.0, 6.0,
    26.0,16.0,11.0, 7.0,41.0,35.0,26.0, 9.0,15.0, 3.0, 
     1.0, 2.0,22.0,27.0,20.0,11.0,12.0,10.0, 9.0,17.0};        
  double [] service = {0,
    10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,
    10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,
    10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,
    10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,
    10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,
    10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,
    10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,
    10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,
    10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,
    10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0,10.0};  
    
  int [] visits = new int[N+1];     
  int [] tour = new int[V+N+1];       
  int tourSize;     
  double sumservice = 0.0; 
  double [][] d = new double[N][N];   // distance matrix
  int [] moves = new int[N];          // static array for successor
  double [] value = new double[N];    // static array for roulette wheel selection 
  double [][] global = new double[N][N];          // global tour policy 
  double [][][] backup = new double[LEVEL+1][N][N]; // backup of global tour 
  long evaluations = 0;      // counting number of runs/rollouts 
  Random random;
  
  class Pair { double score; int [] tour = new int[N+V+1]; }
  
  public VRP() {
    random = new Random(1019281); 
    for (int i = 0 ; i < N; i++) { 
      sumservice += service[i];
      for (int j = 0 ; j < N; j++) {
        d[i][j] = service[i] - 0.0000001 + Math.sqrt(
	       (coordx[i] - coordx[j]) * (coordx[i] - coordx[j]) +
	       (coordy[i] - coordy[j]) * (coordy[i] - coordy[j]));  
      }
    }
    for(int i = 0; i < N; i++) 
      for(int j = 0; j < N; j++) 
        global[i][j] = 0.0;
    for (int i = 0; i<N; i++) 
      for (int j = 0; j<N; j++) 
        if (l[i] + d[i][j] > r[j]) 
          d[i][j] = 1000000.0; 
    for (int i = 0; i< N;i++) {
      double dsum = 0.0;
      int count = 0;
      for (int j = 0; j< N;j++) 
        if (d[i][j] != 1000000.0) {
          dsum += d[i][j];
          count++;
        }
      double aver = dsum/count;
      for (int j = 0; j< N;j++) {
        if (d[i][j]!= 1000000.0) {
          global[i][j] = (aver/d[i][j])/N; 
        }
        if (d[i][j] == 1000000.0) 
          global[i][j] = -40;
      }
    }   
    System.out.println("--------------------- ");
    Pair l = search(LEVEL);
    double eval = l.score;     
    System.out.println("Best " + eval + ", evaluated nodes " + evaluations); 
  }

  boolean check (int i) { 
    //  return ((i<=n/2 || !visits[i-n/2]) && visits[i]); 
    return visits[i] > 0; 
  }
  
  void swap_rows(int [][] split, int [] sizes, boolean [] seen, int i, int j) {
    for (int k=0;k<sizes[j];k++) 
      seen[split[j][k]] = true;
    if (sizes[i] > sizes[j]) { int t=i; i=j; j=t; }
    for (int k=0;k<sizes[i];k++) {
      int t = split[i][k]; split[i][k] = split[j][k]; split[j][k] = t;
    }
    for (int k=sizes[i];k<sizes[j];k++) {
      split[i][k] = split[j][k];
    }
    int t = sizes[i]; sizes[i] = sizes[j]; sizes[j] = t;
  }  
  private double rollout() {
    for (int j=1;j<N;j++) 
      visits[j] = 1;
    visits[0] = V-1; 
    tour[0] = 0;
    tourSize = 1; // start node already visited
    int node = 0, prev = 0;
    double makespan = 0.0, capacity = 0.0;
    double violations = 0.0;
    double cost = -1.0 * sumservice; 
    while(tourSize < N+V-1) {
      double sum = 0;
      int successors = 0;
      for(int i = 0; i < N; i++) {
        if (check(i) && d[node][i] != 1000000.0) {
          moves[successors++] = i;
          for (int j = 0;j < N;j++) {
            if (i != j) {
              if (check(j)) { 
                if (l[i] > r[j] || makespan + d[node][i] > r[j]) {
                  successors--;
                  break;
                }
              }
            }
          }
        }
      }
      if (successors == 0) {
        for(int i = 0; i < N; i++) 
         if(check(i))   
           moves[successors++] = i;
      }
      for(int i=0; i<successors; i++) {
        value[i] = Math.exp(global[node][moves[i]]);
        sum += value[i];
      }
      double mrand = random.nextDouble()*sum;
      int i=0;
      sum = value[0];
      while(sum<mrand) 
        sum += value[++i];
      prev = node;
      node = moves[i];
      tour[tourSize++] = node;
      visits[node]--;
      cost += d[prev][node]; 
      makespan = Math.max(makespan + d[prev][node],l[node]);  
      capacity += w[node];
      if (node == 0) {
        prev = 0;
        makespan = 0;
        capacity = 0;
      }
      if (capacity > max_capacity) violations += 1.0;
      if (makespan > r[node]) violations += 1.0;
    }
    tour[tourSize++] = 0;  
    cost += d[node][0];   
    makespan = Math.max(makespan + d[node][0], l[0]); 
    if (makespan > r[0]) violations += 1.0;  
    return (100000.0 * violations + cost);
  }   
  
  void adapt(int [] tour_param, int level) {    
    for (int j=1;j<N;j++) 
      visits[j] = 1;
    visits[0] = V-1;
    int successors;
    int node = 0;
    for(int j=1; j<N+V-1; j++)  {
      successors = 0;    
      for(int i = 0; i < N; i++) 
        if (check(i)) 
          moves[successors++] = i;
      double factor = 1.0;
      if (node == 0) factor /= V;
      backup[level][node][tour_param[j]] += factor;     
      double z = 0.0;
      for(int i=0; i<successors; i++) 
      z += Math.exp(global[node][moves[i]]);
      for (int i=0; i<successors; i++) 
         backup[level][node][moves[i]] -= factor *
       Math.exp(global[node][moves[i]])/z;
      node = tour_param[j];
      visits[node]--;
     }
  }  
  void print(int [] tour) {  
    for (int l=1;l<V+N-1;l++) {
      System.out.print("[" + tour[l-1] + "," + tour[l] + "],");
      if (tour[l] == 0)
        System.out.println(); 
    }  
    System.out.println("----------------");
  }
  Pair search(int level) {
    Pair best = new Pair();
    best.score = MAXVALUE;
    if (level == 0) {
      double eval = rollout();  // sets tour and visits
      evaluations++;
      best.score = eval;
      for (int j = 0; j < N+V; j++) 
        best.tour[j] = tour[j]; 
    }
    else {
      for(int i = 0; i < N; i++)  
        for(int j = 0; j < N; j++)  backup[level][i][j] = global[i][j];
      for(int i=0; i<ITERATIONS; i++) {  
        Pair r = search(level - 1); 
        if (r.score < best.score) {
          best.score = r.score;     
          for (int j = 0; j < N+V; j++) 
            best.tour[j] = r.tour[j];
          if(level > SHOW) {
            System.out.println(" Level: " + level + "," + i + ", score: " + r.score + ", runs: " + evaluations);
            if(r.score < 100000.0) print(best.tour);
          }
          adapt(best.tour,level);
        }
        r = null;
      }
      for(int i = 0; i < N; i++) 
        for(int j = 0; j < N; j++)  global[i][j] = backup[level][i][j];
    }
    return best;
  }      
}