/* 

  Floyd-Warshall all shortest path algorithm in Picat.

  http://en.wikipedia.org/wiki/Floyd%E2%80%93Warshall_algorithm

  This program implements:
    - Floyd-Warshall's all shortest path algorithm with path reconstruction
    - three variants of detecting cycles in the graph

  This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/

import util.

main => go.


go =>
  % Example from the Wikipedia article
  Edges = [ % from, to, weight
            [1,3,-2],
            [2,1,4],
            [2,3,3],
            [3,4,2],
            [4,2,-1]

            % , [1,2,-14] % make a cycle

          ],
  println(edges=Edges),
  Nodes = [ [I,J] : [I,J,_] in Edges].flatten().sort_remove_dups(),
  println(nodes=Nodes),
  Inf=1000,
  floydWarshall(Nodes,Edges,Inf, Dist,Next),

  println("Dist:"),
  foreach(Row in Dist)
    println(Row.to_list())
  end,
  DistDiagonal = Dist.diagonal1(),
  println(diagonalDist=DistDiagonal),

  % 
  % Cycle detection:
  % If there are any != 0 in the diagonal of the Dist matrix then it
  % contain cycles.
  % (See go2/0 for some more methods to detect cycles.)
  % 
  if [D : D in DistDiagonal, D != 0].length > 0 then
    println("Cycle detected (via Dist diagonal)"),
    println(dist_cycle_detected=DistDiagonal),
    % println(edges_involved_in_cycle=[I : I in 1..DistDiagonal.length, DistDiagonal[I] != 0]),
    % a more finegrained analysis of the cycle
    println(tc_cycles=detect_cycle_tc(Next))
  else 
    println("No cycle detected."),
    foreach(I in Nodes, J in Nodes)
       if Dist[I,J] == Inf then
         printf("There is no path between %d and %d\n", I,J)
       else 
         get_path(Next, I,J) = Path,
         println(path=[I,J]=Path)
       end
    end
  end, 


  nl.

%
% some more tests
%
go2 =>
  % Example from the Wikipedia article
  Edges = [ % from, to, weight
            % [1,3,-2],
            % [2,1,4],
            %  [2,3,3],
            % [3,4,2],
            % [4,2,-1]

            % , [5,4,-13]
            % , [1,2,-14]
            % , [3,1,-4]
            % , [4,3,-11]
            % , [3,5,1]
            % , [5,3,-11]
      

            % a non-connected graph
            % [1,2,2],
            % [1,3,-4],
            % [3,2,1],

            % [4,5,2],
            % [5,6,3],
            % [4,6,-3] % make it cyclic
            % , [6,4,-3]


            % Graph from http://www.geeksforgeeks.org/dynamic-programming-set-16-floyd-warshall-algorithm/
            % [1,2,5],
            % [1,4,10],
            % [2,3,3],
            % [3,4,1]


            % Graph from http://en.algoritmy.net/article/45708/Floyd-Warshall-algorithm
            % [1,2,5],
            % [1,4,2],
            % [2,3,2],
            % [3,1,3],
            % [3,5,7],
            % [4,3,4],
            % [4,5,1],
            % [5,1,1],
            % [5,2,3]
            % % , [3,2,-3] % make cycle


            % Graph from http://faculty.ycp.edu/~dbabcock/PastCourses/cs360/lectures/lecture23.html
            % [1,3,6],
            % [1,4,3],
            % [2,1,3],
            % [3,4,2],
            % [4,2,1],
            % [4,3,1],
            % [5,2,4],
            % [5,4,2]

            % Graph from http://www.cs.utexas.edu/users/djimenez/utsa/cs3343/lecture19.html
            % [1,2,1], 
            % [1,3,1], 
            % [2,4,1], 
            % [2,5,1], 
            % [3,1,1], 
            % [3,6,1], 
            % [4,6,1], 
            % [4,3,1], 
            % [6,5,1]

            % Graph for http://www.puzzlor.com/2010-02_PlanetColonization.html
            [1,2,1],
            [1,8,1],
            [1,5,1],
            [2,1,1],
            [2,10,1],
            [2,3,1],
            [3,4,1],
            [3,12,1],
            [3,2,1],
            [4,3,1],
            [4,14,1],
            [4,5,1],
            [5,4,1],
            [5,1,1],
            [5,6,1],
            [6,7,1],
            [6,5,1],
            [6,15,1],
            [7,6,1],
            [7,8,1],
            [7,17,1],
            [8,7,1],
            [8,1,1],
            [8,9,1],
            [9,8,1],
            [9,10,1],
            [9,20,1],
            [10,9,1],
            [10,2,1],
            [10,11,1],
            [11,10,1],
            [11,19,1],
            [11,12,1],
            [12,11,1],
            [12,3,1],
            [12,13,1],
            [13,12,1],
            [13,14,1],
            [13,18,1],
            [14,13,1],
            [14,4,1],
            [14,15,1],
            [15,14,1],
            [15,16,1],
            [15,6,1],
            [16,15,1],
            [16,18,1],
            [16,17,1],
            [17,16,1],
            [17,20,1],
            [17,7,1],
            [18,16,1],
            [18,19,1],
            [18,13,1],
            [19,18,1],
            [19,11,1],
            [19,20,1],
            [20,19,1],
            [20,17,1],
            [20,9,1]

           % % From https://rosettacode.org/wiki/Floyd-Warshall_algorithm
           % [1, 3, -2], 
           % [2, 1, 4], 
           % [2, 3, 3], 
           % [3, 4, 2], 
           % [4, 2, -1]

          ],

  % Edges2 = (Edges ++ [[J,I,-W] : [I,J,W] in Edges]).sort(),
  % Edges := Edges2,
  println(edges=Edges),
  Nodes = [ [I,J] : [I,J,_] in Edges].flatten().sort_remove_dups(),
  Len = Nodes.length,
  println(nodes=Nodes),
  Inf=1000,
  floydWarshall(Nodes,Edges,Inf, Dist,Next),
  println("Dist:"),
  foreach(Row in Dist)
    println(Row.to_list())
  end,
  DistDiagonal = Dist.diagonal1(),
  println(diagonalDist=DistDiagonal),

  % cycle detection: Dist diagonal  
  if [D : D in DistDiagonal, D != 0].length > 0 then
    println("Cycle detected (Dist diagonal)"),
    println(dist_cycle_detected=DistDiagonal),
    println(edges_involved_in_cycle=[I : I in 1..DistDiagonal.length, DistDiagonal[I] != 0])
  end, 

  println("\nNext:"),
  foreach(Row in Next)
    % println(Row.to_list())
    foreach(R in Row)
       printf("%3d", R)
    end,
    nl
  end,
  NextDiagonal = Next.diagonal1(),
  println(diagonalNext=NextDiagonal),

  % cycle detection: transitive closure of the next diagonal
  TC_Cycles = detect_cycle_tc(Next),
  % println(tc_cycles=TC_Cycles),
  % Cycles = detect_cycle(Next),
  % println(cycles=Cycles),
  if TC_Cycles != [] then
    println("Cycle detected!"),
    println(tc_cycles=TC_Cycles)
  else 
    foreach(I in Nodes, J in Nodes, I != J)
       % println([I,J,dist=Dist[I,J]]),
       if Dist[I,J] == Inf then
         printf("There is no path between %d and %d\n", I,J)
       else 
         get_path(Next, I,J) = Path,
         % println(path=[I,J]=Path)
         printf("%3d,%3d  (dist %2d) %w \n", I,J,Dist[I,J],Path)
       end
    end
  end,

  % another way to detect cycles in the Next diagonal
  Cycles2 = detect_cycles(Next),
  if Cycles2 != [] then
     println("Cycles detected (detect_cycles(Next))"),
     println(cycles2=Cycles2)
  end,

  % checking transitive closure
  EdgesTC = [[I,J,1] : [I,J,_] in Edges],
  floydWarshall(Nodes,EdgesTC,Inf, DistTC,_NextTC),
  IsTC = true,
  foreach(I in 1..Len, J in 1..Len)
    if DistTC[I,J] > Len then
      % println("Not transitive closure"=[I,J]),
      IsTC := false
    end
  end,
  if IsTC then
    println("Is a transitive closure")
  else 
    println("Not a transitive closure")
  end,
  nl.

%
% Floyd-Warshall with path reconstruction
%
% let dist be a |V| × |V| array of minimum distances initialized to ∞ (infinity)
% let next be a |V| × |V| array of vertex indices initialized to null
floydWarshall(Nodes,Edges,Inf, Dist,Next) =>
  % println($floydWarshall(nodes=Nodes,edges=Edges, Dist,Next)),
  Len = Nodes.length,
  Dist = new_array(Len,Len), bind_vars(Dist,Inf),
  % note: cannot use bind_vars(Dist,_) since then all share the same
  %       variable
  Next = new_array(Len,Len), 
  foreach(I in 1..Len, J in 1..Len)
    Next[I,J] := _
  end,
  foreach([U,V,Weight] in Edges)
    Dist[U,V] := Weight,  % the weight of the edge (u,v)
    Next[U,V] := V
  end,
  % for cycle detection
  foreach(I in Nodes)
    Dist[I,I] := 0,
    Next[I,I] := I
  end,
  foreach(K in 1..Len) % standard Floyd-Warshall implementation
     foreach(I in 1..Len)
       foreach(J in 1..Len)
          if Dist[K,J] != Inf, Dist[I,K] != Inf,
             Dist[I,K] + Dist[K,J] < Dist[I,J] 
           then
             Dist[I,J] := Dist[I,K] + Dist[K,J],
             Next[I,J] := Next[I,K]
          end
       end
     end
  end.

% path extraction from the Next structure
get_path(Next, U,V) = Path =>
   if var(Next[U,V]) then
     Path = []
   else 
     Path = [U],
     while (U != V, Next[U,V] != U)
         U := Next[U,V],
         Path := Path ++ [U]
     end
   end.

%
% Detect all cycles in the Next diagonal.
%
detect_cycles(Next) = Cycles => 
  NextDiagonal = Next.diagonal1(),
  Nodes = 1..NextDiagonal.length,
  Cycles = [],
  foreach(I in Nodes)
    D = NextDiagonal[I],
    Seen = new_map([I=1]),
    ThisCycle = [I,D],
    Cycle = [],
    Path = [I,D],
    while(D != I, Cycle = []) 
      Path := Path ++ [D],
      D := NextDiagonal[D],
      ThisCycle := ThisCycle ++ [D],
      Seen.put(D,1),
      if Seen.has_key(I) then
        Cycle := ThisCycle ++ [I]
      end
    end,
    if Cycle != [] then
      Cycles := Cycles ++ [Cycle]
    end
  end.


%
% (Adapted from apl_util.pi)
%
tc(List, Start) = TC => 
   TC = [Start],
   Next = List[Start],
   while (not member(Next,TC))
       TC := TC ++ [Next],
       % added to show the cycle
       if List[Next] == Start then
         TC := TC ++ [Start]
       end,
       Next := List[Next]
   end.

% (from apl_util.pi)
%
% all transitive closures for a list L (Start in 1..L.length)
%
tc_all(List) = [tc(List,I) : I in 1..List.length].   

% detect cycles via transitive cycles
detect_cycle_tc(Next) = Cycles =>
  TC_All = tc_all(Next.diagonal1()),
  Cycles = [],
  foreach({I,TC} in zip(1..length(TC_All),TC_All))
    if length(TC) > 1, I == TC.last() then
      Cycles := Cycles ++ [TC]
    end
  end.