source: branches/HeuristicLab.Problems.Orienteering/HeuristicLab.Problems.Orienteering/3.3/OrienteeringScript.cs @ 11270

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2014 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
 *
 * This file is part of HeuristicLab.
 *
 * HeuristicLab is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * HeuristicLab is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
 */
#endregion

using System;
using System.Collections.Generic;
using System.Globalization;
using System.IO;
using System.Linq;
using System.Threading;
using HeuristicLab.Analysis;
using HeuristicLab.Core;
using HeuristicLab.Optimization;
using HeuristicLab.Random;
using HeuristicLab.Scripting;

namespace HeuristicLab.Problems.Orienteering {

  /* 
   * Based on
   * Implementation of an Pareto Ant Colony
   * Optimization (PACO) metaheuristic and a
   * Pareto Variable Neighborhood Search (PVNS)
   * metaheuristic for the Multi-Objective
   * Orienteering Problem (MOP)
   * Copyright: © 2007 by Michael Schilde
   * 
   * See Metaheuristics for the bi-objective orienteering problem
   * http://link.springer.com/article/10.1007/s11721-009-0029-5
   * Michael Schilde, Karl F. Doerner, Richard F. Hartl, Guenter Kiechle
   */
  public class OrienteeringScript : CSharpScriptBase {
    #region Constants

    private const bool ShowStatus = true;

    private const double PrecisionUsed = 0.0000001;

    private const string VersionRequired = "*tpv1.1* or *tpv1.0*";

    private const int ConstraintTimeDistance = -1;

    private static readonly TimeSpan RunDuration = TimeSpan.FromSeconds(15);

    public enum ErrorCode {
      ErrorNone = 0,
      ErrorOpenInputFile = 10,
      ErrorOpenOutputFile = 11,
      ErrorConvertToLower = 12,
      ErrorRemovingWhitespace = 13,
      ErrorWrongVersion = 14,
      ErrorReadingData = 20,
      ErrorInvalidDefinitionN = 30,
      ErrorInvalidDefinitionC = 31,
      ErrorInvalidDefinitionU = 32,
      ErrorInvalidDefinitionL = 33,
      ErrorInvalidDefinitionX = 34,
      ErrorInvalidDefinitionI = 35,
      ErrorInvalidDefinitionP = 36,
      ErrorInvalidDefinitionS = 37,
      ErrorInvalidDefinitionB = 38,
      ErrorInvalidDefinitionE = 39,
      ErrorInvalidDefinitionD = 40,
      ErrorInvalidDefinitionM = 41,
      ErrorMissingStartingPoint = 50,
      ErrorMissingTerminusPoint = 51,
      ErrorTooFewConstraints = 63,
      ErrorTooFewPoints = 64,
      ErrorTooFewScores = 65,
      ErrorTooManyConstraints = 66,
      ErrorTooManyPoints = 67,
      ErrorTooManyScores = 68
    }

    #endregion

    #region Fields
    private TimeSpan avgTimeDurationVns;
    private List<TimeSpan> timeDurationVns = new List<TimeSpan>();
    private List<int> solutionsFoundVns = new List<int>();
    private List<List<TourData>> TourBackupVNS = new List<List<TourData>>();
    private List<int> solutionsFoundPrVns = new List<int>();
    private List<ProblemData> problems = new List<ProblemData>();
    private double averageScoreVns;
    private double bestScoreVns;
    private TimeSpan timeOfBestScoreVns;
    private double worstScoreVns;

    private ItemList<IItem> results;
    private DataTable scores;

    #endregion

    public override void Main() {
      Thread.CurrentThread.CurrentCulture = new CultureInfo("en-US");

      vars.Clear();

      const int numberOfRuns = 10;
      const string inputFilename = @"C:\Users\P41107\Desktop\moop\1 objective\1_p66\1_p66_t080.txt";
      const string csvFilename = @"C:\temp\results1.csv";
      const string summaryFilename = @"C:\temp\summary.csv";
      const int numberOfIterations = 10000;
      const int acceptFrequency = 100;
      const double tresholdValue = 0.95;
      const double fixedVisitingTime = 0.0;
      const double fixedPenalty = fixedVisitingTime;
      const double maximumDistanceOverride = 0.0;
      const bool performPathRelinkingComparison = true;

      string testInstanceName = inputFilename.Contains(@"\")
        ? inputFilename.Substring(inputFilename.LastIndexOf(@"\") + 1)
        : inputFilename.Substring(inputFilename.LastIndexOf(@"/") + 1);

      bool csvExists = File.Exists(summaryFilename);

      // Open the results CSV file for writing the data
      using (var csvFile = new StreamWriter(new FileStream(csvFilename, FileMode.Create)))
      // Open the summary CSV file for appending the data
      using (var summaryFile = new StreamWriter(new FileStream(summaryFilename, FileMode.OpenOrCreate))) {

        // Remember the overall starting time
        var startOverall = DateTime.Now.Ticks;

        for (int runCounter = 0; runCounter < numberOfRuns; runCounter++) {
          // Print out a status information
          if (ShowStatus)
            Console.WriteLine("Run {0} started.", (runCounter + 1));

          // Create a new problem data container
          var problem = new ProblemData {
            FixedPenalty = fixedVisitingTime,
            Constraints = new List<ConstraintData>(),
            Points = new List<PointData>(),
            Edges = new List<List<EdgeData>>(),
            MTR = new MersenneTwister(),
            VNS = new VnsData {
              NumberOfIterations = numberOfIterations,
              AcceptFrequency = acceptFrequency,
              TresholdValue = tresholdValue,
              InitialTour = new TourData { Path = new List<int>(), Scores = new List<int>() },
              DominantTours = new List<TourData>()
            },
          };

          // Do the initialization
          ErrorCode errorCode = LoadData(inputFilename, problem);
          if (errorCode != 0) {
            Console.WriteLine("Initialization failed. (Code: {0})", errorCode);
            Console.WriteLine(GetErrorDescription(errorCode));
            Console.WriteLine("Exiting.");
            return;
          }

          // Override the distance limit from file by parameter if specified
          if (maximumDistanceOverride > 0.0) {
            problem.Constraints.Clear();
            problem.Constraints.Add(new ConstraintData {
              Type = Constraint.Upperbound,
              Target = -1,
              TargetValue = maximumDistanceOverride
            });
          }

          results = new ItemList<IItem>();
          vars["run" + runCounter] = results;
          scores = new DataTable("Scores");
          results.Add(scores);
          scores.Rows.AddRange(Enumerable.Range(0, problem.ScoreCount).Select(s => new DataRow("Score " + s)));

          if (ShowStatus)
            Console.WriteLine("\tVNS started.");

          // Run the VNS algorithm
          var startPart = DateTime.Now.Ticks;
          var durationPart = RunDuration.Ticks;
          RunVns(problem, startPart, durationPart);
          var endPart = DateTime.Now.Ticks;
          var temp = CalculateDifference(startPart, endPart);
          avgTimeDurationVns += temp;
          timeDurationVns.Add(temp);
          solutionsFoundVns.Add((int)problem.VNS.DominantTours.Count);

          if (ShowStatus)
            Console.WriteLine("\tVNS completed ({0} Tours in {1:g}).", problem.VNS.DominantTours.Count, temp);

          if (performPathRelinkingComparison) {
            // Backup original VNS results
            TourBackupVNS.Add(problem.VNS.DominantTours);
          }

          // If we have two objectives, run path relinking
          if (problem.ScoreCount == 2) {
            // Run the path relinking for the VNS solutions
            if (ShowStatus)
              Console.WriteLine("\tPath relinking for VNS solutions started.");
            int toursFoundPrVns = PathRelinking(problem, problem.VNS.DominantTours, startPart);

            temp = CalculateDifference(startPart, DateTime.Now.Ticks) - temp;
            solutionsFoundPrVns.Add(toursFoundPrVns);
            if (ShowStatus)
              Console.WriteLine("\tPath relinking for VNS solutions completed ({0} tours in {1:g}).", toursFoundPrVns, temp);
          }

          // Store all the information in the run-spanning storage
          problems.Add(problem);
        }

        // If we have more than one objective, calculate metrics
        if (problems[0].ScoreCount == 2) {
          if (ShowStatus)
            Console.WriteLine("Performing 2-objective post-calculations.");
          throw new NotImplementedException();
        } else {
          if (ShowStatus)
            Console.WriteLine("\tPerforming 1-objective post-calculations.");

          // Only one objective value, so just calculate best, average and worst score
          for (int counterRuns = 0; counterRuns < numberOfRuns; counterRuns++) {
            averageScoreVns += problems[counterRuns].VNS.DominantTours[0].Scores[0];

            if (problems[counterRuns].VNS.DominantTours[0].Scores[0] > bestScoreVns) {
              bestScoreVns = problems[counterRuns].VNS.DominantTours[0].Scores[0];
              timeOfBestScoreVns = problems[counterRuns].VNS.DominantTours[0].TimeOfAdding;
            }
            if ((counterRuns == 0) || (problems[counterRuns].VNS.DominantTours[0].Scores[0] < worstScoreVns)) {
              worstScoreVns = problems[counterRuns].VNS.DominantTours[0].Scores[0];
            }
          }

          // Calculate the average computation times
          averageScoreVns /= numberOfRuns;

          if (ShowStatus)
            Console.WriteLine("\tWriting output files.");

          // Write out the CSV results file
          csvFile.WriteLine("Test instance;Tmax;# runs;iter. VNS;accept frequency;treshold");
          csvFile.Write(testInstanceName + ";");
          csvFile.Write("{0};", problems[0].Constraints[0].TargetValue);
          csvFile.Write("{0};", numberOfRuns);
          csvFile.Write("{0};", numberOfIterations);
          csvFile.Write("{0};", acceptFrequency);
          csvFile.Write("{0};", tresholdValue);
          csvFile.WriteLine();
          csvFile.WriteLine();

          // If we have more than one objective, print out the full detail report
          if (problems[0].ScoreCount > 1) {
            throw new NotImplementedException();
          } else {
            // Only one objective, so print out the short report only
            csvFile.WriteLine("Run;CPU VNS;Score VNS");
            for (int counterRuns = 0; counterRuns < numberOfRuns; counterRuns++) {
              csvFile.Write("{0};", counterRuns + 1); // Number of runs
              csvFile.Write("{0};", timeDurationVns[counterRuns]); // CPU VNS
              csvFile.Write("{0};", problems[counterRuns].VNS.DominantTours[0].Scores[0]); // Score VNS
              //csv_file << convert_locale(time_of_best_score_aco) << ";";          // Runtime until finding best tour ACO
              //csv_file << convert_locale(time_of_best_score_vns) << ";";          // Runtime until finding best tour VNS
              csvFile.WriteLine();
            }
            csvFile.WriteLine();
            for (int counterRuns = 0; counterRuns < numberOfRuns; counterRuns++) {
              csvFile.WriteLine("Run;{0}", counterRuns + 1);
              csvFile.WriteLine(";;Length;Time Of Adding;Score;Path");
              csvFile.Write(";;");


              csvFile.WriteLine();
              csvFile.WriteLine();

              csvFile.WriteLine(";Tour VNS");
              csvFile.WriteLine(";;Length;Time Of Adding;Score;Path");
              csvFile.Write(";;");
              csvFile.Write("{0};", problems[counterRuns].VNS.DominantTours[0].Length);
              csvFile.Write("{0};", problems[counterRuns].VNS.DominantTours[0].TimeOfAdding);
              csvFile.Write("{0};", problems[counterRuns].VNS.DominantTours[0].Scores[0]);
              csvFile.Write("{0}", problems[counterRuns].VNS.DominantTours[0].Path[0]);
              for (int counterPoints = 1;
                counterPoints < problems[counterRuns].VNS.DominantTours[0].PathSize;
                counterPoints++) {
                csvFile.Write("-{0}", problems[counterRuns].VNS.DominantTours[0].Path[counterPoints]);
              }
              csvFile.WriteLine();
              csvFile.WriteLine();
            }
            // Append the CSV summary file information
            if (!csvExists) {
              summaryFile.Write("Test instance;Tmax;# runs;iter. VNS;accept frequency;treshold;");
              summaryFile.Write("avg. CPU VNS;max score VNS;avg. score VNS;min score VNS;time of adding best tour VNS");
              summaryFile.WriteLine();
            }
            summaryFile.Write("{0};", testInstanceName);
            summaryFile.Write("{0};", problems[0].Constraints[0].TargetValue);
            summaryFile.Write("{0};", numberOfRuns);
            summaryFile.Write("{0};", numberOfIterations); // number of iterations vns
            summaryFile.Write("{0};", acceptFrequency); // accept frequency
            summaryFile.Write("{0};", tresholdValue); // treshold value
            summaryFile.Write("{0};", avgTimeDurationVns); // CPU VNS
            summaryFile.Write("{0};", bestScoreVns); // Best score VNS
            summaryFile.Write("{0};", averageScoreVns); // Average score VNS
            summaryFile.Write("{0};", worstScoreVns); // Worst score VNS
            summaryFile.Write("{0};", timeOfBestScoreVns); // Time until best score VNS
            summaryFile.WriteLine();
          }
        }
      }
    }

    #region Initialization

    private static ErrorCode LoadData(string fileName, ProblemData problem) {
      StreamReader inputFile;
      try {
        // Open the input file
        inputFile = new StreamReader(new FileStream(fileName, FileMode.OpenOrCreate));
      } catch (IOException) {
        return ErrorCode.ErrorOpenInputFile;
      }

      // Initialize the Data structure
      problem.ConstraintCount = 0;
      problem.PointCount = 0;
      problem.ScoreCount = 0;
      problem.StartingPoint = -1;
      problem.TerminusPoint = -1;
      problem.MatrixGiven = false;

      ErrorCode error = ErrorCode.ErrorNone;

      int lineNumber = 0;

      // Read the lines from the file
      while ((!inputFile.EndOfStream) && (error == ErrorCode.ErrorNone)) {
        string lineBuffer;
        try {
          lineBuffer = inputFile.ReadLine();
        } catch (IOException) {
          return ErrorCode.ErrorReadingData;
        }

        // Remove all whitespace
        lineBuffer = new string(lineBuffer.ToCharArray().Where(c => !char.IsWhiteSpace(c)).ToArray());

        // If the line contains data, continue
        if (lineBuffer.Length > 1) {
          // If the line contains not only a comment, process it
          int commentPosition = lineBuffer.IndexOf("//");
          if (commentPosition != 0) {
            // Count the lines that include data
            lineNumber++;

            // Remove a comment if there is one
            if (commentPosition > 0)
              lineBuffer = lineBuffer.Remove(commentPosition);

            // Convert the line to lowercase characters
            lineBuffer = lineBuffer.ToLower();

            // Make sure the first data line includes the correct version
            if (lineNumber == 1) {
              string versions = VersionRequired;
              bool versionFound = false;
              do {
                if (lineBuffer == versions.Substring(versions.IndexOf(" or ") + 1)) {
                  versionFound = true;
                  break;
                }
                versions = versions.Remove(versions.IndexOf(" or "));
              }
              while (versions.IndexOf(" or ") != -1);

              //if (line_buffer != POPT_Version_Required) return ERROR_WRONG_VERSION;
            } else {
              // Extract the needed information
              error = ExtractData(lineBuffer, problem);
            }
          }
        }
      }

      // Check the data for inconsistencies in the data structure
      if (error == ErrorCode.ErrorNone) {
        if (problem.Constraints.Count < problem.ConstraintCount)
          return ErrorCode.ErrorTooFewConstraints;

        if (problem.Points.Count < problem.PointCount)
          return ErrorCode.ErrorTooFewPoints;

        if (problem.StartingPoint == -1)
          return ErrorCode.ErrorMissingStartingPoint;

        if (problem.TerminusPoint == -1)
          return ErrorCode.ErrorMissingTerminusPoint;

        // Calculate all distances from one point to another 
        // and initialize pheromone values
        int pointCounter1;
        int pointCounter2;
        double maxDistance;
        if (!problem.MatrixGiven) {
          maxDistance = 0.0;
          for (pointCounter1 = (problem.PointCount - 1); pointCounter1 >= 0; --pointCounter1) {
            // Delete the last generated edgeset
            var tempEdgeset = new List<EdgeData>();

            for (pointCounter2 = (problem.PointCount - 1); pointCounter2 >= 0; --pointCounter2) {
              // Calculate the distances
              double distanceX = Math.Abs(problem.Points[pointCounter1].PositionX - problem.Points[pointCounter2].PositionX);
              double distanceY = Math.Abs(problem.Points[pointCounter1].PositionY - problem.Points[pointCounter2].PositionY);
              var tempEdge = new EdgeData {
                Distance = Math.Sqrt(Math.Pow(distanceX, 2.0) + Math.Pow(distanceY, 2.0))
              };

              // Remember the longest distance
              if ((tempEdge.Distance - maxDistance) >= PrecisionUsed) maxDistance = tempEdge.Distance; //(temp_edge.Distance > max_distance)

              // Initialize the edge's pheromone values
              //temp_edge.Pheromone.assign(Problem.Score_Count, Problem.ACO.Initial_Pheromone_Value);

              // Put the now generated edge to the temporal edge set
              tempEdgeset.Insert(0, tempEdge);
            }

            // Put the now generated temporal edge set to the original edge set
            problem.Edges.Insert(0, tempEdgeset);
          }
        } else {
          maxDistance = 0.0;

          // Find the maximum distance in the problem set
          for (pointCounter1 = 0; pointCounter1 < problem.PointCount; pointCounter1++) {
            for (pointCounter2 = 0; pointCounter2 < problem.PointCount; pointCounter2++) {
              if ((problem.Edges[pointCounter1][pointCounter2].Distance - maxDistance) >= PrecisionUsed) { //(Problem.Edges[point_counter1][point_counter2].Distance > max_distance)
                maxDistance = problem.Edges[pointCounter1][pointCounter2].Distance;
              }
            }
          }
        }
      }

      // Close the input file
      inputFile.Close();

      // Return with no error
      return error;
    }

    private static ErrorCode ExtractData(string inputString, ProblemData problem) {
      string temp;

      switch (inputString[0]) {
        case 'p': // Point description
          if (inputString.Length < 8)
            return ErrorCode.ErrorInvalidDefinitionP;

          // Remove the 'p,'
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
          // Get the position_x
          string positionX = inputString.Substring(0, inputString.IndexOf(','));
          // Remove the position_x
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
          // Get the postion_y
          string positionY = inputString.Substring(0, inputString.IndexOf(','));
          // Remove the position_y to get the scores
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);

          var tempData = new PointData {
            PositionX = double.Parse(positionX),
            PositionY = double.Parse(positionY),
            Score = new List<int>()
          };

          // Extract all score values for this point
          while (true) {
            string score;
            if (inputString.IndexOf(',') != -1) {
              score = inputString.Substring(0, inputString.IndexOf(','));
              inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
              tempData.Score.Add(int.Parse(score));
            } else {
              score = inputString;
              tempData.Score.Add(int.Parse(score));
              break;
            }
          }
          problem.Points.Add(tempData);
          break;

        case 'd': // Distance information
          if (inputString.Length < 11)
            return ErrorCode.ErrorInvalidDefinitionD;

          var tempEdges = new List<EdgeData>();

          // Remove the 'd,'
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);

          // Extract all distances for this point
          while (true) {
            string distance;
            if (inputString.IndexOf(',') != -1) {
              distance = inputString.Substring(0, inputString.IndexOf(','));
              inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
              var tempEdge = new EdgeData {
                Distance = double.Parse(distance)
              };
              tempEdges.Add(tempEdge);
            } else {
              distance = inputString;
              var tempEdge = new EdgeData {
                Distance = double.Parse(distance)
              };
              tempEdges.Add(tempEdge);
              break;
            }
          }
          problem.Edges.Add(tempEdges);
          break;

        case 'n': // Number of points
          if (inputString.Length < 3)
            return ErrorCode.ErrorInvalidDefinitionN;

          // Remove the 'n,'
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
          // Extract the number of points
          problem.PointCount = int.Parse(inputString);
          break;

        case 'm': // 1 if Distance-Matrix is given
          if (inputString.Length < 3)
            return ErrorCode.ErrorInvalidDefinitionM;

          // Remove the 'm,'
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
          // Extract the number of points
          if (int.Parse(inputString) != 0) problem.MatrixGiven = true;
          break;

        case 'c': // Number of constraints
          if (inputString.Length < 3)
            return ErrorCode.ErrorInvalidDefinitionC;

          // Remove the 'c,'
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
          problem.ConstraintCount = int.Parse(inputString);
          break;

        case 's': // Number of scores per point
          if (inputString.Length < 3)
            return ErrorCode.ErrorInvalidDefinitionS;

          // Remove the 's,'
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
          problem.ScoreCount = int.Parse(inputString);
          break;

        case 'b': // Index of the starting point
          if (inputString.Length < 3)
            return ErrorCode.ErrorInvalidDefinitionB;

          // Remove the 'b,'
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
          problem.StartingPoint = int.Parse(inputString);
          break;

        case 'e': // Index of the terminus point
          if (inputString.Length < 3)
            return ErrorCode.ErrorInvalidDefinitionE;

          // Remove the 'e,'
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
          problem.TerminusPoint = int.Parse(inputString);
          break;

        case 'u': // Upper bound constraint
          if (inputString.Length < 5)
            return ErrorCode.ErrorInvalidDefinitionU;

          // Remove the 'u,'
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
          // Get the target
          temp = inputString.Substring(0, inputString.IndexOf(','));
          // Remove the target to get the target value
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);

          problem.Constraints.Add(new ConstraintData {
            Type = Constraint.Upperbound,
            Target = int.Parse(temp),
            TargetValue = double.Parse(inputString)
          });
          break;

        case 'l': // Lower bound constraint
          if (inputString.Length < 5)
            return ErrorCode.ErrorInvalidDefinitionL;

          // Remove the 'l,'
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
          // Get the target
          temp = inputString.Substring(0, inputString.IndexOf(','));
          // Remove the target to get the target value
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);

          problem.Constraints.Add(new ConstraintData {
            Type = Constraint.Lowerbound,
            Target = int.Parse(temp),
            TargetValue = double.Parse(inputString)
          });
          break;

        case 'x': // Maximization constraint
          if (inputString.Length < 5)
            return ErrorCode.ErrorInvalidDefinitionX;

          // Remove the 'x,'
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
          // Get the target
          temp = inputString.Substring(0, inputString.IndexOf(','));
          // Remove the target to get the target value
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);

          problem.Constraints.Add(new ConstraintData {
            Type = Constraint.Maximize,
            Target = int.Parse(temp),
            TargetValue = Double.Parse(inputString)
          });
          break;

        case 'i': // Minimization constraint
          if (inputString.Length < 5)
            return ErrorCode.ErrorInvalidDefinitionI;

          // Remove the 'i,'
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);
          // Get the target
          temp = inputString.Substring(0, inputString.IndexOf(','));
          // Remove the target to get the target value
          inputString = inputString.Remove(0, inputString.IndexOf(',') + 1);

          problem.Constraints.Add(new ConstraintData {
            Type = Constraint.Minimize,
            Target = int.Parse(temp),
            TargetValue = double.Parse(inputString)
          });
          break;
      }
      return ErrorCode.ErrorNone;
    }

    private static string GetErrorDescription(ErrorCode errorCode) {
      switch (errorCode) {
        case ErrorCode.ErrorOpenInputFile: return "The input file could not be opened.";
        case ErrorCode.ErrorOpenOutputFile: return "The output file could not be opened.";
        case ErrorCode.ErrorConvertToLower: return "An error occoured at the conversion to lowercase.";
        case ErrorCode.ErrorRemovingWhitespace: return "An error occoured at the removal of whitespace.";
        case ErrorCode.ErrorWrongVersion: return ("The input file has the wrong version. Version " + VersionRequired + " is required.");
        case ErrorCode.ErrorReadingData: return "An error occoured while reading from the input file.";
        case ErrorCode.ErrorInvalidDefinitionN: return "Check your definition for the number of points.";
        case ErrorCode.ErrorInvalidDefinitionC: return "Check your definition for the number of constraints.";
        case ErrorCode.ErrorInvalidDefinitionU: return "Check your definition for upper bound constraints.";
        case ErrorCode.ErrorInvalidDefinitionL: return "Check your definition for lower bound constraints.";
        case ErrorCode.ErrorInvalidDefinitionX: return "Check your definition for maximization constraints.";
        case ErrorCode.ErrorInvalidDefinitionI: return "Check your definition for minimization constraints.";
        case ErrorCode.ErrorInvalidDefinitionP: return "At least one point definition is missing something.";
        case ErrorCode.ErrorInvalidDefinitionS: return "Check your definition for the number of scores.";
        case ErrorCode.ErrorInvalidDefinitionB: return "Check your definition for the starting point.";
        case ErrorCode.ErrorInvalidDefinitionE: return "Check your definition for the terminus point.";
        case ErrorCode.ErrorInvalidDefinitionD: return "Check your definition for the distance matrix.";
        case ErrorCode.ErrorInvalidDefinitionM: return "Check your definition for the distance matrix flag.";
        case ErrorCode.ErrorMissingStartingPoint: return "The definition of the starting point is missing.";
        case ErrorCode.ErrorMissingTerminusPoint: return "The definition of the terminus point is missing.";
        case ErrorCode.ErrorTooFewConstraints: return "There are too few constraints on the list.";
        case ErrorCode.ErrorTooFewPoints: return "There are too few points on the list.";
        case ErrorCode.ErrorTooFewScores: return "There are too few scores on the list of at least one point.";
        case ErrorCode.ErrorTooManyConstraints: return "There are too many constraints on the list.";
        case ErrorCode.ErrorTooManyPoints: return "There are too many points on the list.";
        case ErrorCode.ErrorTooManyScores: return "There are too many scores on the list of at least one point.";
        default: return "Unknown error.";
      }
    }

    #endregion

    #region VNS

    private void RunVns(ProblemData problem, long startTime, long maxRuntime) {
      long iterationCounter = 0;
      int neighborhoodMax;
      long iterationOfLastAcception = 0;
      long iterationOfLastFinding = 0;
      var preferences = new List<double>();
      bool increasingPreference = false;
      TourData newSolution;

      // Define the maximum neighborhood distance
      //neighborhoodMax = 9;

      // Generate an initial solution
      GenerateGreedyTour(problem);

      // Bring this tour to a local optimum
      OptimizePath(problem, problem.VNS.InitialTour);

      // Save this tour as at the moment it seems efficient
      SaveDominantToursVns(problem, problem.VNS.InitialTour, startTime);

      // Create the starting preference for the shaking
      preferences.Add(1);
      preferences.Add(0);

      // Perform the VNS
      while ((DateTime.Now.Ticks - startTime) < maxRuntime) {
        int neighborhoodCounter = 1;
        //while ((neighborhoodCounter < neighborhoodMax) && ((DateTime.Now.Ticks - start_time) < max_runtime))
        while ((neighborhoodCounter < problem.VNS.InitialTour.PathSize - 2) && ((DateTime.Now.Ticks - startTime) < maxRuntime)) {
          iterationCounter++;

          // Generate a random tour in the kth neighborhood
          if ((iterationCounter - iterationOfLastFinding) > 1000) {
            // Create the next preference vector for shaking
            if (problem.ScoreCount == 2) {
              double temp = problem.MTR.NextDouble() / 1000.0;
              if (increasingPreference) {
                preferences[0] += temp;
                if (preferences[0] >= 1.0) {
                  preferences[0] = 1.0;
                  increasingPreference = false;
                }
                preferences[1] = 1.0 - preferences[0];
              } else {
                preferences[1] += temp;
                if (preferences[1] >= 1.0) {
                  preferences[1] = 1.0;
                  increasingPreference = true;
                }
                preferences[0] = 1.0 - preferences[1];
              }
            }
          }

          // Shake forcing one of the scores
          ShakeTour(problem, neighborhoodCounter, preferences);

          // Bring the tour back to be feasible 
          CleanupTour(problem);

          // Bring this tour to a local optimum
          OptimizePath(problem, problem.VNS.ActualTour);


          for (int i = 0; i < problem.ScoreCount; i++)
            scores.Rows["Score " + i].Values.Add(problem.VNS.ActualTour.Scores[i]);

          // Save the tour if it is pareto dominant
          if (SaveDominantToursVns(problem, problem.VNS.ActualTour, startTime)) {
            // Move to the better solution and return to the nearest neighborhood
            problem.VNS.InitialTour = new TourData(problem.VNS.ActualTour);
            iterationOfLastFinding = iterationCounter;
            neighborhoodCounter = 0;
          } else if ((iterationCounter - iterationOfLastAcception) >= problem.VNS.AcceptFrequency) {
            // Sometimes accept a solution that is not more than a given ratio below
            // any of the dominant solutions found so far
            bool allow = true;
            int tourMax = problem.VNS.DominantTours.Count;
            for (int tourCounter = 0; tourCounter < tourMax; tourCounter++) {
              for (int scoreCounter = 0; scoreCounter < problem.ScoreCount; scoreCounter++) {
                if ((problem.VNS.DominantTours[tourCounter].Scores[scoreCounter] * problem.VNS.TresholdValue) > problem.VNS.ActualTour.Scores[scoreCounter]) {
                  allow = false;
                  break;
                }
              }
              if (!allow) break;
            }
            if (allow) {
              // Move to the worse solution and return to the nearest neighborhood
              problem.VNS.InitialTour = new TourData(problem.VNS.ActualTour);
              iterationOfLastAcception = iterationCounter;
              neighborhoodCounter = 0;
            }
          }

          // Increase the neighborhood to search
          neighborhoodCounter++;
        }
      }

      Console.WriteLine("Time Of Adding of Best Tour: {0}", problem.VNS.DominantTours[0].TimeOfAdding.TotalSeconds);
      //printf("DEBUG: iteration_counter = %d\n", iteration_counter);
    }

    private static void GenerateGreedyTour(ProblemData problem) {
      int constraintCounter;
      double maxDistance = 0.0;
      var feasiblePoints = new List<FeasiblePointVns>();

      // Find the maximum tour length
      for (constraintCounter = 0; constraintCounter < problem.ConstraintCount; constraintCounter++) {
        if ((problem.Constraints[constraintCounter].Type == Constraint.Upperbound) && (problem.Constraints[constraintCounter].Target == ConstraintTimeDistance)) {
          maxDistance = problem.Constraints[constraintCounter].TargetValue;
          break;
        }
      }

      int numberFeasiblePoints = 0;
      // Find all points within the maximum distance allowed (ellipse)
      for (int pointCounter = 0; pointCounter < problem.PointCount; pointCounter++) {
        double distance = problem.Edges[problem.StartingPoint][pointCounter].Distance + problem.Edges[pointCounter][problem.TerminusPoint].Distance + problem.FixedPenalty;
        // DEBUG
        //printf("Distance = %.9f, max_dist = %.9f\n", distance, max_distance);
        //if ((distance <= max_distance) && (point_counter != Problem.Starting_Point) && (point_counter != Problem.Terminus_Point))
        if (((maxDistance - distance) >= PrecisionUsed) && (pointCounter != problem.StartingPoint) && (pointCounter != problem.TerminusPoint)) {
          var tempPoint = new FeasiblePointVns {
            Distance = distance,
            Index = pointCounter,
            Scores = new List<int>(problem.Points[pointCounter].Score)
          };

          // Insert the point at the corresponding position
          int sortMax = feasiblePoints.Count;
          int sortCounter;
          for (sortCounter = 0; sortCounter < sortMax; sortCounter++) {
            int totalScoreNew = 0;
            int totalScoreOld = 0;
            for (int counterScores = 0; counterScores < problem.ScoreCount; counterScores++) {
              totalScoreNew += feasiblePoints[sortCounter].Scores[counterScores];
              totalScoreOld += tempPoint.Scores[counterScores];
            }
            if (totalScoreNew < totalScoreOld) {
              feasiblePoints.Insert(sortCounter, tempPoint);
              numberFeasiblePoints++;
              distance = -1.0;
              break;
            }
          }

          // If the point was not yet inserted, append it
          if (distance > 0) {
            feasiblePoints.Add(tempPoint);
            numberFeasiblePoints++;
          }
        }
      }

      // Add the starting and terminus point
      problem.VNS.InitialTour.Path = new List<int> {
        problem.StartingPoint, 
        problem.TerminusPoint
      };
      problem.VNS.InitialTour.PathSize = 2;
      problem.VNS.InitialTour.Length = CalcLength(problem, problem.VNS.InitialTour.Path);
      problem.VNS.InitialTour.Scores = new List<int>(Enumerable.Repeat(0, problem.ScoreCount));

      // Add points in a greedy way
      bool insertionPerformed = true;
      while (insertionPerformed) {
        insertionPerformed = false;

        for (int pointCounter = 0; pointCounter < feasiblePoints.Count; pointCounter++) {
          for (int positionCounter = 1; positionCounter < problem.VNS.InitialTour.Path.Count; positionCounter++) {
            // Create the candidate tour
            TourData tempTour = new TourData(problem.VNS.InitialTour);
            tempTour.Path.Insert(positionCounter, feasiblePoints[pointCounter].Index);
            tempTour.Length = CalcLength(problem, tempTour.Path);

            // If the insertion would be feasible, perform it
            if ((maxDistance - tempTour.Length) >= PrecisionUsed) {//(Temp_Tour.Length <= max_distance)
              problem.VNS.InitialTour.Length = tempTour.Length;
              problem.VNS.InitialTour.Path = new List<int>(tempTour.Path);
              problem.VNS.InitialTour.PathSize++;
              for (int scoreCounter = 0; scoreCounter < problem.ScoreCount; scoreCounter++) {
                problem.VNS.InitialTour.Scores[scoreCounter] += feasiblePoints[pointCounter].Scores[scoreCounter];
              }
              feasiblePoints.RemoveAt(pointCounter);
              insertionPerformed = true;
              break;
            }
          }
          if (insertionPerformed) break;
        }
      }
    }

    private static void ShakeTour(ProblemData problem, int maxNeighborhood, List<double> preferences) {
      var visitablePoints = new List<int>();

      // Limit the neighborhood to the tour length
      int limit = problem.VNS.InitialTour.PathSize - 3;
      int neighborhood = problem.MTR.Next((limit > maxNeighborhood) ? maxNeighborhood : limit) + 1;

      // Determine the maximum travelling distance
      double maxDistance = problem.Constraints[0].TargetValue;

      // Elect the starting index of the part to be replaced
      int tourSize = problem.VNS.InitialTour.PathSize;
      int randomPosition = problem.MTR.Next(tourSize - neighborhood - 2) + 1;

      // Initialize the new tour
      problem.VNS.ActualTour = new TourData {
        Length = 0.0,
        Path = new List<int> { problem.StartingPoint },
        PathSize = 1,
        Scores = new List<int>(Enumerable.Repeat(0, problem.ScoreCount))
      };

      // Find all points that are not yet included in the tour and are
      // within the maximum distance allowed (ellipse) and sort them with
      // regard to their utility 
      visitablePoints.Clear();
      int numberVisitablePoints = 0;
      int pathSize = problem.VNS.ActualTour.PathSize;
      for (int pointCounter = 0; pointCounter < problem.PointCount; pointCounter++) {
        // Calculate the distance when going from the starting
        // point to this point and then to the end point
        double distance = problem.Edges[problem.StartingPoint][pointCounter].Distance + problem.Edges[pointCounter][problem.TerminusPoint].Distance + problem.FixedPenalty;

        // If this distance is feasible and the point is
        // neither starting nor ending point, check the point
        if (((maxDistance - distance) >= PrecisionUsed) && (pointCounter != problem.StartingPoint) && (pointCounter != problem.TerminusPoint)) {
          bool alreadyVisited = false;
          int tourPointCounter;
          for (tourPointCounter = 0; tourPointCounter < tourSize; tourPointCounter++) {
            // Check if the point is included in the initial tour
            if (pointCounter == problem.VNS.InitialTour.Path[tourPointCounter]) {
              alreadyVisited = true;
              break;
            }
          }

          // The point was not yet visited, so add it to the candidate list
          if (!alreadyVisited) {
            double utilityPoint = 0;
            bool inserted = false;

            // Calculate the current point's utility
            for (int counterScore = 0; counterScore < problem.ScoreCount; counterScore++) {
              utilityPoint += problem.Points[pointCounter].Score[counterScore] * preferences[counterScore];
            }

            // Add the point at the right position according to its utility (increasing)
            for (int counterPosition = 0; counterPosition < numberVisitablePoints; counterPosition++) {
              // Calculate the utility of the point at this position
              double utilityPosition = 0;
              for (int counterScore = 0; counterScore < problem.ScoreCount; counterScore++) {
                utilityPosition += problem.Points[visitablePoints[counterPosition]].Score[counterScore] * preferences[counterScore];
              }

              // If the utility of the point is lower, insert it here
              if (utilityPosition > utilityPoint) {
                visitablePoints.Insert(counterPosition, pointCounter);
                inserted = true;
                break;
              }
            }

            // If no bigger utility exists up to now, add the point at the end
            if (!inserted) {
              visitablePoints.Add(pointCounter);
            }
            numberVisitablePoints++;
          }
        }
      }

      // Perform the insertions according to the utility sorting
      for (int position = 1; position < tourSize; position++) {
        int scoreCounter;
        if ((position < randomPosition) || (position > (randomPosition + neighborhood - 1))) {
          // Add the waypoints to the new solution
          problem.VNS.ActualTour.Path.Add(problem.VNS.InitialTour.Path[position]);
          problem.VNS.ActualTour.PathSize++;
          problem.VNS.ActualTour.Length = CalcLength(problem, problem.VNS.ActualTour.Path);
          for (scoreCounter = 0; scoreCounter < problem.ScoreCount; scoreCounter++) {
            problem.VNS.ActualTour.Scores[scoreCounter] += problem.Points[problem.VNS.InitialTour.Path[position]].Score[scoreCounter];
          }

          // Delete this point from the candidate list
          for (int pointCounter = 0; pointCounter < numberVisitablePoints; pointCounter++) {
            if (problem.VNS.InitialTour.Path[position] == visitablePoints[pointCounter]) {
              visitablePoints.RemoveAt(pointCounter);
              numberVisitablePoints--;
              break;
            }
          }
        } else {
          if (numberVisitablePoints > 0) {
            // Add the point with the highest utility from the candidate list
            int randomFactor = problem.MTR.Next(3);
            int insertionIndex = numberVisitablePoints - 1;
            if (numberVisitablePoints > 4) insertionIndex -= randomFactor;
            problem.VNS.ActualTour.Path.Add(visitablePoints[insertionIndex]); //[number_visitable_points - 1]);
            problem.VNS.ActualTour.PathSize++;
            problem.VNS.ActualTour.Length = CalcLength(problem, problem.VNS.ActualTour.Path);
            for (scoreCounter = 0; scoreCounter < problem.ScoreCount; scoreCounter++) {
              //Problem.VNS.Actual_Tour.Scores[score_counter] += Problem.Points[visitable_points[number_visitable_points - 1]].Score[score_counter];
              problem.VNS.ActualTour.Scores[scoreCounter] += problem.Points[visitablePoints[insertionIndex]].Score[scoreCounter];
            }
            //visitable_points.pop_back();
            visitablePoints.RemoveAt(insertionIndex);
            numberVisitablePoints--;
          } else {
            // We don't have any points left that could be inserted so we can only re-insert
            // the removed and not already re-inserted points in a random order
            for (int reinsertCounter = randomPosition; reinsertCounter < (randomPosition + neighborhood); reinsertCounter++) {
              bool alreadyReinserted = false;
              for (int checkCounter = 0; checkCounter < problem.VNS.ActualTour.PathSize; checkCounter++) {
                if (problem.VNS.InitialTour.Path[reinsertCounter] == problem.VNS.ActualTour.Path[checkCounter]) {
                  alreadyReinserted = true;
                  break;
                }
              }
              if (!alreadyReinserted) {
                visitablePoints.Add(problem.VNS.InitialTour.Path[reinsertCounter]);
                numberVisitablePoints++;
              }
            }

            int randomIndex = problem.MTR.Next(numberVisitablePoints - 1);
            problem.VNS.ActualTour.Path.Add(visitablePoints[randomIndex]);
            problem.VNS.ActualTour.PathSize++;
            problem.VNS.ActualTour.Length = CalcLength(problem, problem.VNS.ActualTour.Path);
            for (scoreCounter = 0; scoreCounter < problem.ScoreCount; scoreCounter++) {
              problem.VNS.ActualTour.Scores[scoreCounter] += problem.Points[visitablePoints[randomIndex]].Score[scoreCounter];
            }
            visitablePoints.Clear();
            numberVisitablePoints = 0;
          }
        }
      }
    }

    private static void CleanupTour(ProblemData problem) {
      int tourPointCounter, sortCounter;
      var distanceSavingsList = new List<double>();
      var sortedPoints = new List<int>();

      // Find the maximum distance restriction
      double maxLength = problem.Constraints[0].TargetValue;

      distanceSavingsList.Clear();
      sortedPoints.Clear();

      // Sort the points on the tour according to their average score
      int tourPointMax = problem.VNS.ActualTour.PathSize - 1;
      for (tourPointCounter = 1; tourPointCounter < tourPointMax; tourPointCounter++) {
        var tempPath = new List<int>();

        double distanceSaving = problem.Edges[problem.VNS.ActualTour.Path[tourPointCounter - 1]][problem.VNS.ActualTour.Path[tourPointCounter]].Distance;
        distanceSaving += problem.Edges[problem.VNS.ActualTour.Path[tourPointCounter]][problem.VNS.ActualTour.Path[tourPointCounter + 1]].Distance;
        distanceSaving -= problem.Edges[problem.VNS.ActualTour.Path[tourPointCounter - 1]][problem.VNS.ActualTour.Path[tourPointCounter + 1]].Distance;

        bool inserted = false;
        int sortMax = sortedPoints.Count;
        for (sortCounter = 0; sortCounter < sortMax; sortCounter++) {
          if (distanceSaving > distanceSavingsList[sortCounter]) {
            distanceSavingsList.Insert(sortCounter, distanceSaving);
            sortedPoints.Insert(sortCounter, tourPointCounter);
            inserted = true;
            break;
          }
        }
        if (!inserted) {
          distanceSavingsList.Add(distanceSaving);
          sortedPoints.Add(tourPointCounter);
        }
      }

      // As long as the created path is infeasible, remove elements
      while ((problem.VNS.ActualTour.Length - maxLength) >= PrecisionUsed) { //(Problem.VNS.Actual_Tour.Length > max_length)
        // Remove the point that frees the largest distance
        int scoreCounter;
        for (scoreCounter = 0; scoreCounter < problem.ScoreCount; scoreCounter++) {
          problem.VNS.ActualTour.Scores[scoreCounter] -= problem.Points[problem.VNS.ActualTour.Path[sortedPoints[0]]].Score[scoreCounter];
        }
        problem.VNS.ActualTour.Path.RemoveAt(sortedPoints[0]);
        problem.VNS.ActualTour.PathSize--;
        problem.VNS.ActualTour.Length = CalcLength(problem, problem.VNS.ActualTour.Path);
        int sortMax = sortedPoints.Count;
        for (sortCounter = 1; sortCounter < sortMax; sortCounter++) {
          if (sortedPoints[sortCounter] > sortedPoints[0]) {
            sortedPoints[sortCounter]--;
          }
        }
        sortedPoints.RemoveAt(0);
        distanceSavingsList.RemoveAt(0);
      }
    }

    private static bool SaveDominantToursVns(ProblemData problem, TourData tour, long startTime) {
      int tourCounter;

      bool processed = false;
      int tourCounterMax = problem.VNS.DominantTours.Count - 1;

      for (tourCounter = tourCounterMax; tourCounter >= 0; --tourCounter) {
        bool scoreGt = false;
        bool scoreLt = false;
        int scoreCounter = 0;
        for (scoreCounter = (problem.ScoreCount - 1); scoreCounter >= 0; --scoreCounter) {
          if (tour.Scores[scoreCounter] > problem.VNS.DominantTours[tourCounter].Scores[scoreCounter]) {
            scoreGt = true;
          } else if (tour.Scores[scoreCounter] < problem.VNS.DominantTours[tourCounter].Scores[scoreCounter]) {
            scoreLt = true;
          }
        }

        // If this tour equals an already existing one, skip it
        if (tour.Path == problem.VNS.DominantTours[tourCounter].Path) {
          processed = true;
          break;
        }

        // If the tour is dominated by an existing, skip it
        if (scoreLt && !scoreGt) {
          processed = true;
          break;
        }

        // If the tour dominates an existing one, eliminate the existing one
        if (!scoreLt && scoreGt) {
          problem.VNS.DominantTours.RemoveRange(tourCounter, 1);
          tourCounterMax--;
        }
      }
      if (!processed) {
        // The tour either dominates an existing one or is just not dominated 
        // by an existing tour, so we assume it to be dominant
        tour.TimeOfAdding = CalculateDifference(startTime, DateTime.Now.Ticks);
        problem.VNS.DominantTours.Add(tour);
        return true;
      } else {
        return false;
      }
    }

    #endregion

    #region Optimization

    private static void ShortenPath(ProblemData problem, TourData tour) {
      bool optimizationDone = true;
      int pathSize = tour.PathSize;
      int maxBlockLength = (pathSize > 31) ? 30 : (pathSize - 2);

      // DEBUG
      //printf("Shorten:");

      // Perform a 2-opt
      while (optimizationDone) {
        optimizationDone = false;

        int blockLength;
        for (blockLength = 2; blockLength < maxBlockLength; blockLength++) { //(path_size - 1); ++block_length)
          // If an optimization has been done, start from the beginning
          if (optimizationDone) break;

          // DEBUG
          //printf("\tBlock length = %d, Max block length = %d.\n", block_length, max_block_length);

          int position;
          for (position = 1; position < (pathSize - blockLength); position++) {
            // If an optimization has been done, start from the beginning
            if (optimizationDone) break;

            // DEBUG
            //printf("\t\tPosition = %d, Path size = %d.\n", position, path_size);

            double newLength = tour.Length;
            for (int counterLength = (position - 1); counterLength < (position + blockLength); counterLength++) {
              newLength -= problem.Edges[tour.Path[counterLength]][tour.Path[counterLength + 1]].Distance;
            }
            for (int counterLength = (position + blockLength - 1); counterLength > (position); counterLength--) {
              newLength += problem.Edges[tour.Path[counterLength]][tour.Path[counterLength - 1]].Distance;
            }
            newLength += problem.Edges[tour.Path[position - 1]][tour.Path[position + blockLength - 1]].Distance;
            newLength += problem.Edges[tour.Path[position]][tour.Path[position + blockLength]].Distance;


            if (newLength <= (tour.Length - 0.00001)) { // Avoid cycling caused by precision
              var parthPath = tour.Path.GetRange(position, blockLength);

              tour.Path.RemoveRange(position, blockLength);
              tour.Path.InsertRange(position, parthPath.AsEnumerable().Reverse());

              tour.Length = CalcLength(problem, tour.Path);

              // Re-run the optimization
              optimizationDone = true;
            }
          }
        }
      }
      // DEBUG
      //printf("Shorten done.\n");
    }

    private static void OptimizePath(ProblemData problem, TourData tour, int optimizationPreference = -1) {
      bool optimizationDone = true;
      int acceptNonParetoExchange = 0;
      bool isFeasible = true;
      var visitablePoints = new List<int>();

      // Find the maximum distance restriction
      double maxLength = problem.Constraints[0].TargetValue;

      // Check if the tour can be improved by adding or replacing points
      while (optimizationDone) {
        optimizationDone = false;

        // Try to shorten the path
        ShortenPath(problem, tour);

        // Determine the number of points visited so far
        int tourSize = tour.PathSize; //, constraint_counter;

        // Determine all points that have not yet been visited by this tour
        visitablePoints.Clear();
        int numberVisitablePoints = 0; //, constraint_counter;
        int tourPointCounter; //, constraint_counter;
        int pointCounter; //, constraint_counter;
        for (pointCounter = 0; pointCounter < problem.PointCount; pointCounter++) {
          bool alreadyVisited = false;
          for (tourPointCounter = 0; tourPointCounter < tourSize; tourPointCounter++) {
            if (pointCounter == tour.Path[tourPointCounter]) {
              alreadyVisited = true;
              break;
            }
          }
          if (!alreadyVisited) {
            visitablePoints.Add(pointCounter);
            numberVisitablePoints++;
          }
        }

        // Determine if any of the visitable points can be included at any
        // position within the tour
        int scoreCounter; //, constraint_counter;
        for (tourPointCounter = 1; tourPointCounter < tourSize; tourPointCounter++) {
          // If an optimization has been done, start from the beginning
          if (optimizationDone) break;

          for (pointCounter = 0; pointCounter < numberVisitablePoints; pointCounter++) {
            // If an optimization has been done, start from the beginning
            if (optimizationDone) break;

            double newLength = tour.Length;
            newLength -= problem.Edges[tour.Path[tourPointCounter - 1]][tour.Path[tourPointCounter]].Distance;
            newLength += problem.Edges[tour.Path[tourPointCounter - 1]][visitablePoints[pointCounter]].Distance;
            newLength += problem.Edges[visitablePoints[pointCounter]][tour.Path[tourPointCounter]].Distance;
            newLength += problem.FixedPenalty;

            // Determine if including the point does not violate any constraint
            if (newLength <= maxLength) {
              // We do not need to check if the point increases the score - it can only increase!

              // Update the scores achieved by visiting this point
              for (scoreCounter = (problem.ScoreCount - 1); scoreCounter >= 0; --scoreCounter) {
                tour.Scores[scoreCounter] += problem.Points[visitablePoints[pointCounter]].Score[scoreCounter];
              }

              // Insert the new point at this position
              tour.Path.Insert(tourPointCounter, visitablePoints[pointCounter]);
              tour.PathSize++;

              // Update the overall tour length
              tour.Length = CalcLength(problem, tour.Path);

              // Re-run this optimization
              optimizationDone = true;
            }
          }
        }

        // Determine if any of the visitable points can take the place of an
        // already visited point in the tour to improve the scores
        for (tourPointCounter = 1; tourPointCounter < (tourSize - 1); tourPointCounter++) {
          // If an optimization has been done, start from the beginning
          if (optimizationDone) break;

          for (pointCounter = 0; pointCounter < numberVisitablePoints; pointCounter++) {
            // If an optimization has been done, start from the beginning
            if (optimizationDone) break;

            double newLength = tour.Length;
            newLength -= problem.Edges[tour.Path[tourPointCounter - 1]][tour.Path[tourPointCounter]].Distance;
            newLength -= problem.Edges[tour.Path[tourPointCounter]][tour.Path[tourPointCounter + 1]].Distance;
            newLength += problem.Edges[tour.Path[tourPointCounter - 1]][visitablePoints[pointCounter]].Distance;
            newLength += problem.Edges[visitablePoints[pointCounter]][tour.Path[tourPointCounter + 1]].Distance;

            if (newLength <= maxLength) {
              // Check if the point dominated the one it would replace
              bool scoreGt = false;
              bool scoreLt = false;
              if (optimizationPreference == -1) {
                for (scoreCounter = (problem.ScoreCount - 1); scoreCounter >= 0; --scoreCounter) {
                  if (problem.Points[visitablePoints[pointCounter]].Score[scoreCounter] > problem.Points[tour.Path[tourPointCounter]].Score[scoreCounter])
                    scoreGt = true;
                  else if (problem.Points[visitablePoints[pointCounter]].Score[scoreCounter] < problem.Points[tour.Path[tourPointCounter]].Score[scoreCounter])
                    scoreLt = true;
                }
              } else {
                if (problem.Points[visitablePoints[pointCounter]].Score[optimizationPreference] > problem.Points[tour.Path[tourPointCounter]].Score[optimizationPreference])
                  scoreGt = true;
                else if (problem.Points[visitablePoints[pointCounter]].Score[optimizationPreference] < problem.Points[tour.Path[tourPointCounter]].Score[optimizationPreference])
                  scoreLt = true;
              }

              if (!scoreLt && scoreGt) {
                // Update the overall tour length
                // Tour.Length = new_length;
                // Update the scores achieved by visiting this point
                for (scoreCounter = (problem.ScoreCount - 1); scoreCounter >= 0; --scoreCounter) {
                  tour.Scores[scoreCounter] -= problem.Points[tour.Path[tourPointCounter]].Score[scoreCounter];
                  tour.Scores[scoreCounter] += problem.Points[visitablePoints[pointCounter]].Score[scoreCounter];
                }

                // Replace the old point by the new one
                tour.Path.RemoveAt(tourPointCounter);
                tour.Path.Insert(tourPointCounter, visitablePoints[pointCounter]);

                // Update the overall tour length
                tour.Length = CalcLength(problem, tour.Path);

                // Re-run this optimization
                optimizationDone = true;
              }

              // If no dominant exchange could be performed, accept a non dominant exchange
              if (scoreGt && (acceptNonParetoExchange > 0) && (acceptNonParetoExchange < 6) && (optimizationPreference == -1)) {
                // Update the overall tour length
                //Tour.Length = new_length;

                // Update the scores achieved by visiting this point
                for (scoreCounter = (problem.ScoreCount - 1); scoreCounter >= 0; --scoreCounter) {
                  tour.Scores[scoreCounter] -= problem.Points[tour.Path[tourPointCounter]].Score[scoreCounter];
                  tour.Scores[scoreCounter] += problem.Points[visitablePoints[pointCounter]].Score[scoreCounter];
                }

                // Replace the old point by the new one
                tour.Path.RemoveAt(tourPointCounter);
                tour.Path.Insert(tourPointCounter, visitablePoints[pointCounter]);

                // Update the overall tour length
                tour.Length = CalcLength(problem, tour.Path);

                // Re-run this optimization
                optimizationDone = true;
              }
            }
          }
        }

        if (!optimizationDone && (acceptNonParetoExchange == 0)) {
          // If no optimization has been done, accept non dominant exchanges
          // Non dominant exchanges will only be accepted 5 times to prevent infinite loops
          optimizationDone = true;
          acceptNonParetoExchange++;
        } else if (optimizationDone && (acceptNonParetoExchange > 0)) {
          // If an optimization has been done, do not accept non dominant exchange
          acceptNonParetoExchange++;
        }
      }
    }

    private static double CalcLength(ProblemData problem, List<int> path) {
      double length = 0.0;

      for (int pointCounter = 1; pointCounter < path.Count; pointCounter++) {
        length += problem.Edges[path[pointCounter - 1]][path[pointCounter]].Distance;
      }
      // Add the fixed penalty for every vertex except for the starting and ending vertex
      if (path.Count > 1) {
        if (path[path.Count - 1] == problem.TerminusPoint) length += ((path.Count() - 2.0) * problem.FixedPenalty);
        else length += ((path.Count - 1.0) * problem.FixedPenalty);
      }

      return length;
    }

    #endregion

    #region Path Relinking

    private static int PathRelinking(ProblemData problem, List<TourData> tours, long startTime) {
      int numberOfToursFound = 0;

      // Sort the tours with regard to their scores
      SortTours(tours);

      for (int counterGuidingTour = 1; counterGuidingTour < tours.Count; counterGuidingTour++) {
        TourData currentTour = tours[counterGuidingTour - 1];
        TourData guidingTour = tours[counterGuidingTour];
        numberOfToursFound += PathRelinkingStep(0, problem, tours, currentTour, guidingTour, startTime);
      }

      return numberOfToursFound;
    }

    private static void SortTours(List<TourData> tours) {
      int numberTours = tours.Count;

      // Sort the tours with regard to their scores
      for (int counterTours = (numberTours - 1); counterTours >= 0; counterTours--) {
        for (int counterSearch = 1; counterSearch <= counterTours; counterSearch++) {
          if (tours[counterSearch].Scores[0] < tours[counterSearch - 1].Scores[0]) {
            tours.Insert(counterSearch - 1, tours[counterSearch]);
            tours.RemoveAt(counterSearch + 1);
          }
        }
      }
    }

    private static int PathRelinkingStep(int level, ProblemData problem, List<TourData> tours, TourData currentTour, TourData guidingTour, long startTime) {
      int numberOfToursFound = 0;
      List<int> nodesToRemove = new List<int>();
      List<int> nodesToInsert = new List<int>();
      int numberOfNodesToRemove = FindNodesToRemove(currentTour, guidingTour, nodesToRemove);
      int numberOfNodesToInsert = FindNodesToInsert(currentTour, guidingTour, nodesToInsert);

      // If no other transformations are possible without reaching the guiding solution, return
      if ((numberOfNodesToRemove <= 1) && (numberOfNodesToInsert <= 1)) {
        return 0;
      }

      if (numberOfNodesToRemove != 0 && numberOfNodesToInsert != 0) {
        // If too many combinations can be performed on this layer, draw one randomly
        if ((numberOfNodesToRemove * numberOfNodesToInsert) < 16) {
          // Nodes for removing and inserting are available
          for (int counterRemove = 0; counterRemove < numberOfNodesToRemove; counterRemove++) {
            for (int counterInsert = 0; counterInsert < numberOfNodesToInsert; counterInsert++) {
              // Make a working copy of the current tour
              TourData candidateTour = new TourData(currentTour);

              // Perform the transformation step
              for (int counterScore = 0; counterScore < problem.ScoreCount; counterScore++) {
                candidateTour.Scores[counterScore] -= problem.Points[candidateTour.Path[nodesToRemove[counterRemove]]].Score[counterScore];
                candidateTour.Scores[counterScore] += problem.Points[guidingTour.Path[nodesToInsert[counterInsert]]].Score[counterScore];
              }
              candidateTour.Path.RemoveAt(nodesToRemove[counterRemove]);
              candidateTour.Path.Insert(nodesToRemove[counterRemove], guidingTour.Path[nodesToInsert[counterInsert]]);
              candidateTour.Length = CalcLength(problem, candidateTour.Path);

              // Locally optimize the modified tour
              TourData tempTour = new TourData(candidateTour);
              OptimizePath(problem, tempTour);

              // If the tour is efficient, save it and increase the number of tours found
              if (SaveCandidateTour(problem, tours, tempTour, startTime)) numberOfToursFound++;

              // Recursively perform all remaining steps after this one
              if ((numberOfNodesToRemove > 1) || (numberOfNodesToInsert > 1)) {
                numberOfToursFound += PathRelinkingStep((level + 1), problem, tours, candidateTour, guidingTour, startTime);
              }
            }
          }
        } else {
          // Draw one combination randomly
          int nodeToRemove = problem.MTR.Next(numberOfNodesToRemove - 1);
          int nodeToInsert = problem.MTR.Next(numberOfNodesToInsert - 1);

          // Make a working copy of the current tour
          TourData candidateTour = currentTour;

          // Perform the transformation step
          for (int counterScore = 0; counterScore < problem.ScoreCount; counterScore++) {
            candidateTour.Scores[counterScore] -= problem.Points[candidateTour.Path[nodesToRemove[nodeToRemove]]].Score[counterScore];
            candidateTour.Scores[counterScore] += problem.Points[guidingTour.Path[nodesToInsert[nodeToInsert]]].Score[counterScore];
          }
          candidateTour.Path.RemoveAt(nodesToRemove[nodeToRemove]);
          candidateTour.Path.Insert(nodesToRemove[nodeToRemove], guidingTour.Path[nodesToInsert[nodeToInsert]]);
          candidateTour.Length = CalcLength(problem, candidateTour.Path);

          // Locally optimize the modified tour
          TourData tempTour = new TourData(candidateTour);
          OptimizePath(problem, tempTour);

          // If the tour is efficient, save it and increase the number of tours found
          if (SaveCandidateTour(problem, tours, tempTour, startTime)) numberOfToursFound++;

          // Recursively perform all remaining steps after this one
          if ((numberOfNodesToRemove > 1) || (numberOfNodesToInsert > 1)) {
            numberOfToursFound += PathRelinkingStep((level + 1), problem, tours, candidateTour, guidingTour, startTime);
          }
        }
      } else if (numberOfNodesToRemove != 0) {
        // Only removing nodes is possible
        for (int counter_remove = 0; counter_remove < numberOfNodesToRemove; counter_remove++) {
          // Make a working copy of the current tour
          TourData candidateTour = new TourData(currentTour);

          // Perform the transformation step
          for (int counterScore = 0; counterScore < problem.ScoreCount; counterScore++) {
            candidateTour.Scores[counterScore] -= problem.Points[candidateTour.Path[nodesToRemove[counter_remove]]].Score[counterScore];
          }
          candidateTour.Path.RemoveAt(nodesToRemove[counter_remove]);
          candidateTour.PathSize--;
          candidateTour.Length = CalcLength(problem, currentTour.Path);

          // Locally optimize the modified tour
          TourData tempTour = new TourData(candidateTour);
          OptimizePath(problem, tempTour);

          // If the tour is efficient, save it and increase the number of tours found
          if (SaveCandidateTour(problem, tours, tempTour, startTime)) numberOfToursFound++;
        }
      } else if (numberOfNodesToInsert != 0) {
        // Only inserting nodes is possible
        for (int counterInsert = 0; counterInsert < numberOfNodesToInsert; counterInsert++) {
          // Make a working copy of the current tour
          TourData candidateTour = new TourData(currentTour);

          // Perform the transformation step
          for (int counterScore = 0; counterScore < problem.ScoreCount; counterScore++) {
            candidateTour.Scores[counterScore] += problem.Points[guidingTour.Path[nodesToInsert[counterInsert]]].Score[counterScore];
          }
          candidateTour.Path.Insert((candidateTour.PathSize / 2), guidingTour.Path[nodesToInsert[counterInsert]]);
          candidateTour.PathSize++;
          candidateTour.Length = CalcLength(problem, candidateTour.Path);

          // Locally optimize the modified tour
          TourData tempTour = new TourData(candidateTour);
          OptimizePath(problem, tempTour);

          // If the tour is efficient, save it and increase the number of tours found
          if (SaveCandidateTour(problem, tours, tempTour, startTime)) numberOfToursFound++;
        }
      }

      return numberOfToursFound;
    }

    private static int FindNodesToRemove(TourData currentTour, TourData guidingTour, List<int> nodesToRemove) {
      // Find all node to remove from the current tour
      for (int counterNodes = 1; counterNodes < currentTour.PathSize; counterNodes++) {
        bool nodeFound = false;
        for (int counterSearch = 1; counterSearch < guidingTour.PathSize; counterSearch++) {
          if (guidingTour.Path[counterSearch] == currentTour.Path[counterNodes]) {
            nodeFound = true;
            break;
          }
        }
        if (!nodeFound) {
          // Remember the position of the node that has to be removed
          nodesToRemove.Add(counterNodes);
        }
      }
      return nodesToRemove.Count;
    }

    private static int FindNodesToInsert(TourData currentTour, TourData guidingTour, List<int> nodesToInsert) {
      // Find all node to insert into the current tour
      for (int counterNodes = 1; counterNodes < guidingTour.PathSize; counterNodes++) {
        bool nodeFound = false;
        for (int counterSearch = 1; counterSearch < currentTour.PathSize; counterSearch++) {
          if (currentTour.Path[counterSearch] == guidingTour.Path[counterNodes]) {
            nodeFound = true;
            break;
          }
        }
        if (!nodeFound) {
          // Remember the position of the node that has to be removed
          nodesToInsert.Add(counterNodes);
        }
      }
      return nodesToInsert.Count;
    }

    private static bool SaveCandidateTour(ProblemData problem, List<TourData> tours, TourData candidateTour, long startTime) {
      bool candidateSaved = false;

      if (candidateTour.Length <= problem.Constraints[0].TargetValue) {
        bool processed = false;
        int tourCounterMax = tours.Count - 1;
        for (int tourCounter = tourCounterMax; tourCounter >= 0; --tourCounter) {
          bool scoreGt = false;
          bool scoreLt = false;
          for (int scoreCounter = 0; scoreCounter < problem.ScoreCount; scoreCounter++) {
            if (candidateTour.Scores[scoreCounter] > tours[tourCounter].Scores[scoreCounter]) {
              scoreGt = true;
            } else if (candidateTour.Scores[scoreCounter] < tours[tourCounter].Scores[scoreCounter]) {
              scoreLt = true;
            }
          }

          // If this tour equals an already existing one, skip it
          if (candidateTour.Path == tours[tourCounter].Path) {
            processed = true;
            break;
          }

          // If the tour is dominated by an existing, skip it
          if (scoreLt && !scoreGt) {
            processed = true;
            break;
          }

          // If the tour yields the same scores as an existing one, skip it
          if (!scoreLt && !scoreGt) {
            processed = true;
            break;
          }

          // If the tour dominates an existing one, eliminate the existing one
          if (!scoreLt && scoreGt) {
            tours.RemoveAt(tourCounter);
          }
        }
        if (!processed) {
          // The tour is neither dominated by an existing one nor does it
          // dominate an existing tour, so we assume it to be dominant
          candidateTour.TimeOfAdding = CalculateDifference(startTime, DateTime.Now.Ticks);
          tours.Add(candidateTour);
          candidateSaved = true;
        }
      }
      return candidateSaved;
    }

    #endregion

    #region Helper
    private static TimeSpan CalculateDifference(long startValue, long endValue) {
      return TimeSpan.FromTicks(endValue - startValue);
    }

    #endregion
  }

  #region Data Structures

  public class FeasiblePointVns {
    public int Index;
    public double Distance;
    public List<int> Scores;
  }

  public enum Constraint {
    Upperbound = 0,
    Lowerbound = 1,
    Maximize = 2,
    Minimize = 3
  }

  public class ConstraintData {
    public Constraint Type;           // 0...upper bound
    // 1...lower bound
    // 2...maximize
    // 3...minimize
    public int Target;            // Which score is affected by 
    // this constraint:
    // 0...Time/Distance
    // #...number of the score
    public double TargetValue;        // The value that should
    // (not) be reached
  }

  public class PointData {
    public double PositionX;    // X-Position of the point
    public double PositionY;    // Y-Position of the point

    public List<int> Score;   // List of scores for this point
    public List<int> Neighbors; // List of all points in the neighborhood
  }

  public class EdgeData {
    public double Distance;       // Length of this edge

    //vector<double> Pheromone;   // List of pheromones deposited
    // on this edge
  }

  public class TourData {
    public TourData() { }
    public TourData(TourData other) {
      Length = other.Length;
      PathSize = other.PathSize;
      Scores = new List<int>(other.Scores);
      Path = new List<int>(other.Path);
      TimeOfAdding = other.TimeOfAdding;
    }

    public double Length;       // Length of this tour
    public int PathSize;

    public List<int> Scores;    // Scores achieved by this tour
    public List<int> Path;      // Indices of the vertices visited by this tour
    public TimeSpan TimeOfAdding; // Time at which the solution was added
  }

  public class VnsData {
    public int NumberOfIterations;        // Number of iterations to run
    public int AcceptFrequency;           // How often "bad" solutions will be accepted
    public double TresholdValue;          // Ratio that allows the system to accept
    // Solutions that are not dominant but do
    // not differ more than this from the actual
    // solution
    public TourData InitialTour;          // Initially created solution tour
    public TourData ActualTour;           // Actually created neighborhood tour

    public List<TourData> DominantTours;  // List containing all dominant tours
  }

  public class ProblemData {
    public int StartingPoint;                 // Index of the starting point
    public int TerminusPoint;                 // Index of the ending point
    public int ConstraintCount;               // Number of constraints
    public int ScoreCount;                    // Number of scores in the problem
    public int PointCount;                    // Number of points in the problem
    public bool MatrixGiven;                  // true if the distance matrix is given

    public List<ConstraintData> Constraints;  // List of constraints
    public List<PointData> Points;            // List of points
    public List<List<EdgeData>> Edges;        // List of all edges (First index = origin, second index = target)

    public MersenneTwister MTR;               // Mersenne Twister PRNG

    //ACO_Data ACO;                           // All data needed by the ACO algorithm
    public VnsData VNS;                       // All data needed by the VNS algorithm
    public double FixedPenalty;               // User specified penalty for each visited vertex
  }

  #endregion
}