source: branches/2457_ExpertSystem/WalkExporter/DirectedWalk.cs @ 17103

using System;
using System.Collections.Generic;
using System.Linq;
using System.Threading;
using HeuristicLab.Analysis.FitnessLandscape;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Encodings.PermutationEncoding;
using HeuristicLab.Problems.Instances;
using HeuristicLab.Problems.Instances.QAPLIB;
using HeuristicLab.Problems.QuadraticAssignment;
using HeuristicLab.Random;
using static HeuristicLab.Analysis.FitnessLandscape.QAPDirectedWalk;

namespace WalkExporter {
  class DirectedWalk {
    public int? Seed;
    public QuadraticAssignmentProblem Problem;
    public WalkType Type;
    public int Paths;
    public double EvaluatedSolutionEquiv;
    public bool BestImprovement;


    public static (Experiment experiment, Dictionary<int, Knowledgebase> training, Dictionary<int, Knowledgebase> test)
        PerformExperiment(WalkType type, bool bestImprovement = true) {

      string typeStr;
      switch (type) {
        case WalkType.RandomRandom: typeStr = "(rr)-dw"; break;
        case WalkType.RandomLocal: typeStr = "(rl)-dw"; break;
        case WalkType.LocalLocal: typeStr = "(ll)-dw"; break;
        case WalkType.LocalInverse: typeStr = "(li)-dw"; break;
        default: throw new ArgumentException(string.Format("unknown type {0}", type), "type");
      }
      var experiment = new Experiment() { Algorithm = typeStr, Trials = new List<Exploration>() };
      var training = new[] { 1, 2, 5, 10, 20, 50, 100, 200, 500 }.ToDictionary(x => x, x => new Knowledgebase() { Problems = new List<ProblemInstanceDescriptor>() });
      var test = new[] { 1, 2, 5, 10, 20, 50, 100, 200, 500 }.ToDictionary(x => x, x => new Knowledgebase() { Problems = new List<ProblemInstanceDescriptor>() });

      foreach (var dimension in new[] { 20, 30, 40 }) {
        var provider = new OneSizeInstanceProvider(dimension);
        foreach (var desc in provider.GetDataDescriptors()) {
          var qapData = provider.LoadData(desc);

          var qap = new QuadraticAssignmentProblem();
          qap.Load(qapData);
          var calculator = new DirectedWalk() {
            Problem = qap,
            BestImprovement = bestImprovement,
            Paths = 1000,
            Seed = 42,
            Type = type
          };
          var solutions = calculator.CalculateRelinkingPoints().ToList();
          var walk = calculator.Run(solutions).ToList();

          var exploration = new Exploration() { Problem = qapData.Name, Dimension = qapData.Dimension, Walks = new List<Walk>() };
          foreach (var w in walk) {
            exploration.Walks.Add(new Walk() {
              QualityTrail = w.Select(x => x.Item2).ToList()
            });
          }
          experiment.Trials.Add(exploration);
          AddToKnowledgeBase(training, qapData, calculator, walk, 0);
          AddToKnowledgeBase(test, qapData, calculator, walk, 500);
        }
      }
      return (experiment, training, test);
    }

    private static void AddToKnowledgeBase(Dictionary<int, Knowledgebase> training, QAPData qapData, DirectedWalk calculator, List<List<Tuple<Permutation, double>>> walk, int start) {
      foreach (var t in training.Keys) {
        var kb = training[t];
        var sbf = CurveAnalysis<Permutation>.GetCharacteristics(walk.Skip(start).Take(t), (left, right) => (1.0 - HammingSimilarityCalculator.CalculateSimilarity(left, right)) * left.Length);
        var trail = walk.Skip(start).Take(t).SelectMany(x => x).Select(x => x.Item2).ToArray();
        var ia = new InformationAnalysis(trail);
        var clen = RuggednessCalculator.CalculateCorrelationLength(trail, out double[] acf);

        var pid = new ProblemInstanceDescriptor() {
          Name = qapData.Name,
          Class = Util.GetGeneratorClass(qapData.Name),
          Dimension = qapData.Dimension,
          DescriptionEffort = calculator.EvaluatedSolutionEquiv / 1000 * t,
          Features = new List<KeyValue>()
        };
        pid.Features.Add(new KeyValue { Key = "AC1", ContinuousValue = acf[1] });
        pid.Features.Add(new KeyValue { Key = "CorrelationLength", DiscreteValue = clen });
        foreach (var f in CurveAnalysisResult.AllFeatures.Zip(sbf.GetValues(), (a, b) => new { Feature = a, Value = b })) {
          pid.Features.Add(new KeyValue { Key = f.Feature.ToString(), ContinuousValue = f.Value });
        }
        foreach (var f in ia.GetFeatures()) {
          if (f.Item2 is double d) pid.Features.Add(new KeyValue { Key = f.Item1, ContinuousValue = d });
          else if (f.Item2 is int i) pid.Features.Add(new KeyValue { Key = f.Item1, DiscreteValue = i });
        }

        kb.Problems.Add(pid);
      }
    }

    public IEnumerable<Permutation> CalculateRelinkingPoints() {
      IRandom random = Seed.HasValue ? new MersenneTwister((uint)Seed.Value) : new MersenneTwister();
      var pathCount = Paths;

      var perm = new Permutation(PermutationTypes.Absolute, Problem.Weights.Rows, random);
      var startLocal = Type == WalkType.LocalGlobal || Type == WalkType.LocalLocal;
      var targetLocal = Type == WalkType.LocalLocal || Type == WalkType.RandomLocal;
      var targetGlobal = Type == WalkType.LocalGlobal || Type == WalkType.RandomGlobal;
      var inverseWalk = Type == WalkType.LocalInverse;

      if (Type == WalkType.LocalInverse) pathCount--; // inverse walks equal one walk per solution
      for (var i = 0; i <= pathCount; i++) {
        var start = i % 2 == 0;
        var target = i % 2 == 1;
        if (inverseWalk || start && startLocal) {
          perm = new Permutation(PermutationTypes.Absolute, Problem.Weights.Rows, random);
          var fit = new DoubleValue(QAPEvaluator.Apply(perm, Problem.Weights, Problem.Distances));
          EvaluatedSolutionEquiv++;
          var evalSol = new IntValue(0);
          QAPExhaustiveSwap2LocalImprovement.ImproveFast(perm, Problem.Weights, Problem.Distances, fit, new IntValue(0), evalSol, Problem.Maximization.Value, int.MaxValue, CancellationToken.None);
          EvaluatedSolutionEquiv += evalSol.Value;
        } else if (target && (targetLocal || targetGlobal)) {
          if (targetGlobal) {
            perm = Problem.BestKnownSolutions.SampleRandom(random);
          } else {
            perm = new Permutation(PermutationTypes.Absolute, Problem.Weights.Rows, random);
            var fit = new DoubleValue(QAPEvaluator.Apply(perm, Problem.Weights, Problem.Distances));
            EvaluatedSolutionEquiv++;
            var evalSol = new IntValue(0);
            QAPExhaustiveSwap2LocalImprovement.ImproveFast(perm, Problem.Weights, Problem.Distances, fit, new IntValue(0), evalSol, Problem.Maximization.Value, int.MaxValue, CancellationToken.None);
            EvaluatedSolutionEquiv += evalSol.Value;
          }
        } else { // random
          BiasedShuffle(perm, random);
        }
        yield return (Permutation)perm.Clone();
      }
    }

    public IEnumerable<List<Tuple<Permutation, double>>> Run(IEnumerable<Permutation> permutations) {
      if (Type == WalkType.LocalInverse) return RunInverse(permutations);
      return RunRegular(permutations);
    }

    private IEnumerable<List<Tuple<Permutation, double>>> RunInverse(IEnumerable<Permutation> permutations) {
      var qap = (QuadraticAssignmentProblem)Problem;
      var min = qap.LowerBound.Value;
      var max = qap.AverageQuality.Value;
      var bestImprovement = BestImprovement;

      foreach (var start in permutations) {
        var startFitness = QAPEvaluator.Apply(start, qap.Weights, qap.Distances);
        EvaluatedSolutionEquiv++;

        // inverse walks are applied to all solutions
        Func<Tuple<Permutation, double>, IEnumerable<Tuple<Permutation, double>>> inverseNeighborFunc = (p) => RestrictedInverseNeighborhood(qap, p.Item1, p.Item2, start);
        var walk = DirectedWalk<Permutation>.WalkRestricted(qap.Maximization.Value, inverseNeighborFunc, start, startFitness, !bestImprovement);
        yield return walk.Select(x => Tuple.Create(x.Item1, (x.Item2 - min) / (max - min))).ToList();
      } // end paths
    }

    private IEnumerable<List<Tuple<Permutation, double>>> RunRegular(IEnumerable<Permutation> permutations) {
      var iter = permutations.GetEnumerator();
      if (!iter.MoveNext()) yield break;
      var bestImprovement = BestImprovement;

      var qap = (QuadraticAssignmentProblem)Problem;
      var min = qap.LowerBound.Value;
      var max = qap.AverageQuality.Value;

      var start = iter.Current;

      while (iter.MoveNext()) {
        var startFitness = QAPEvaluator.Apply(start, qap.Weights, qap.Distances);
        EvaluatedSolutionEquiv++;
        var end = iter.Current;

        var invSol = new int[start.Length];
        Func<Tuple<Permutation, double>, IEnumerable<Tuple<Permutation, double>>> regularNeighborFunc = (p) => RestrictedNeighborhood(qap, p.Item1, p.Item2, end, invSol);
        var walk = DirectedWalk<Permutation>.WalkRestricted(qap.Maximization.Value, regularNeighborFunc, start, startFitness, !bestImprovement);

        var trail = new List<Tuple<Permutation, double>>();
        foreach (var step in walk) {
          for (var i = 0; i < invSol.Length; i++)
            invSol[step.Item1[i]] = i;
          trail.Add(Tuple.Create(step.Item1, (step.Item2 - min) / (max - min)));
        }
        yield return trail;

        start = end;
      } // end paths
    }

    private IEnumerable<Tuple<Permutation, double>> RestrictedInverseNeighborhood(QuadraticAssignmentProblem qap, Permutation current, double currentFit, Permutation start) {
      var evalSolPerMove = 4.0 / current.Length;
      var N = current.Length;
      for (var index = 0; index < N; index++) {
        if (current[index] != start[index]) continue;
        for (var k = 0; k < N; k++) {
          if (k == index) continue;
          var fit = QAPSwap2MoveEvaluator.Apply(current, new Swap2Move(index, k), qap.Weights, qap.Distances);
          EvaluatedSolutionEquiv += evalSolPerMove;
          current.Swap(index, k);
          yield return Tuple.Create(current, currentFit + fit);
          current.Swap(index, k); // reverse
        }
      }
    }

    private IEnumerable<Tuple<Permutation, double>> RestrictedNeighborhood(QuadraticAssignmentProblem qap, Permutation current, double currentFit, Permutation target, int[] inverse) {
      var evalSolPerMove = 4.0 / current.Length;
      for (var index = 0; index < current.Length; index++) {
        if (current[index] == target[index]) continue;
        var k = inverse[target[index]];
        var fit = QAPSwap2MoveEvaluator.Apply(current, new Swap2Move(index, k), qap.Weights, qap.Distances);
        EvaluatedSolutionEquiv += evalSolPerMove;
        current.Swap(index, k);
        yield return Tuple.Create(current, currentFit + fit);
        current.Swap(index, k); // reverse
      }
    }

    // permutation must be strictly different in every position
    private static void BiasedShuffle(Permutation p, IRandom random) {
      for (var i = p.Length - 1; i > 0; i--) {
        // Swap element "i" with a random earlier element (excluding itself)
        var swapIndex = random.Next(i);
        var h = p[swapIndex];
        p[swapIndex] = p[i];
        p[i] = h;
      }
    }
  }
}