diff --git a/include/CXXGraph/Graph/Algorithm/Pow_impl.hpp b/include/CXXGraph/Graph/Algorithm/Pow_impl.hpp new file mode 100644 index 00000000..d299f6f0 --- /dev/null +++ b/include/CXXGraph/Graph/Algorithm/Pow_impl.hpp @@ -0,0 +1,221 @@ +/***********************************************************/ +/*** ______ ____ ______ _ ***/ +/*** / ___\ \/ /\ \/ / ___|_ __ __ _ _ __ | |__ ***/ +/*** | | \ / \ / | _| '__/ _` | '_ \| '_ \ ***/ +/*** | |___ / \ / \ |_| | | | (_| | |_) | | | | ***/ +/*** \____/_/\_\/_/\_\____|_| \__,_| .__/|_| |_| ***/ +/*** |_| ***/ +/***********************************************************/ +/*** Header-Only C++ Library for Graph ***/ +/*** Representation and Algorithms ***/ +/***********************************************************/ +/*** Author: ZigRazor ***/ +/*** E-Mail: zigrazor@gmail.com ***/ +/***********************************************************/ +/*** Collaboration: ----------- ***/ +/***********************************************************/ +/*** License: AGPL v3.0 ***/ +/***********************************************************/ + +#ifndef __CXXGRAPH_POW_IMPL_H__ +#define __CXXGRAPH_POW_IMPL_H__ + +#pragma once + +#include "CXXGraph/Graph/Graph_decl.h" + +/** + * @brief simple matrix multiplication of two matrices + * A and B + * @param a matrix A + * @param b matrix B + * @return A times B + */ +template +std::vector> matMult(const std::vector> &a, + const std::vector> &b) { + static_assert(std::is_arithmetic::value, + "Type T must be an arithmetic type"); + + // two square matrices both of size N x N where N > 0 + if (a.empty() || a[0].size() != b.size() || a.size() != a[0].size() || + b.size() != b[0].size()) { + throw std::invalid_argument( + "Matrix must have valid dimensions and be at least 1x1."); + } + + int n = static_cast(a.size()); // N x N matrix + std::vector> res(n, std::vector(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; +} + +/** + * @brief exponentation takes a matrix of arithmetic type and + * raises it to the power of k. + * @param mat a square matrix + * @param k a nonnegative integer + * @return M, the result of mat^k + */ +template +std::vector> exponentiation(std::vector> mat, + unsigned int k) { + static_assert(std::is_arithmetic::value, + "Type T must be any arithmetic type"); + + // validate size and shape of matrix + if (mat.size() == 0 || mat.size() != mat[0].size()) { + throw std::invalid_argument("Matrix must be square and at least 1x1."); + } + + int n = static_cast(mat.size()); + std::vector> res(n, std::vector(n, 0)); + + // build identity matrix + for (int i = 0; i < n; i++) res[i][i] = 1; + + // fast exponentiation! + while (k) { + if (k & 1) res = matMult(res, mat); + mat = matMult(mat, mat); + k >>= 1; + } + + return res; +} + +namespace CXXGraph { + +/** + * @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 + */ +template +const PowAdjResult matrixPow(const shared> &adj, + unsigned int k) { + PowAdjResult result; + result.success = false; + result.errorMessage = ""; + + // convert back and forth between user ids and index in temporary adj matrix + std::unordered_map userIdToIdx; + std::unordered_map idxToUserId; + + int n = 0; + for (const auto &[node, _] : *adj) { + userIdToIdx[node->getUserId()] = n; + idxToUserId[n] = node->getUserId(); + n++; + } + + // adj int will store the temporary (integer) adjacency matrix + std::vector> tempIntAdj( + n, std::vector(n, 0)); + + // populate temporary adjacency matrix w/ edges + // can handle both directed and undirected + for (const auto &[_, edges] : *adj) { + for (const auto &[_, edge] : edges) { + const auto &[u, v] = edge->getNodePair(); + const int uIdx = userIdToIdx[u->getUserId()]; + const int vIdx = userIdToIdx[v->getUserId()]; + + // if undirected, add both sides + if (!(edge->isDirected().has_value() && edge->isDirected().value())) + tempIntAdj[vIdx][uIdx] = 1; + tempIntAdj[uIdx][vIdx] = 1; + } + } + + // calculate the power matrix + const auto powerMatrix = exponentiation(tempIntAdj, k); + + // remap values + std::unordered_map, unsigned long long, + CXXGraph::pair_hash> + powAdj; + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + powAdj[std::make_pair(idxToUserId[i], idxToUserId[j])] = + powerMatrix[i][j]; + } + } + + result.result = std::move(powAdj); + result.success = true; + + return result; +} + +/** + * @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 + */ +template +const PowTransResult matrixPow(const shared> &trans, + unsigned int k) { + PowTransResult result; + result.success = false; + result.errorMessage = ""; + + // get a map between index in adj matrix + // and userId + std::unordered_map userIdToIdx; + std::unordered_map idxToUserId; + + int n = 0; + for (const auto &[node, _] : *trans) { + userIdToIdx[node->getUserId()] = n; + idxToUserId[n] = node->getUserId(); + n++; + } + + std::vector> stochasticMatrix(n, + std::vector(n, 0.0)); + + // given transition matrix, convert it to + // stochastic matrix + for (const auto &[u, edges] : *trans) { + const int uIdx = userIdToIdx[u->getUserId()]; + + for (const auto &[v, weight] : edges) { + const int vIdx = userIdToIdx[v->getUserId()]; + stochasticMatrix[uIdx][vIdx] = weight; + } + } + + // exponentiate stochastic matrix + auto powerTransMatrix = exponentiation(stochasticMatrix, k); + + // turn back into a map between nodes + std::unordered_map, double, + CXXGraph::pair_hash> + powTrans; + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + powTrans[std::make_pair(idxToUserId[i], idxToUserId[j])] = + powerTransMatrix[i][j]; + } + } + + result.result = std::move(powTrans); + result.success = true; + + return result; +} +} // namespace CXXGraph + +#endif //_CXXGRAPH_POW_IMPL_H__ \ No newline at end of file diff --git a/include/CXXGraph/Graph/Graph.h b/include/CXXGraph/Graph/Graph.h index aedb55ee..2ed02e8d 100644 --- a/include/CXXGraph/Graph/Graph.h +++ b/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" diff --git a/include/CXXGraph/Utility/Typedef.hpp b/include/CXXGraph/Utility/Typedef.hpp index b26440b5..cd1297f2 100644 --- a/include/CXXGraph/Utility/Typedef.hpp +++ b/include/CXXGraph/Utility/Typedef.hpp @@ -237,6 +237,26 @@ struct BestFirstSearchResult_struct { template using BestFirstSearchResult = BestFirstSearchResult_struct; +/// Struct that contains the results from an adjacency matrix exponentiation +/// results +struct PowAdjResult_struct { + bool success = false; + std::string errorMessage = ""; + std::unordered_map, 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, double, pair_hash> + result = {}; +}; +typedef PowTransResult_struct PowTransResult; + /////////////////////////////////////////////////////////////////////////////////// // Using Definition // /////////////////////////////////////////////////////////////// diff --git a/test/PowTest.cpp b/test/PowTest.cpp new file mode 100644 index 00000000..20bf3287 --- /dev/null +++ b/test/PowTest.cpp @@ -0,0 +1,162 @@ + +#include + +#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 +using unique = std::unique_ptr; +template +using shared = std::shared_ptr; + +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 a("a", 1); + CXXGraph::Node b("b", 1); + CXXGraph::Node c("c", 1); + CXXGraph::Node d("d", 1); + + CXXGraph::DirectedEdge e1(0, a, b); + CXXGraph::DirectedEdge e2(1, a, d); + CXXGraph::DirectedEdge e3(2, b, c); + CXXGraph::DirectedEdge e4(3, b, d); + CXXGraph::DirectedEdge e5(4, c, b); + CXXGraph::DirectedEdge e6(5, d, c); + + CXXGraph::Graph graph; + graph.addEdge(&e1); + graph.addEdge(&e2); + graph.addEdge(&e3); + graph.addEdge(&e4); + graph.addEdge(&e5); + graph.addEdge(&e6); + + CXXGraph::PowAdjResult res = matrixPow(graph.getAdjMatrix(), 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); +} + +TEST(PowAdjTest, triangle) { + CXXGraph::Node a("a", 1); + CXXGraph::Node b("b", 1); + CXXGraph::Node c("c", 1); + + CXXGraph::UndirectedEdge e1(0, a, b); + CXXGraph::UndirectedEdge e2(1, b, c); + CXXGraph::UndirectedEdge e3(2, c, a); + + CXXGraph::Graph graph; + graph.addEdge(&e1); + graph.addEdge(&e2); + graph.addEdge(&e3); + + CXXGraph::PowAdjResult res = matrixPow(graph.getAdjMatrix(), 5); + + ASSERT_TRUE(res.success); + ASSERT_TRUE(res.result[std::make_pair("a", "a")] == 10); + ASSERT_TRUE(res.result[std::make_pair("b", "b")] == 10); + ASSERT_TRUE(res.result[std::make_pair("c", "c")] == 10); + ASSERT_TRUE(res.result[std::make_pair("a", "b")] == 11); + ASSERT_TRUE(res.result[std::make_pair("b", "a")] == 11); + ASSERT_TRUE(res.result[std::make_pair("b", "c")] == 11); + ASSERT_TRUE(res.result[std::make_pair("c", "b")] == 11); + ASSERT_TRUE(res.result[std::make_pair("a", "c")] == 11); + ASSERT_TRUE(res.result[std::make_pair("c", "a")] == 11); +} + +// https://docs.dgl.ai/generated/dgl.khop_adj.html#dgl.khop_adj +TEST(PowAdjTest, dgl) { + CXXGraph::Node a("a", 1); + CXXGraph::Node b("b", 1); + CXXGraph::Node c("c", 1); + CXXGraph::Node d("d", 1); + CXXGraph::Node e("e", 1); + + CXXGraph::DirectedEdge e1(0, a, a); + CXXGraph::DirectedEdge e2(1, b, b); + CXXGraph::DirectedEdge e3(2, c, c); + CXXGraph::DirectedEdge e4(3, d, d); + CXXGraph::DirectedEdge e5(4, e, e); + CXXGraph::DirectedEdge e6(5, a, b); + CXXGraph::DirectedEdge e7(6, b, c); + CXXGraph::DirectedEdge e8(7, c, d); + CXXGraph::DirectedEdge e9(8, d, e); + CXXGraph::DirectedEdge e10(9, d, a); + + CXXGraph::Graph 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); + + CXXGraph::PowAdjResult res = matrixPow(graph.getAdjMatrix(), 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 n1("a", 1); // deg 2 + CXXGraph::Node n2("b", 1); // get 3 + CXXGraph::Node n3("c", 1); // deg 2 + CXXGraph::Node n4("d", 1); // get 0 + + CXXGraph::UndirectedEdge e1(0, n1, n2); + CXXGraph::UndirectedEdge e2(1, n1, n3); + CXXGraph::UndirectedEdge e3(2, n2, n3); + CXXGraph::DirectedEdge e4(3, n2, n4); + CXXGraph::DirectedEdge e5(4, n4, n1); + CXXGraph::DirectedEdge e6(5, n4, n2); + CXXGraph::DirectedEdge e7(6, n4, n3); + + CXXGraph::Graph graph; + graph.addEdge(&e1); + graph.addEdge(&e2); + graph.addEdge(&e3); + graph.addEdge(&e4); + graph.addEdge(&e5); + graph.addEdge(&e6); + graph.addEdge(&e7); + + CXXGraph::PowTransResult res = matrixPow(graph.getTransitionMatrix(), 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); +}