add src/centrality

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: dodgr
 Title: Distances on Directed Graphs
-Version: 0.2.5.002
+Version: 0.2.5.003
 Authors@R: c(
     person("Mark", "Padgham", email="mark.padgham@email.com", role=c("aut", "cre")),
     person("Andreas", "Petutschnig", role="aut"),

R/RcppExports.R

@@ -1,6 +1,13 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+#' rcpp_centrality
+#'
+#' @noRd
+rcpp_centrality <- function(graph, vert_map_in, heap_type) {
+    .Call(`_dodgr_rcpp_centrality`, graph, vert_map_in, heap_type)
+}
+
 #' Make unordered_set of all new edge names
 #' @noRd
 NULL

codemeta.json

@@ -10,7 +10,7 @@
   "codeRepository": "https://github.com/ATFutures/dodgr",
   "issueTracker": "https://github.com/ATFutures/dodgr/issues",
   "license": "https://spdx.org/licenses/GPL-3.0",
-  "version": "0.2.5.2",
+  "version": "0.2.5.3",
   "programmingLanguage": {
     "@type": "ComputerLanguage",
     "name": "R",
@@ -348,14 +348,6 @@
         "url": "https://cran.r-project.org"
       },
       "sameAs": "https://CRAN.R-project.org/package=RcppParallel"
-    },
-    {
-      "@type": "SoftwareApplication",
-      "identifier": "https://sysreqs.r-hub.io/get/gnumake"
-    },
-    {
-      "@type": "SoftwareApplication",
-      "identifier": "https://sysreqs.r-hub.io/get/cxx11"
     }
   ],
   "releaseNotes": "https://github.com/ATFutures/dodgr/blob/master/NEWS.md",
@@ -400,5 +392,5 @@
     }
   ],
   "relatedLink": "https://CRAN.R-project.org/package=dodgr",
-  "fileSize": "638.249KB"
+  "fileSize": "639.187KB"
 }

src/Makevars

@@ -6,7 +6,7 @@ PKG_LIBS += $(shell ${R_HOME}/bin/Rscript -e "RcppParallel::RcppParallelLibs()")
 
 OBJ_HEAPS = heaps/bheap.o heaps/fheap.o heaps/heap23.o \
     heaps/radixheap.o heaps/triheap_ext.o heaps/triheap.o
-OBJ_SRC = dgraph.o pathfinders.o dodgr-to-sf.o flows.o fund-cycles.o \
+OBJ_SRC = centrality.o dgraph.o pathfinders.o dodgr-to-sf.o flows.o fund-cycles.o \
     graph-contract.o graph.o graph-sample.o RcppExports.o run_sp.o \
     sc-as-network.o sf-as-network.o turn_penalty.o
 OBJECTS = $(OBJ_HEAPS) $(OBJ_SRC)

src/Makevars.win

@@ -7,7 +7,7 @@ PKG_LIBS += $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" \
 
 OBJ_HEAPS = heaps/bheap.o heaps/fheap.o heaps/heap23.o \
     heaps/radixheap.o heaps/triheap_ext.o heaps/triheap.o
-OBJ_SRC = dgraph.o pathfinders.o dodgr-to-sf.o flows.o fund-cycles.o \
+OBJ_SRC = centrality.o dgraph.o pathfinders.o dodgr-to-sf.o flows.o fund-cycles.o \
     graph-contract.o graph.o graph-sample.o RcppExports.o run_sp.o \
     sc-as-network.o sf-as-network.o turn_penalty.o
 OBJECTS = $(OBJ_HEAPS) $(OBJ_SRC)

src/RcppExports.cpp

@@ -5,6 +5,19 @@
 
 using namespace Rcpp;
 
+// rcpp_centrality
+Rcpp::NumericVector rcpp_centrality(const Rcpp::DataFrame graph, const Rcpp::DataFrame vert_map_in, const std::string& heap_type);
+RcppExport SEXP _dodgr_rcpp_centrality(SEXP graphSEXP, SEXP vert_map_inSEXP, SEXP heap_typeSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const Rcpp::DataFrame >::type graph(graphSEXP);
+    Rcpp::traits::input_parameter< const Rcpp::DataFrame >::type vert_map_in(vert_map_inSEXP);
+    Rcpp::traits::input_parameter< const std::string& >::type heap_type(heap_typeSEXP);
+    rcpp_result_gen = Rcpp::wrap(rcpp_centrality(graph, vert_map_in, heap_type));
+    return rcpp_result_gen;
+END_RCPP
+}
 // rcpp_aggregate_to_sf
 Rcpp::List rcpp_aggregate_to_sf(const Rcpp::DataFrame& graph_full, const Rcpp::DataFrame& graph_contr, const Rcpp::DataFrame& edge_map);
 RcppExport SEXP _dodgr_rcpp_aggregate_to_sf(SEXP graph_fullSEXP, SEXP graph_contrSEXP, SEXP edge_mapSEXP) {
@@ -234,6 +247,7 @@ END_RCPP
 }
 
 static const R_CallMethodDef CallEntries[] = {
+    {"_dodgr_rcpp_centrality", (DL_FUNC) &_dodgr_rcpp_centrality, 3},
     {"_dodgr_rcpp_aggregate_to_sf", (DL_FUNC) &_dodgr_rcpp_aggregate_to_sf, 3},
     {"_dodgr_rcpp_aggregate_files", (DL_FUNC) &_dodgr_rcpp_aggregate_files, 2},
     {"_dodgr_rcpp_flows_aggregate_par", (DL_FUNC) &_dodgr_rcpp_flows_aggregate_par, 8},

src/centrality.cpp

@@ -0,0 +1,164 @@
+#include "centrality.h"
+#include "run_sp.h"
+#include "heaps/heap.h"
+
+#include <algorithm> // std::fill
+
+// # nocov start
+template <typename T>
+void inst_graph (std::shared_ptr<DGraph> g, unsigned int nedges,
+        const std::map <std::string, unsigned int>& vert_map,
+        const std::vector <std::string>& from,
+        const std::vector <std::string>& to,
+        const std::vector <T>& dist,
+        const std::vector <T>& wt)
+{
+    for (unsigned int i = 0; i < nedges; ++i)
+    {
+        unsigned int fromi = vert_map.at(from [i]);
+        unsigned int toi = vert_map.at(to [i]);
+        g->addNewEdge (fromi, toi, dist [i], wt [i]);
+    }
+}
+// # nocov end
+
+// igraph code starts at
+// https://github.com/igraph/igraph/blob/master/src/centrality.c#L1633
+
+void PF::PathFinder::Centrality (
+        std::vector <double>& w,
+        std::vector <double>& cent,
+        const unsigned int s)
+{
+    const DGraphEdge *edge;
+
+    const unsigned int n = m_graph->nVertices();
+    const std::vector <DGraphVertex>& vertices = m_graph->vertices();
+
+    std::deque <unsigned int> v_stack;
+
+    std::fill (w.begin (), w.end (), 0.0);
+    w [s] = 1.0;
+
+    m_heap->insert (s, -1.0);
+
+    std::vector <int> sigma (n, 0);
+    sigma [s] = 1;
+
+    std::vector <std::vector <int> > prev_arr (n);
+
+    while (m_heap->nItems() > 0) {
+        unsigned int v = m_heap->deleteMin();
+
+        v_stack.push_back (v);
+
+        edge = vertices [v].outHead;
+        while (edge) {
+
+            unsigned int et = edge->target;
+            double wt = w [v] + edge->dist;
+
+            std::unordered_set <int> w_set;
+            std::vector <int> w_vec;
+
+            if (w [et] == 0.0) // first connection to et
+            {
+                w_vec.resize (1);
+                w_vec [0] = v;
+                prev_arr [et] = w_vec;
+
+                sigma [et] = sigma [v];
+                w [et] = wt;
+                m_heap->insert (et, wt);
+            }  else if (wt < w [et])
+            {
+                w_vec.resize (1);
+                w_vec [0] = v;
+                prev_arr [et] = w_vec;
+
+                sigma [et] = sigma [v];
+                w [et] = wt;
+                m_heap->decreaseKey (et, wt);
+            } else if (wt == w [et])
+            {
+                w_vec = prev_arr [et];
+                w_vec.resize (w_vec.size () + 1);
+                w_vec [w_vec.size () - 1] = v;
+                prev_arr [et] = w_vec;
+
+                sigma [et] += sigma [v];
+            }
+
+            edge = edge->nextOut;
+        }
+    } // end while nItems > 0
+
+    // Then read from the stack and count centrality paths
+    std::vector <double> delta (n, 0.0);
+    while (!v_stack.empty ())
+    {
+        const unsigned int v = v_stack.back ();
+        v_stack.pop_back ();
+        std::vector <int> w_vec = prev_arr [v];
+        double tempd = (1.0 + delta [v]) / sigma [v];
+        for (auto ws: w_vec)
+        {
+            delta [ws] += sigma [ws] * tempd;
+        }
+        if (v != s)
+            cent [v] += delta [v];
+    }
+}
+
+
+//' rcpp_centrality
+//'
+//' @noRd
+// [[Rcpp::export]]
+Rcpp::NumericVector rcpp_centrality (const Rcpp::DataFrame graph,
+        const Rcpp::DataFrame vert_map_in,
+        const std::string& heap_type)
+{
+    std::vector <std::string> from = graph ["from"];
+    std::vector <std::string> to = graph ["to"];
+    std::vector <double> dist = graph ["d"];
+    std::vector <double> wt = graph ["d_weighted"];
+
+    unsigned int nedges = static_cast <unsigned int> (graph.nrow ());
+    std::map <std::string, unsigned int> vert_map;
+    std::vector <std::string> vert_map_id = vert_map_in ["vert"];
+    std::vector <unsigned int> vert_map_n = vert_map_in ["id"];
+    size_t nverts = run_sp::make_vert_map (vert_map_in, vert_map_id,
+            vert_map_n, vert_map);
+
+    std::shared_ptr <DGraph> g = std::make_shared <DGraph> (nverts);
+    inst_graph (g, nedges, vert_map, from, to, dist, wt);
+
+    Rcpp::NumericVector na_vec = Rcpp::NumericVector (nverts,
+            Rcpp::NumericVector::get_na ());
+    //Rcpp::NumericVector dout (static_cast <int> (nverts), na_vec.begin ());
+
+    // Create parallel worker
+    /*
+    OneIso one_iso (fromi, nverts, g, dlim, heap_type, dout);
+
+    RcppParallel::parallelFor (0, static_cast <size_t> (fromi.length ()),
+            one_iso);
+    */
+    std::vector <double> w (nverts), cent (nverts, 0.0);
+    for (unsigned int i = 0; i < nverts; i++)
+    {
+        std::shared_ptr<PF::PathFinder> pathfinder =
+            std::make_shared <PF::PathFinder> (nverts,
+                *run_sp::getHeapImpl(heap_type), g);
+
+        pathfinder->init (g); // specify the graph
+
+        pathfinder->Centrality (w, cent, i);
+    }
+
+    Rcpp::NumericVector dout = Rcpp::wrap (cent);
+
+    return (dout);
+}
+

src/centrality.h

@@ -0,0 +1,36 @@
+#pragma once
+
+#include <memory>
+#include <vector>
+#include <algorithm> // std::fill, std::reverse
+#include <iostream>
+#include <fstream>
+#include <deque> // used in centrality
+
+#include <Rcpp.h>
+// [[Rcpp::depends(RcppParallel)]]
+#include <RcppParallel.h>
+
+#include "pathfinders.h"
+#include "centrality.h"
+
+class DGraph;
+class PathFinder;
+
+struct by_w
+{
+    template <class T1, class T2>
+    bool operator () (const std::pair <T1, T2>& lhs,
+            const std::pair <T1, T2>& rhs)
+    {
+        if (lhs.second == rhs.second)
+            return (lhs.first > rhs.first);
+        else
+            return lhs.second > rhs.second;
+    }
+}; 
+
+Rcpp::NumericVector rcpp_centrality (const Rcpp::DataFrame graph,
+        const Rcpp::DataFrame vert_map_in,
+        const std::string& heap_type);
+

src/pathfinders.h

@@ -109,6 +109,10 @@ class PathFinder {
                 std::vector<double>& w,
                 std::vector<int>& prev,
                 unsigned int v0);
+        void Centrality (
+                std::vector <double>& w,
+                std::vector <double>& cent,
+                const unsigned int v0);
 
     private:
         Heap *m_heap;        // pointer: heap