Defining pow functions for adjacency and transition matrices

include/CXXGraph/Graph/Algorithm/Pow_impl.hpp

@@ -0,0 +1,189 @@
+/***********************************************************/
+/***      ______  ____  ______                 _         ***/
+/***     / ___\ \/ /\ \/ / ___|_ __ __ _ _ __ | |__	     ***/
+/***    | |    \  /  \  / |  _| '__/ _` | '_ \| '_ \	 ***/
+/***    | |___ /  \  /  \ |_| | | | (_| | |_) | | | |    ***/
+/***     \____/_/\_\/_/\_\____|_|  \__,_| .__/|_| |_|    ***/
+/***                                    |_|			     ***/
+/***********************************************************/
+/***     Header-Only C++ Library for Graph			     ***/
+/***	 Representation and Algorithms				     ***/
+/***********************************************************/
+/***     Author: ZigRazor ***/
+/***	 E-Mail: [email protected] 				     ***/
+/***********************************************************/
+/***	 Collaboration: ----------- 				     ***/
+/***********************************************************/
+/***	 License: AGPL v3.0 ***/
+/***********************************************************/
+
+#ifndef __CXXGRAPH_POW_IMPL_H__
+#define __CXXGRAPH_POW_IMPL_H__
+
+#pragma once
+
+#include "CXXGraph/Graph/Graph_decl.h"
+
+template <typename T>
+std::vector<std::vector<T>> matMult(std::vector<std::vector<T>> &a, std::vector<std::vector<T>> &b) {
+    static_assert(std::is_same<T, unsigned long long>::value || std::is_same<T, double>::value,
+                  "Type T must be either unsigned long long or double");
+
+    int n = (int)a.size(); // N x N matrix
+    std::vector<std::vector<T>> res(n, std::vector<T>(n, 0));
+
+    // O(n^3) matrix multiplication
+    for (int i = 0; i < n; i++) {
+        for (int j = 0; j < n; j++) {
+            for (int k = 0; k < n; k++) {
+                res[i][j] += a[i][k] * b[k][j];
+            }
+        }
+    }
+
+    return res;
+}
+
+template <typename T>
+std::vector<std::vector<T>> exponentiation(std::vector<std::vector<T>> &mat, unsigned int k) {
+    // typechecking
+    static_assert(std::is_same<T, unsigned long long>::value || std::is_same<T, double>::value,
+                  "Type T must be either unsigned long long or double");
+
+    int n = (int)mat.size();
+    std::vector<std::vector<T>> res(n, std::vector<T>(n, 0)); 
+
+    // build identity matrix    
+    for (int i = 0; i < n; i++) res[i][i] = 1;
+
+    // fast exponentation!
+    while (k) {
+        if (k & 1) res = matMult(res, mat);
+        mat = matMult(mat, mat);
+        k >>= 1;
+    }
+
+    return res;
+}
+
+namespace CXXGraph {
+
+template <typename T>
+const PowAdjResult Graph<T>::powAdjMatrix(unsigned int k) const {
+    PowAdjResult result;
+    result.success = false;
+    result.errorMessage = "";
+    std::unordered_map<std::pair<std::string, std::string>, unsigned long long,
+                        CXXGraph::pair_hash> powAdj;
+
+    // get graph sets
+    const auto &nodeSet = Graph<T>::getNodeSet();
+    const auto &edgeSet = Graph<T>::getEdgeSet();
+
+    int n = (int)nodeSet.size();
+
+    // convert back and forth between user id's and index in temporary adj matrix
+    std::unordered_map<std::string, int> userIdToIdx;
+    std::unordered_map<int, std::string> idxToUserId;
+
+    // adj int will store the temporary (integer) adjacency matrix
+    std::vector<std::vector<unsigned long long>> tempIntAdj(n, std::vector<unsigned long long>(n, 0));
+
+    // write forwards and backwards mapping
+    int i = 0;
+    for (const auto &node : nodeSet) { 
+        userIdToIdx[node->getUserId()] = i;
+        idxToUserId[i] = node->getUserId();
+        i++;
+    }
+
+    // populated temporary adjacency matrix w/ edges
+    // can handle both directed and undirected
+    for (const auto &e : edgeSet) {
+        const auto &edge = e->getNodePair();
+
+        const auto firstId = edge.first->getUserId();
+        const auto secondId = edge.second->getUserId();
+
+        // if undirected, add both sides
+        if (e->isDirected() == false)
+            tempIntAdj[userIdToIdx[secondId]][userIdToIdx[firstId]] = 1;
+        tempIntAdj[userIdToIdx[firstId]][userIdToIdx[secondId]] = 1;
+    } 
+
+    // calculate the power matrix
+    auto powerMatrix = exponentiation(tempIntAdj, k);
+
+    for (int i = 0; i < n; i++) {
+        for (int j = 0; j < n; j++) {
+            auto pr = std::make_pair(idxToUserId[i], idxToUserId[j]);
+            powAdj[pr] = powerMatrix[i][j];
+        }
+    }
+
+    result.result = std::move(powAdj);
+    result.success = true;
+
+    return result;
+}
+
+template <typename T>
+const PowTransResult Graph<T>::powTransitionMatrix(unsigned int k) const {
+    PowTransResult result;
+    result.success = false;
+    result.errorMessage = "";
+    std::unordered_map<std::pair<std::string, std::string>, double,
+                        CXXGraph::pair_hash> powTrans;
+
+    // get graph sets
+    const auto &nodeSet = Graph<T>::getNodeSet();
+    const auto &edgeSet = Graph<T>::getEdgeSet();
+
+    int n = (int)nodeSet.size();
+
+    const auto &transMatrix = Graph<T>::getTransitionMatrix();
+    std::unordered_map<std::string, int> userIdToIdx;
+    std::unordered_map<int, std::string> idxToUserId;
+
+    std::vector<std::vector<double>> tempDoubleTrans(n, std::vector<double>(n, 0));
+
+    // get a map between index in adj matrix
+    // and userId
+    int i = 0;
+    for (const auto &node : nodeSet) { 
+        userIdToIdx[node->getUserId()] = i;
+        idxToUserId[i] = node->getUserId();
+        i++;
+    }
+
+    // given transition matrix, convert it to 
+    // stochastic matrix
+    for (const auto &it : *transMatrix) {
+        const auto f = it.first;
+        const auto idx = userIdToIdx[f->getUserId()];
+
+        for (const auto &e : it.second) {
+            const auto idx2 = userIdToIdx[e.first->getUserId()];
+            tempDoubleTrans[idx][idx2] = e.second;
+        }
+    }
+
+    // exponentiate stochastic matrix
+    auto powerTransMatrix = exponentiation(tempDoubleTrans, k);
+
+    // turn back into a map between nodes
+    for (int i = 0; i < n; i++) {
+        for (int j = 0; j < n; j++) {
+            auto pr = std::make_pair(idxToUserId[i], idxToUserId[j]);
+            powTrans[pr] = powerTransMatrix[i][j];
+        }
+    }
+
+    result.result = std::move(powTrans);
+    result.success = true;
+
+    return result;
+}
+}
+
+#endif //_CXXGRAPH_POW_IMPL_H__

include/CXXGraph/Graph/Graph.h

@@ -37,6 +37,7 @@
 #include "CXXGraph/Graph/Algorithm/Kahn_impl.hpp"
 #include "CXXGraph/Graph/Algorithm/Kosaraju_impl.hpp"
 #include "CXXGraph/Graph/Algorithm/Kruskal_impl.hpp"
+#include "CXXGraph/Graph/Algorithm/Pow_impl.hpp"
 #include "CXXGraph/Graph/Algorithm/Prim_impl.hpp"
 #include "CXXGraph/Graph/Algorithm/Tarjan_impl.hpp"
 #include "CXXGraph/Graph/Algorithm/TopologicalSort_impl.hpp"

include/CXXGraph/Graph/Graph_decl.h

@@ -534,6 +534,22 @@ class Graph {
    */
   virtual std::shared_ptr<std::vector<Node<T>>> eulerianPath() const;
 
+  /**
+   * @brief This function raises the adjacency matrix to some k. 
+   * Best used for counting number of k-length walks from i to j.
+   * @param k value by which to raise matrix
+   * @return (success, errorMessage, matrix): where matrix is equivalent to A^k
+   */
+  virtual const PowAdjResult powAdjMatrix(unsigned int k) const;
+
+  /**
+   * @brief This function raises a transition matrix to some k. 
+   * Best used for finding equilibrium. 
+   * @param k value by which to raise matrix
+   * @return (success, errorMessage, matrix): where matrix is equivalent to S^k
+   */
+  virtual const PowTransResult powTransitionMatrix(unsigned int k) const;
+
   /**
    * @brief Function runs the dijkstra algorithm for some source node and
    * target node in the graph and returns the shortest distance of target

include/CXXGraph/Utility/Typedef.hpp

@@ -237,6 +237,26 @@ struct BestFirstSearchResult_struct {
 template <typename T>
 using BestFirstSearchResult = BestFirstSearchResult_struct<T>;
 
+/// Struct that contains the results from an adjacency matrix exponentiation
+/// results
+struct PowAdjResult_struct {
+  bool success = false;
+  std::string errorMessage = "";
+  std::unordered_map<std::pair<std::string, std::string>, unsigned long long, pair_hash>
+      result = {};
+};
+typedef PowAdjResult_struct PowAdjResult;
+
+/// Struct that contains the results from a transition matrix exponentiation
+/// results
+struct PowTransResult_struct {
+  bool success = false;
+  std::string errorMessage = "";
+  std::unordered_map<std::pair<std::string, std::string>, double, pair_hash>
+      result = {};
+};
+typedef PowTransResult_struct PowTransResult;
+
 ///////////////////////////////////////////////////////////////////////////////////
 // Using Definition
 // ///////////////////////////////////////////////////////////////

test/PowTest.cpp

@@ -0,0 +1,134 @@
+
+#include <memory>
+
+#include "CXXGraph/CXXGraph.hpp"
+#include "gtest/gtest.h"
+
+// include it to check that the static asserts don't fail
+#include "TypeTraitsTest.hpp"
+
+// Smart pointers alias
+template <typename T>
+using unique = std::unique_ptr<T>;
+template <typename T>
+using shared = std::shared_ptr<T>;
+
+using std::make_shared;
+using std::make_unique;
+
+// https://eng.libretexts.org/Bookshelves/Computer_Science/Programming_and_Computation_Fundamentals/Mathematics_for_Computer_Science_(Lehman_Leighton_and_Meyer)/02%3A_Structures/09%3A_Directed_graphs_and_Partial_Orders/9.03%3A_Adjacency_Matrices
+TEST(PowAdjTest, libre_texts) {
+  CXXGraph::Node<int> a("a", 1);
+  CXXGraph::Node<int> b("b", 1);
+  CXXGraph::Node<int> c("c", 1);
+  CXXGraph::Node<int> d("d", 1);
+
+  CXXGraph::DirectedEdge<int> e1(0, a, b);
+  CXXGraph::DirectedEdge<int> e2(1, a, d);
+  CXXGraph::DirectedEdge<int> e3(2, b, c);
+  CXXGraph::DirectedEdge<int> e4(3, b, d);
+  CXXGraph::DirectedEdge<int> e5(4, c, b);
+  CXXGraph::DirectedEdge<int> e6(5, d, c);
+
+  CXXGraph::Graph<int> graph;
+  graph.addEdge(&e1);
+  graph.addEdge(&e2);
+  graph.addEdge(&e3);
+  graph.addEdge(&e4);
+  graph.addEdge(&e5);
+  graph.addEdge(&e6);
+
+  auto res = graph.powAdjMatrix(2);
+
+  ASSERT_TRUE(res.success);
+  ASSERT_TRUE(res.result[std::make_pair("a", "c")] == 2);
+  ASSERT_TRUE(res.result[std::make_pair("a", "d")] == 1);
+  ASSERT_TRUE(res.result[std::make_pair("b", "b")] == 1);
+  ASSERT_TRUE(res.result[std::make_pair("b", "c")] == 1);
+  ASSERT_TRUE(res.result[std::make_pair("c", "c")] == 1);
+  ASSERT_TRUE(res.result[std::make_pair("c", "d")] == 1);
+  ASSERT_TRUE(res.result[std::make_pair("d", "b")] == 1);
+  ASSERT_TRUE(res.result[std::make_pair("c", "a")] == 0);
+  ASSERT_TRUE(res.result[std::make_pair("d", "c")] == 0);
+}
+
+// https://docs.dgl.ai/generated/dgl.khop_adj.html#dgl.khop_adj
+TEST(PowAdjTest, dgl) {
+  CXXGraph::Node<int> a("a", 1);
+  CXXGraph::Node<int> b("b", 1);
+  CXXGraph::Node<int> c("c", 1);
+  CXXGraph::Node<int> d("d", 1);
+  CXXGraph::Node<int> e("e", 1);
+
+  CXXGraph::DirectedEdge<int> e1(0, a, a);
+  CXXGraph::DirectedEdge<int> e2(1, b, b);
+  CXXGraph::DirectedEdge<int> e3(2, c, c);
+  CXXGraph::DirectedEdge<int> e4(3, d, d);
+  CXXGraph::DirectedEdge<int> e5(4, e, e);
+  CXXGraph::DirectedEdge<int> e6(5, a, b);
+  CXXGraph::DirectedEdge<int> e7(6, b, c);
+  CXXGraph::DirectedEdge<int> e8(7, c, d);
+  CXXGraph::DirectedEdge<int> e9(8, d, e);
+  CXXGraph::DirectedEdge<int> e10(9, d, a);
+
+  CXXGraph::Graph<int> graph;
+  graph.addEdge(&e1);
+  graph.addEdge(&e2);
+  graph.addEdge(&e3);
+  graph.addEdge(&e4);
+  graph.addEdge(&e5);
+  graph.addEdge(&e6);
+  graph.addEdge(&e7);
+  graph.addEdge(&e8);
+  graph.addEdge(&e9);
+  graph.addEdge(&e10);
+
+  auto res = graph.powAdjMatrix(3);
+
+  ASSERT_TRUE(res.success);
+  ASSERT_TRUE(res.result[std::make_pair("a", "a")] == 1);
+  ASSERT_TRUE(res.result[std::make_pair("b", "b")] == 1);
+  ASSERT_TRUE(res.result[std::make_pair("c", "c")] == 1);
+  ASSERT_TRUE(res.result[std::make_pair("d", "d")] == 1);
+  ASSERT_TRUE(res.result[std::make_pair("e", "e")] == 1);
+  ASSERT_TRUE(res.result[std::make_pair("a", "b")] == 3);
+  ASSERT_TRUE(res.result[std::make_pair("c", "d")] == 3);
+  ASSERT_TRUE(res.result[std::make_pair("c", "e")] == 3);
+  ASSERT_TRUE(res.result[std::make_pair("e", "d")] == 0);
+}
+
+TEST(PowTransTest, transition_matrix) {
+  CXXGraph::Node<int> n1("a", 1); // deg 2
+  CXXGraph::Node<int> n2("b", 1); // get 3
+  CXXGraph::Node<int> n3("c", 1); // deg 2
+  CXXGraph::Node<int> n4("d", 1); // get 0
+
+  CXXGraph::UndirectedEdge<int> e1(0, n1, n2);
+  CXXGraph::UndirectedEdge<int> e2(1, n1, n3);
+  CXXGraph::UndirectedEdge<int> e3(2, n2, n3);
+  CXXGraph::DirectedEdge<int> e4(3, n2, n4);
+  CXXGraph::DirectedEdge<int> e5(4, n4, n1);
+  CXXGraph::DirectedEdge<int> e6(5, n4, n2);
+  CXXGraph::DirectedEdge<int> e7(6, n4, n3);
+
+  CXXGraph::Graph<int> graph;
+  graph.addEdge(&e1);
+  graph.addEdge(&e2);
+  graph.addEdge(&e3);
+  graph.addEdge(&e4);
+  graph.addEdge(&e5);
+  graph.addEdge(&e6);
+  graph.addEdge(&e7);
+
+  auto res = graph.powTransitionMatrix(10);
+
+  const double threshold = 1e-3;
+
+  ASSERT_TRUE(res.success);
+  ASSERT_NEAR(res.result[std::make_pair("a", "a")], 0.286, threshold);
+  ASSERT_NEAR(res.result[std::make_pair("a", "b")], 0.321, threshold);
+  ASSERT_NEAR(res.result[std::make_pair("c", "d")], 0.107, threshold);
+  ASSERT_NEAR(res.result[std::make_pair("c", "c")], 0.286, threshold);
+  ASSERT_NEAR(res.result[std::make_pair("d", "a")], 0.286, threshold);
+  ASSERT_NEAR(res.result[std::make_pair("a", "d")], 0.107, threshold);
+}