15.4.774  src/lib/graph/johnsons-all-pairs-shortest-paths-g.pkg

# johnsons-all-pairs-shortest-paths-g.pkg
#
# This is Johnson's algorithm for computing all pairs shortest paths.
# Good for sparse graphs.
# -- Allen Leung

# Compiled by:
#     src/lib/graph/graphs.lib

# See also:
#     src/lib/compiler/back/low/doc/latex/graphs.tex
#     src/lib/graph/test5.pkg



###           "Do you believe then that the sciences
###            would ever have arisen and become great
###            if there had not before hand been
###            magicians, alchemists, astrologers and wizards,
###            who thirsted and hungered after
###            abscondite and forbidden powers?"
###
###                         -- Friedrich Nietzsche 1886



stipulate
    package odg =  oop_digraph;                                         # oop_digraph                   is from   src/lib/graph/oop-digraph.pkg
    package mat =  rw_matrix;                                           # rw_matrix                     is from   src/lib/std/src/rw-matrix.pkg
    package rwv =  rw_vector;                                           # rw_vector                     is from   src/lib/std/src/rw-vector.pkg
    package ugv =  union_graph_view;                                    # union_graph_view              is from   src/lib/graph/uniongraph.pkg
herein
    generic package johnsons_all_pairs_shortest_paths_g (
        #
        num: Abelian_Group_With_Infinity                                # Abelian_Group_With_Infinity   is from   src/lib/graph/group.api
    )
    : (weak)     api {  include api All_Pairs_Shortest_Paths;           # All_Pairs_Shortest_Paths      is from   src/lib/graph/shortest-paths.api
                        exception NEGATIVE_CYCLE;
                     }
    {
        package num =  num;                                             # Export for client packages.

        stipulate
            package dsp =  dijkstras_single_source_shortest_paths_g(     num );
            package bsp =  bellman_fords_single_source_shortest_paths_g( num );
            #
            package meg =  digraph_by_adjacency_list_g( sparse_rw_vector );     # "meg" == "make_empty_graph" -- the only visible call in our "object-oriented" inteface (all other calls are via fns in the returned record).
        herein

            exception NEGATIVE_CYCLE =   bsp::NEGATIVE_CYCLE;

            fun all_pairs_shortest_paths
                  {   weight,   graph=> ggg as odg::DIGRAPH g:  odg::Digraph( N, E, G)   }
                =
                {   nnn   =  g.capacity ();
                    #
                    dist  =  mat::make_rw_matrix ((nnn, nnn), num::inf);
                    prior =  mat::make_rw_matrix ((nnn, nnn), -1);

                    exception EDGE  E;
                    exception NODE  N;
                    exception EMPTY;

                    fun arb_edge ()
                        = 
                        {   g.forall_edges (\\ (_, _, e) = raise exception EDGE e);
                            raise exception EMPTY;
                        }
                        except
                            EDGE e = e;

                    fun arb_node ()
                        = 
                        {   g.forall_nodes (\\ (_, n) = raise exception NODE n);
                            raise exception EMPTY;
                        }
                        except
                            NODE n = n;

                    {   e    = arb_edge ();
                        n    = arb_node ();

                        my  ggg' as odg::DIGRAPH g'
                            =
                            meg::make_empty_graph
                              {
                                graph_name          =>  "dummy source",         # Arbitrary client name for graph, for human-display purposes.
                                graph_info          =>  g.graph_info,           # Arbitrary client value to associate with graph.
                                expected_node_count  => 1                       # Hint for initial sizing of internal graph vectors.  This is not a hard limit. 
                              };

                        ggg''  = ugv::union_view (\\ (a, b) = a) (ggg, ggg');

                        my (+) = num::(+) ;
                        my (-) = num::(-) ;

                        s    = nnn;

                        g.forall_nodes
                           (\\ (v, _) = g'.add_edge (s, v, e));

                        g'.add_node (s, n);

                        fun weight'(u, v, e)
                            =
                            if (u == s)   num::zero;
                            else          weight (u, v, e);
                            fi;

                        my { dist=>h, ... }
                            =
                            dsp::single_source_shortest_paths
                                               { graph=>ggg'', s, weight=>weight'};

                        fun weight''(u, v, e)
                            =
                            weight (u, v, e) + ((rwv::get (h, u)) - (rwv::get (h, v)));

                        g.forall_nodes
                            (    \\ (u, _)
                                     =
                                     {    my { dist=>d, prior=>p }
                                               =
                                               bsp::single_source_shortest_paths
                                                   { graph=>ggg, s=>u, weight=>weight''};

                                          h_u =  rwv::get (h, u);

                                          g.forall_nodes
                                              (    \\ (v, _)
                                                       =
                                                       {   mat::set (dist,  (u, v), rwv::get (d, v) + rwv::get (h, v) - h_u);
                                                           mat::set (prior, (u, v), rwv::get (p, v));
                                                       }
                                              );
                                     }
                            );
                    }
                    except
                        EMPTY = ();

                    { dist, prior };
                }; 
        end;
    };
end;