source: branches/PTSP/HeuristicLab.Problems.PTSP/3.3/EstimatedPTSP.cs @ 12269

using System.Text;
using System.Threading.Tasks;
using HeuristicLab.Optimization;
using HeuristicLab.PluginInfrastructure;
using HeuristicLab.Core;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Problems.Instances;
using HeuristicLab.Encodings.PermutationEncoding;
using HeuristicLab.Common;
using HeuristicLab.Parameters;
using HeuristicLab.Data;
using HeuristicLab.Random;
using System;
using System.Linq;

namespace HeuristicLab.Problems.PTSP {
  [Item("Estimated Probabilistic Traveling Salesman Problem", "Represents an estimated Probabilistic Traveling Salesman Problem.")]
  [Creatable("Problems")]
  [StorableClass]
  public sealed class EstimatedProbabilisticTravelingSalesmanProblem : ProbabilisticTravelingSalesmanProblem, IStorableContent,
  IProblemInstanceConsumer<TSPData> {

    private ItemList<ItemList<IntValue>> realizations;

    [StorableConstructor]
    private EstimatedProbabilisticTravelingSalesmanProblem(bool deserializing) : base(deserializing) { }
    private EstimatedProbabilisticTravelingSalesmanProblem(EstimatedProbabilisticTravelingSalesmanProblem original, Cloner cloner)
      : base(original, cloner) {
    }

    public override IDeepCloneable Clone(Cloner cloner) {
      return new EstimatedProbabilisticTravelingSalesmanProblem(this, cloner);
    }

    public override double Evaluate(Individual individual, IRandom random) {
      Permutation p = individual.Permutation();
      // Estimation-based evaluation
      MersenneTwister r = new MersenneTwister();
      double estimatedSum = 0;
      for (int i = 0; i < realizations.Count; i++) {
        int singleRealization = -1, firstNode = -1;
        for (int j = 0; j < realizations[i].Count; j++) {
          if (realizations[i][p[j]].Value == 1) {
            if (singleRealization != -1) {
              estimatedSum += DistanceMatrix[singleRealization, p[j]];
            } else {
              firstNode = p[j];
            }
            singleRealization = p[j];
          }
        }
        if (singleRealization != -1) {
          estimatedSum += DistanceMatrix[singleRealization, firstNode];
        }
      }
      return estimatedSum / SampleSize.Value;
    }

    public double[] EvaluateWithParams(DistanceMatrix distances, DoubleArray probabilities, ItemList<ItemList<IntValue>> realizations, Permutation individual ) {
      // Estimation-based evaluation
      MersenneTwister r = new MersenneTwister();
      double estimatedSum = 0;
      double[] partialSums = new double[realizations.Count];
      for (int i = 0; i < realizations.Count; i++) {
        partialSums[i] = 0;
        int singleRealization = -1, firstNode = -1;
        for (int j = 0; j < realizations[i].Count; j++) {
          if (realizations[i][individual[j]].Value == 1) {
            if (singleRealization != -1) {
              partialSums[i] += distances[singleRealization, individual[j]];
            } else {
              firstNode = individual[j];
            }
            singleRealization = individual[j];
          }
        }
        if (singleRealization != -1) {
          partialSums[i] += distances[singleRealization, firstNode];
        }
        estimatedSum += partialSums[i];
      }
      double mean = estimatedSum / realizations.Count;
      double variance = 0;
      for (int i = 0; i < realizations.Count; i++) {
        variance += Math.Pow((partialSums[i] - mean), 2);
      }
      variance = variance / realizations.Count;
      return new double[] {mean,variance};
    }

    public EstimatedProbabilisticTravelingSalesmanProblem() {
      Parameters.Add(new ValueParameter<IntValue>("SampleSize", "Size of the sample for the estimation-based evaluation"));
      Operators.Add(new PTSPEstimatedInversionMovePathEvaluator());
      Operators.Add(new PTSPEstimatedInsertionEvaluator());
      Operators.Add(new PTSPExhaustiveInversionLocalImprovement());
      Operators.Add(new PTSPExhaustiveInsertionLocalImprovement());
      Encoding.ConfigureOperators(Operators.OfType<IOperator>());
    }

    private int Ttest(int ProblemSize) {
      MersenneTwister random = new MersenneTwister();
      Permutation p1 = new Permutation(PermutationTypes.RelativeUndirected, ProblemSize, random);
      Permutation p2 = new Permutation(PermutationTypes.RelativeUndirected, ProblemSize, random);
      ItemList<ItemList<IntValue>> realizations = new ItemList<ItemList<IntValue>>();
      int index = -1;
      while (true) {
        for (int i = index+1; i < index+5; i++) {
          realizations.Add(new ItemList<IntValue>());
          for (int j = 0; j < ProblemSize; j++) {
            if (ProbabilityMatrix[j] > random.NextDouble()) {
              realizations.ElementAt(i).Add(new IntValue(1));
            } else {
              realizations.ElementAt(i).Add(new IntValue(0));
            }
          }
        }
        index += 4;
        double[] eval1 = EvaluateWithParams(DistanceMatrix, ProbabilityMatrix, realizations, p1);
        double[] eval2 = EvaluateWithParams(DistanceMatrix, ProbabilityMatrix, realizations, p2);
        double sx1x2 = Math.Sqrt((eval1[1]+eval2[1])/2);
        int degrees = 2 * realizations.Count - 2;
        double t = (eval1[0]-eval2[0])/(sx1x2*Math.Sqrt(2.0/(double)realizations.Count));
      }
    }

    public override void Load(TSPData data) {
      base.Load(data);
      // For now uses sample size of 20 but should use Student's t-test
      //Ttest(data.Dimension);
      SampleSize = new IntValue(100);
      MersenneTwister r = new MersenneTwister();
      int countOnes = 0;
      realizations = new ItemList<ItemList<IntValue>>(SampleSize.Value);
      for (int i = 0; i < SampleSize.Value; i++) {
        ItemList<IntValue> newRealization = new ItemList<IntValue>();
        while (countOnes < 4) { //only generate realizations with at least 4 cities visited
          countOnes = 0;
          newRealization.Clear();
          for (int j = 0; j < data.Dimension; j++) {
            if (ProbabilityMatrix[j] > r.NextDouble()) {
              newRealization.Add(new IntValue(1));
              countOnes++;
            } else {
              newRealization.Add(new IntValue(0));
            }
          }
        }
        countOnes = 0;
        realizations.Add(newRealization);
      }

      foreach (var op in Operators.OfType<PTSPPathMoveEvaluator>()) {
        op.RealizationsParameter.Value = realizations;
      }
      foreach (var op in Operators.OfType<PTSPExhaustiveInversionLocalImprovement>()) {
        op.RealizationsParameter.Value = realizations;
      }
    }
  }
}