Continuous Traveling Salesman Problem

tsp.adb

with Ada.Text_IO;
use Ada.Text_IO;

with Ada.Containers.Vectors;
with Ada.Numerics.Discrete_Random;
with Ada.Numerics.Float_Random;
with Ada.Numerics.Generic_Elementary_Functions;
with Ada.Numerics.Elementary_Functions;

procedure TSP is
  type Location is range 0 .. 100;
  type Point is record
    X: Location;
    Y: Location;
  end record;

  type Action is (Wiggle, Add, Delete);

  package Random is new Ada.Numerics.Discrete_Random(Location);
  package RandomAction is new Ada.Numerics.Discrete_Random(Action);

  Rng: Random.Generator;
  ARng: RandomAction.Generator;
  FRng: Ada.Numerics.Float_Random.Generator;

  package Point_Vectors is new Ada.Containers.Vectors(
    Index_Type => Natural,
    Element_Type => Point
  );

  procedure Put_Vector ( Vec : Point_Vectors.Vector ) is
  begin
    for E of Vec loop
      Put("[" & Location'Image(E.X) & ", " & Location'Image(E.Y) & "], ");
    end loop;
    New_Line;
  end Put_Vector;

  function Random_Point return Point is
    P : Point;
  begin
    P.X := Random.Random(Rng);
    P.Y := Random.Random(Rng);
    return P;
  end Random_Point;

  function Hypot(P1: Point; P2: Point) return Float is
    DX: Float;
    DY: Float;
  begin
    DX := Float(P1.X)-Float(P2.X);
    DY := Float(P1.Y)-Float(P2.Y);
    return Ada.Numerics.Elementary_Functions.sqrt(DX*DX + DY*DY);
  end Hypot;

  function Closest_Distance(P: Point; Locations: Point_Vectors.Vector) return Float is
    min: Float := Float'Last;
  begin
    for E of Locations loop
      if Hypot(P, E) < min then
        min := Hypot(P, E);
      end if;
    end loop;
    return min;
  end Closest_Distance;

  function Path_Error(Path: Point_Vectors.Vector; Locations: Point_Vectors.Vector) return Float is
    Error: Float := 0.0;
  begin
    for E of Locations loop
      Error := Error + Closest_distance(E, Path);
    end loop;
    return Error;
  end Path_Error;

  function Path_Length(Path: Point_Vectors.Vector) return Float is
    Prev : Point := Path.Last_Element;
    Length: Float := 0.0;
  begin
    for E of Path loop
      Length := Length + Hypot(Prev, E);
      Prev := E;
    end loop;
    return length;
  end Path_Length;

  function Temperature(ratio: Float) return Float is
  begin
    return 1.0-ratio;
  end Temperature;

  function Neighbour(Path: Point_Vectors.Vector) return Point_Vectors.Vector is
    New_Path: Point_Vectors.Vector := Path;
    Idx: Natural := Natural(Float(Path.Length) * Ada.Numerics.Float_Random.Random(FRng)) mod Natural(Path.Length);
    procedure Wiggle(P: in out Point) is
      DX: Integer := Integer(Ada.Numerics.Float_Random.Random(FRng)*10.0-5.0);
      DY: Integer := Integer(Ada.Numerics.Float_Random.Random(FRng)*10.0-5.0);
    begin
      if Integer(P.X) + DX < Integer(Location'Last) and Integer(P.X) + DX > Integer(Location'First) then
        P.X := Location(Integer(P.X) + DX);
      end if;
      if Integer(P.Y) + DY < Integer(Location'Last) and Integer(P.Y) + DY > Integer(Location'First) then
        P.Y := Location(Integer(P.Y) + DY);
      end if;
    end;
  begin
    case RandomAction.Random(ARng) is
      when Wiggle =>
        New_Path.Update_Element(Idx, Wiggle'Access);
      when Add =>
        New_Path.Insert(Idx, Random_Point);
      When Delete =>
        if Integer(New_Path.Length) > 2 then 
          New_Path.Delete(Idx);
        end if;
    end case;
    --New_Path.Replace_Element(Idx, Random_Point); -- TODO more smart
    return New_Path;
  end Neighbour;

  function Energy(
    Path: Point_Vectors.Vector;
    Locations: Point_Vectors.Vector;
    Temperature: Float
  ) return Float is
    Length: Float := Path_Length(Path);
    Error: Float := Path_Error(Path, Locations);
  begin
    return Length + Error/(1.0+Temperature);
  end Energy;

  function Probability(
    Energy: Float;
    New_Energy: Float;
    Temperature: Float) return Float is
  begin
    if New_Energy < Energy then
      return 1.0;
    else
      return Temperature; -- TODO more smart
    end if;
  end Probability;

  Locations : Point_Vectors.Vector;
  Path : Point_Vectors.Vector;
begin
  Random.Reset(Rng);
  RandomAction.Reset(ARng);
  Ada.Numerics.Float_Random.Reset(FRng);

  for I in 1 .. 20 loop
    Locations.append(Random_Point);
  end loop;

  for I in 1 .. 3 loop
    Path.append(Random_Point);
  end loop;

  for I in 0 .. 1_000_000 loop
    declare
      T: Float := Temperature(Float(I)/1_000_000.0);
      New_Path: Point_Vectors.Vector := Neighbour(Path);
      E: Float := Energy(Path, Locations, T);
      Ep: Float := Energy(New_Path, Locations, T);
    begin
      if Probability(E, Ep, T) > Ada.Numerics.Float_Random.Random(FRng) then
        Path := New_Path;
      end if;
      --Put_Line("T: " & Float'Image(T) & " E: " & Float'Image(E) & " E': " & Float'Image(Ep) & " P: " & Float'Image(Probability(E, Ep, T)) & " Length: " & Float'Image(Path_Length(Path)) & " Error:  " & Float'Image(Path_Error(Path, Locations)));
    end;
  end loop;
  Put_Line("loc = list(zip(");
  Put_Vector(Locations);
  Put_Line("))");
  Put_Line("path = list(zip(");
  Put_Vector(Path);
  Put_Line("))");
end TSP;