/*-----------------------------------------------------------------------------

Copyright (C) 2020

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA

This program 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 2 of the License, or
(at your option) any later version.

This program 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 this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

-------------------------------------------------------------------------------

Function      dijkstra - implements Dijkstra's algorithm

Syntax        #include "libscl.h"

              realmat dijkstra(INTEGER nodes, INTEGER edges, INTEGER start, 
                const std::vector<edge>& graph);

              realmat dijkstra(INTEGER nodes, INTEGER edges, INTEGER start, 
                const std::vector<edge>& graph, std::vector<intvec>& route);

              Edge is a struct declared in libscl.h as follows

              struct edge {
                INTEGER first;
                INTEGER second;
                REAL distance;
                edge();
                edge(INTEGER f, INTEGER s, REAL d);
                edge(const edge& e);
                edge& operator=(const edge& e);
              };

Prototype in  libscl.h

Description   Dijkstra's algorithm finds the shortest path to each node
              in a graph from the start node. If there is no route from
              start to node, REAL_MAX is the returned value for that
              node. The intvec at route[u] is the path from start to node
              u with start as its first element. On return the size of
              route is nodes+1; route is indexed as route[i] for = 1,..,nodes.

Remarks       In graph it does not matter if an edge is entered twice as
              edge(1,2,2.5) and edge(2,1,2.5) or just once as either. If
              entered twice and the distances do not agree it will be
              the second that occurs in graph that gets used.

Return value  A realmat dimensioned (nodes,1) containing the distance
              from start to each node. Start is one of the nodes and 
              will have a distance to itself of zero. If there is no
              route from start to node, REAL_MAX is the returned
              value for that node. 

References:   http://zobayer.blogspot.com/2009/12/dijkstras-algorithm-in-c.html
              https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm
              Dijkstra, E. W. (1959), "A Note on Two Problems in Connexion 
              with Graphs," Numerische Mathematik 1 269--271

              
Functions     Library: vector, fabs, sqrt
called        libscl:  intvec, realmat

------------------------------------------------------------------------------*/

#include "libscl.h"
#include <queue>

using namespace std;

namespace scl {

typedef pair<INTEGER,REAL> IRpair;

struct compIRpair {
   bool operator() (const IRpair& a, const IRpair& b) {
       return a.second > b.second;
   }
};

realmat dijkstra(INTEGER nodes, INTEGER edges, INTEGER start, 
        const vector<edge>& graph, vector<intvec>& path)
{
  INTEGER u, v, sz;
  REAL w;

  size_t size = edges > nodes ? edges : nodes;

  ++size;

  if (size > INT_MAX) {
   error("Error, dijkstra, either too many edges or too many nodes");
  }
  
  priority_queue< IRpair, vector< IRpair >, compIRpair > Q;
  vector< vector< IRpair > > G(size);
  realmat D(nodes,1);
  vector<bool> F(size);
  path.resize(nodes+1);

  --size;

  for (size_t i=1; i<=size; ++i) F[i] = false;

  for (vector<edge>::const_iterator i=graph.begin(); i!=graph.end(); i++) {
    u = i->first;
    v = i->second;
    w = i->distance;
    G[u].push_back(IRpair(v, w));
    G[v].push_back(IRpair(u, w)); // for undirected
  }

  // initialize distance vector
  for (INTEGER i=1; i<=nodes; i++) D[i] = DBL_MAX;
  D[start] = 0.0;
  Q.push(IRpair(start,0.0));

  // dijkstra
  while (!Q.empty()) {
    u = Q.top().first; 
    Q.pop();
    if (F[u]) continue;
    sz = G[u].size(); 
    for (INTEGER i=0; i<sz; i++) {
      v = G[u][i].first;
      w = G[u][i].second;
      if (!F[v] && D[u] + w < D[v]) {
        D[v] = D[u] + w;
        //D.check(v) = D.check(u) + w;
        Q.push(IRpair(v,D[v])); 
        intvec prev = path[u];
        prev.push_back(u);
        path[v] = prev;
      }
    }
    F[u] = true; // done with u 
  }

  for (INTEGER v=1; v<=nodes; ++v) {
    if (path[v].size() != 0) path[v].push_back(v);
  }
  path[start].resize(2,start);

  return D; 
}

// Below is above with eight lines involving path commented out

realmat dijkstra(INTEGER nodes, INTEGER edges, INTEGER start, 
        const vector<edge>& graph)
{
  INTEGER u, v, sz;
  REAL w;

  size_t size = edges > nodes ? edges : nodes;

  ++size;

  if (size > INT_MAX) {
   error("Error, dijkstra, either too many edges or too many nodes");
  }
  
  priority_queue< IRpair, vector< IRpair >, compIRpair > Q;
  vector< vector< IRpair > > G(size);
  realmat D(nodes,1);
  vector <bool> F(size);

  --size;

  /*
  path.resize(nodes+1);
  */

  for (size_t i=1; i<=size; ++i) F[i] = false;

  for (vector<edge>::const_iterator i=graph.begin(); i!=graph.end(); i++) {
    u = i->first;
    v = i->second;
    w = i->distance;
    G[u].push_back(IRpair(v, w));
    G[v].push_back(IRpair(u, w)); // for undirected
  }

  // initialize distance vector
  for (INTEGER i=1; i<=nodes; i++) D[i] = DBL_MAX;
  D[start] = 0.0;
  Q.push(IRpair(start,0.0));

  // dijkstra
  while (!Q.empty()) {
    u = Q.top().first; 
    Q.pop();
    if (F[u]) continue;
    sz = G[u].size(); 
    for (INTEGER i=0; i<sz; i++) {
      v = G[u][i].first;
      w = G[u][i].second;
      if (!F[v] && D[u] + w < D[v]) {
        D[v] = D[u] + w;
        //D.check(v) = D.check(u) + w;
        Q.push(IRpair(v,D[v])); 
        /*
        intvec prev = path[u];
        prev.push_back(u);
        path[v] = prev;
        */
      }
    }
    F[u] = true; // done with u 
  }

  /*
  for (INTEGER v=1; v<=nodes; ++v) {
    if (path[v].size() != 0) path[v].push_back(v);
  }
  path[start].resize(2,start);
  */

  return D; 
}
}
