From 1ffa9abdd66c5c80c8825e3beda302167e969e79 Mon Sep 17 00:00:00 2001 From: RonakGSahu Date: Sat, 22 Jul 2023 18:21:50 -0400 Subject: [PATCH 1/2] Changed from back tracking algorithm to forward tracking algorithm --- src/paths/steiner.cpp | 302 ++++++------------------- tests/unit/igraph_steiner_tree_fpt.c | 21 +- tests/unit/igraph_steiner_tree_fpt.out | 5 + 3 files changed, 96 insertions(+), 232 deletions(-) diff --git a/src/paths/steiner.cpp b/src/paths/steiner.cpp index cf697cb247..c1ca012f89 100644 --- a/src/paths/steiner.cpp +++ b/src/paths/steiner.cpp @@ -128,194 +128,6 @@ static int_set fetchSetsBasedonIndex(igraph_integer_t index, const dictionary &s IGRAPH_FATAL("The index that you tried to find doesn't exist. Hence the code won't run."); } -/* - * Calculating factorial of a number. - */ -static igraph_integer_t factorial(igraph_integer_t n) { - igraph_integer_t answer = 1; - for (igraph_integer_t i = 1; i <= n; i++) { - answer *= i; - } - return answer; -} - -/* - * Finding number of combinations nCr - * Used to determine number of elements to be scanned in DP table - * for a particular subset during generation of Steiner Tree. - */ -static igraph_integer_t Combination(igraph_integer_t n, igraph_integer_t r) { - return factorial(n) / (factorial(n - r) * factorial(r)); -} - -/** - * Iterates through a particular row of the DP table - * find the value of column where distance (q,k) + distance(Sk(D)) - * where D is a subset. This signifies a vertex that is part - * of minimum path from vertex q to D - * - * \param dp_cache The DP table. - * \param indexD Index of the subset D. - * \param q vertex from who minimum path needs tp be calculated - */ -static igraph_integer_t findMinimumK(igraph_matrix_t *dp_cache, igraph_integer_t indexD, igraph_integer_t q, const dictionary &subsetMap) { - - igraph_integer_t min_col_num = -1; - igraph_real_t min_sum_for_col; - - int_set D = fetchSetsBasedonIndex(indexD, subsetMap); - - for (igraph_integer_t i = 0; i < dp_cache->ncol; i++) { - if (q != i) { - - if (min_col_num == -1) { - min_col_num = i; - min_sum_for_col = MATRIX(*dp_cache, q, i) + MATRIX(*dp_cache, indexD, i); - } else if (MATRIX(*dp_cache, q, i) + MATRIX(*dp_cache, indexD, i) < min_sum_for_col) { - min_col_num = i; - min_sum_for_col = MATRIX(*dp_cache, q, i) + MATRIX(*dp_cache, indexD, i); - } else if ((MATRIX(*dp_cache, q, i) + MATRIX(*dp_cache, indexD, i) == min_sum_for_col) && (D.find(i) != D.end())) { - min_sum_for_col = MATRIX(*dp_cache, q, i) + MATRIX(*dp_cache, indexD, i); - min_col_num = q; - } - } - } - - return min_col_num; -} - -/** - * \function generate_steiner_tree_exact - * - * Generation of the Steiner Tree based on calculation of minimum distances in DP table. - * - * \param graph The graph object. - * \param weights The edge weights. All edge weights must be - * non-negative. Additionally, no edge weight may be NaN. - * \param dp_cache The DP table. - * \param indexD The index of subset D in DP table. - * \param q The vertex that was remmoved from steiner terminals. - * \param vertexlist_all The vector to capture vertices in resultant Steiner Tree. - * \param edgelist_all The vector to capture edges in resultant Steiner Tree. - */ -static igraph_error_t generate_steiner_tree_exact(const igraph_t *graph, const igraph_vector_t *weights, - igraph_matrix_t *dp_cache, igraph_integer_t indexD, igraph_integer_t q, igraph_vector_int_t *edgelist_all, const dictionary &subsetMap) { - - int_set C = fetchSetsBasedonIndex(indexD, subsetMap); - // Initially the value of m is the vertex that was removed from Steiner Terminals - igraph_integer_t m = q; - int_set D = C; - std::set edgelist_all_set; - while (D.size() > 1) { - std::set::iterator itr; - indexD = fetchIndexofMapofSets(D, subsetMap); - - // Finding the bridge vertex from m to subset D which is part of shortest path. - igraph_integer_t k = findMinimumK(dp_cache, indexD, m, subsetMap); - - igraph_vector_int_t edgelist; - IGRAPH_VECTOR_INT_INIT_FINALLY(&edgelist, 0); - - IGRAPH_CHECK(igraph_get_shortest_path_dijkstra(graph, nullptr, &edgelist, m, k, weights, IGRAPH_ALL)); - - for (int i = 0; i < igraph_vector_int_size(&edgelist); i++) { - edgelist_all_set.insert(VECTOR(edgelist)[i]); - } - - igraph_integer_t min_E_value = IGRAPH_INTEGER_MAX; - int_set min_F; - /* - * When the size of subset is > 2 we need to split the subset into E and F - * where E is singleton and F is subset of size (D.size() - 1) - * and repeat same process - */ - if (D.size() > 2) { - - /* - We iterate through the subset D and split it into E and F where E is singleton set during every iteration. - This process leads to find out value where distance (E,k) + (F,k) is minimum. - */ - igraph_real_t min_value = IGRAPH_INFINITY; - igraph_integer_t D_size = D.size(); - for (igraph_integer_t i = 0; i < D_size; i++) { - igraph_integer_t E_raw = *next(D.begin(), i); - int_set F; - int_set E; - E.insert(E_raw); - std::set_difference(D.begin(), D.end(), E.begin(), E.end(), std::inserter(F, F.end())); - igraph_integer_t indexF = fetchIndexofMapofSets(D, subsetMap); - igraph_real_t temp_value = MATRIX(*dp_cache, *E.begin(), k) + MATRIX(*dp_cache, indexF, k); - if (temp_value < min_value) { - min_value = temp_value; - min_E_value = *E.begin(); - min_F = F; - } - } - - igraph_vector_int_t edgelist_1; - IGRAPH_VECTOR_INT_INIT_FINALLY(&edgelist_1, 0); - - IGRAPH_CHECK(igraph_get_shortest_path_dijkstra(graph, nullptr, &edgelist_1, k, min_E_value, weights, IGRAPH_ALL)); - - for (int i = 0; i < igraph_vector_int_size(&edgelist_1); i++) { - edgelist_all_set.insert(VECTOR(edgelist_1)[i]); - } - - igraph_vector_int_destroy(&edgelist_1); - IGRAPH_FINALLY_CLEAN(1); - } else { - /* - If the size of subset is 2 then just shortest path from - k to first element of subset and k to second element of subset is sufficient. - */ - igraph_integer_t E1, F1; - - E1 = *D.begin(); - F1 = *next(D.begin(), 1); - - - igraph_vector_int_t edgelist_1; - IGRAPH_VECTOR_INT_INIT_FINALLY(&edgelist_1, 0); - - IGRAPH_CHECK(igraph_get_shortest_path_dijkstra(graph, nullptr, &edgelist_1, k, E1, weights, IGRAPH_ALL)); - - igraph_vector_int_t edgelist_2; - IGRAPH_VECTOR_INT_INIT_FINALLY(&edgelist_2, 0); - - IGRAPH_CHECK(igraph_get_shortest_path_dijkstra(graph, nullptr, &edgelist_2, k, F1, weights, IGRAPH_ALL)); - - for (int i = 0; i < igraph_vector_int_size(&edgelist_1); i++) { - edgelist_all_set.insert(VECTOR(edgelist_1)[i]); - } - - for (int i = 0; i < igraph_vector_int_size(&edgelist_2); i++) { - edgelist_all_set.insert(VECTOR(edgelist_2)[i]); - } - - - igraph_vector_int_destroy(&edgelist_1); - igraph_vector_int_destroy(&edgelist_2); - - IGRAPH_FINALLY_CLEAN(2); - - min_F.insert(F1); - } - - m = k; - D = min_F; - - igraph_vector_int_destroy(&edgelist); - - IGRAPH_FINALLY_CLEAN(1); - } - - for (auto i : edgelist_all_set) { - igraph_vector_int_push_back(edgelist_all, i); - } - - return IGRAPH_SUCCESS; -} - /** * \function igraph_steiner_dreyfus_wagner @@ -348,7 +160,7 @@ static igraph_error_t generate_steiner_tree_exact(const igraph_t *graph, const i * that are part of Steiner tree. * \return Error code. * - * Time complexity: O( 3^k ∗ V + 2^k ∗ V^2 + V∗(V+E) ∗ log(V) ) = O(3^k) + * Time complexity: O( 3^k ∗ V + 2^k ∗ V^2 + V∗(V+E) ∗ log(V) ) * where V and E are the number of vertices and edges * and k is the number of Steiner terminals. * It is recommended that V <= 50 and k < 11. @@ -574,15 +386,26 @@ igraph_error_t igraph_steiner_dreyfus_wagner( IGRAPH_CHECK(igraph_matrix_init(&dp_cache, no_of_nodes + pow(2, igraph_vector_int_size(&steiner_terminals_copy)) - (igraph_vector_int_size(&steiner_terminals_copy) + 1), no_of_nodes)); IGRAPH_FINALLY(igraph_matrix_destroy, &dp_cache); + std::vector > distance_matrix = std::vector >(no_of_nodes + pow(2, igraph_vector_int_size(&steiner_terminals_copy)) - (igraph_vector_int_size(&steiner_terminals_copy) + 1), std::vector(no_of_nodes, int_set())); + /* Addition to the dp_cache where at the same indices of the matrix we are storing the set of edge ids that are the shortest path to it + */ igraph_matrix_fill(&dp_cache, IGRAPH_INFINITY); /* * for singleton value the distance in dp cahce is just the * distance between same vertices in distance matrix + * and the set of edges are just the edge between the vertices that we get with dijkstra algorithm */ for (igraph_integer_t i = 0; i < no_of_nodes; i++) { IGRAPH_ALLOW_INTERRUPTION(); for (igraph_integer_t j = 0; j < no_of_nodes; j++) { MATRIX(dp_cache, i, j) = MATRIX(distance, i, j); + igraph_vector_int_t edgelist; + IGRAPH_VECTOR_INT_INIT_FINALLY(&edgelist, 0); + igraph_get_shortest_path_dijkstra(graph, nullptr, &edgelist, i, j, pweights, IGRAPH_ALL); + for (int k = 0 ; k < igraph_vector_int_size(&edgelist) ; k++) { + distance_matrix[i][j].insert(igraph_vector_int_get(&edgelist, k)); + } + IGRAPH_FINALLY_CLEAN(1); } } @@ -612,36 +435,47 @@ igraph_error_t igraph_steiner_dreyfus_wagner( MATRIX(dp_cache, indexOfSubsetD, j) = IGRAPH_INFINITY; } + int_set D_prime = D; + D_prime.erase(*D.begin()); + std::set Subsets = std::set(); + generateD_E(D_prime, Subsets, *D.begin()); + // E are subsets of D such that D[1] is in E and E is subset of D where E != D + for (igraph_integer_t j = 0; j < no_of_nodes; j++) { igraph_real_t distance1 = IGRAPH_INFINITY; - - int_set D_prime = D; - D_prime.erase(*D.begin()); - - std::set Subsets = std::set(); - ; - generateD_E(D_prime, Subsets, *D.begin()); - // E are subsets of D such that D[1] is in E + int_set set_u; for (auto E : Subsets) { - igraph_integer_t indexOfSubsetE = (E.size() == 1) ? *E.begin() : fetchIndexofMapofSets(E, subsetMap); igraph_real_t distanceEJ = MATRIX(dp_cache, indexOfSubsetE, j); - int_set DMinusE = D; + int_set DMinusE; - for (auto elem : E) { - // can be modified with a set difference function as we did in tree generation - DMinusE.erase(elem); - } + std::set_difference(D.begin(), D.end(), E.begin(), E.end(), std::inserter(DMinusE, DMinusE.end())); igraph_integer_t indexOfSubsetDMinusE = (DMinusE.size() == 1) ? *DMinusE.begin() : fetchIndexofMapofSets(DMinusE, subsetMap); - distance1 = std::min(distance1, distanceEJ + MATRIX(dp_cache, indexOfSubsetDMinusE, j)); + if ((distanceEJ + MATRIX(dp_cache, indexOfSubsetDMinusE, j)) < distance1) { + distance1 = distanceEJ + MATRIX(dp_cache, indexOfSubsetDMinusE, j); + + //get the reference to this set and combine them + auto& setEJ = distance_matrix[indexOfSubsetE][j]; + auto& setDMinusEJ = distance_matrix[indexOfSubsetDMinusE][j]; + set_u.clear(); + std::set_union(setEJ.begin(), setEJ.end(), setDMinusEJ.begin(), setDMinusEJ.end(), std::inserter(set_u, set_u.end())); + } } for (igraph_integer_t k = 0; k < no_of_nodes; k++) { - MATRIX(dp_cache, indexOfSubsetD, k) = std::min(MATRIX(dp_cache, indexOfSubsetD, k), MATRIX(distance, k, j) + distance1); + if ((MATRIX(distance, k, j) + distance1) < MATRIX(dp_cache, indexOfSubsetD, k)) { + MATRIX(dp_cache, indexOfSubsetD, k) = MATRIX(distance, k, j) + distance1; + + //get the reference to this set and combine them + auto& setKJ = distance_matrix[k][j]; + distance_matrix[indexOfSubsetD][k].clear(); + std::set_union(setKJ.begin(), setKJ.end(), set_u.begin(), set_u.end(), std::inserter(distance_matrix[indexOfSubsetD][k], distance_matrix[indexOfSubsetD][k].end())); + + } } } } @@ -651,48 +485,59 @@ igraph_error_t igraph_steiner_dreyfus_wagner( std::set E_subsets = std::set(); int_set C_prime; int_set C; - igraph_integer_t C_1 = igraph_vector_int_get(&steiner_terminals_copy,0); - for (int i = 0 ; i < igraph_vector_int_size(&steiner_terminals_copy) ; i++){ - C.insert(igraph_vector_int_get(&steiner_terminals_copy,i)); - if (i!=0) - C_prime.insert(igraph_vector_int_get(&steiner_terminals_copy,i)); + igraph_integer_t C_1 = igraph_vector_int_get(&steiner_terminals_copy, 0); + for (int i = 0 ; i < igraph_vector_int_size(&steiner_terminals_copy) ; i++) { + C.insert(igraph_vector_int_get(&steiner_terminals_copy, i)); + if (i != 0) { + C_prime.insert(igraph_vector_int_get(&steiner_terminals_copy, i)); + } } - generateD_E(C_prime, E_subsets,C_1); + generateD_E(C_prime, E_subsets, C_1); // E are subsets of C such that C[1] is in E and E is subset of C where E != C - igraph_real_t distance2 = IGRAPH_INFINITY; + igraph_real_t distance2 = IGRAPH_INFINITY; + int_set final_set; for (igraph_integer_t j = 0; j < no_of_nodes; j++) { igraph_real_t distance1 = IGRAPH_INFINITY; IGRAPH_ALLOW_INTERRUPTION(); + + int_set set_u; + for (auto E : E_subsets) { - igraph_integer_t indexE = (E.size() == 1) ? *E.begin() : fetchIndexofMapofSets(E,subsetMap); + igraph_integer_t indexE = (E.size() == 1) ? *E.begin() : fetchIndexofMapofSets(E, subsetMap); + int_set CMinusE; - std::set_difference(C.begin(),C.end(),E.begin(),E.end(),std::inserter(CMinusE,CMinusE.end())); - igraph_integer_t indexC_E = (CMinusE.size() == 1) ? *CMinusE.begin() : fetchIndexofMapofSets(CMinusE,subsetMap); - igraph_real_t distanceFJ = MATRIX(dp_cache,indexE,j); + std::set_difference(C.begin(), C.end(), E.begin(), E.end(), std::inserter(CMinusE, CMinusE.end())); + igraph_integer_t indexC_E = (CMinusE.size() == 1) ? *CMinusE.begin() : fetchIndexofMapofSets(CMinusE, subsetMap); + igraph_real_t distanceEJ = MATRIX(dp_cache, indexE, j); + if (((distanceEJ + (MATRIX(dp_cache, indexC_E, j))) < distance1)) { - if (((distanceFJ + (MATRIX(dp_cache, indexC_E, j))) < distance1)) { - distance1 = distanceFJ + (MATRIX(dp_cache, indexC_E, j)); + distance1 = distanceEJ + (MATRIX(dp_cache, indexC_E, j)); + + //get the reference to this set and combine them + set_u.clear(); + auto& setEJ = distance_matrix[indexE][j]; + auto& setCMinusEJ = distance_matrix[indexC_E][j]; + std::set_union(setEJ.begin(), setEJ.end(), setCMinusEJ.begin(), setCMinusEJ.end(), std::inserter(set_u, set_u.end())); } } if ((MATRIX(distance, q, j) + distance1) < distance2) { distance2 = MATRIX(distance, q, j) + distance1; + //get the reference to this set and combine them + final_set.clear(); + auto& setQJ = distance_matrix[q][j]; + std::set_union(setQJ.begin(), setQJ.end(), set_u.begin(), set_u.end(), std::inserter(final_set, final_set.end())); } } *res = distance2; - - int_set newSet; - - for (igraph_integer_t i = 0; i < steiner_terminals_copy_size; i++) { - newSet.insert(VECTOR(steiner_terminals_copy)[i]); - } - igraph_integer_t indexD = fetchIndexofMapofSets(newSet, subsetMap); - + //put the edges in the res_tree and move the final value to res which was stored in distance2 for the scope of this function igraph_vector_int_clear(res_tree); - IGRAPH_CHECK(generate_steiner_tree_exact(graph, pweights, &dp_cache, indexD, q, res_tree, subsetMap)); + for (auto elem : final_set) { + igraph_vector_int_push_back(res_tree, elem); + } igraph_matrix_destroy(&distance); igraph_vector_int_destroy(&steiner_terminals_copy); igraph_matrix_destroy(&dp_cache); @@ -704,7 +549,6 @@ igraph_error_t igraph_steiner_dreyfus_wagner( IGRAPH_FINALLY_CLEAN(1); return IGRAPH_SUCCESS; } - IGRAPH_HANDLE_EXCEPTIONS_END; return IGRAPH_SUCCESS; -} +} \ No newline at end of file diff --git a/tests/unit/igraph_steiner_tree_fpt.c b/tests/unit/igraph_steiner_tree_fpt.c index 6cd0202de2..58e4463b81 100644 --- a/tests/unit/igraph_steiner_tree_fpt.c +++ b/tests/unit/igraph_steiner_tree_fpt.c @@ -385,16 +385,31 @@ int main(void) { igraph_vector_int_print(&tree_edges_new); printf("value: %f\n", value_new_1); + printf("\nA simple square graph with few more edges outside the square\n"); + igraph_t g_new_2; + igraph_real_t value_new_2; + igraph_small(&g_new_2, 6, IGRAPH_UNDIRECTED, + 0, 1, + 1, 2, 1, 3, + 3, 4, + 2, 4, + 4, 5, + -1 + ); + igraph_vector_int_init_int(&terminals_new, 4, 0, 3, 4, 5); + igraph_vector_int_init(&tree_edges_new, 0); + igraph_steiner_dreyfus_wagner(&g_new_2, &terminals_new, NULL, &value_new_2, &tree_edges_new); + printf("Tree edges:\n"); + igraph_vector_int_print(&tree_edges_new); + printf("value: %f\n", value_new_2); igraph_destroy(&g_new); igraph_destroy(&g_new_1); + igraph_destroy(&g_new_2); igraph_vector_int_destroy(&tree_edges_new); igraph_vector_int_destroy(&terminals_new); - - - VERIFY_FINALLY_STACK(); return 0; diff --git a/tests/unit/igraph_steiner_tree_fpt.out b/tests/unit/igraph_steiner_tree_fpt.out index 2e4bd0add2..fd0ed2a4b9 100644 --- a/tests/unit/igraph_steiner_tree_fpt.out +++ b/tests/unit/igraph_steiner_tree_fpt.out @@ -69,3 +69,8 @@ A graph with different structure that previous Tree edges: 0 1 2 4 6 7 value: 6.000000 + +A simple square graph with few more edges outside the square +Tree edges: +0 2 3 5 +value: 4.000000 \ No newline at end of file From e7ca410e3ba7ada82e23bf1a486e36951b38236d Mon Sep 17 00:00:00 2001 From: RonakGSahu Date: Sat, 22 Jul 2023 18:35:54 -0400 Subject: [PATCH 2/2] Added the testcase that was on PR which was previously failing --- tests/unit/igraph_steiner_tree_fpt.c | 28 ++++++++++++++++++++++++++ tests/unit/igraph_steiner_tree_fpt.out | 7 ++++++- 2 files changed, 34 insertions(+), 1 deletion(-) diff --git a/tests/unit/igraph_steiner_tree_fpt.c b/tests/unit/igraph_steiner_tree_fpt.c index 58e4463b81..7f0e1b460b 100644 --- a/tests/unit/igraph_steiner_tree_fpt.c +++ b/tests/unit/igraph_steiner_tree_fpt.c @@ -385,6 +385,9 @@ int main(void) { igraph_vector_int_print(&tree_edges_new); printf("value: %f\n", value_new_1); + igraph_vector_int_destroy(&tree_edges_new); + igraph_vector_int_destroy(&terminals_new); + printf("\nA simple square graph with few more edges outside the square\n"); igraph_t g_new_2; igraph_real_t value_new_2; @@ -404,11 +407,36 @@ int main(void) { printf("value: %f\n", value_new_2); + igraph_vector_int_destroy(&tree_edges_new); + igraph_vector_int_destroy(&terminals_new); + + printf("\nA 3 vertices square graph with weights\n"); + igraph_t g_new_3; + igraph_vector_t weights; + igraph_vector_init(&weights, 12); + igraph_vector_fill(&weights, 100); + igraph_small(&g_new_3, 9, IGRAPH_UNDIRECTED, + 0, 1, 1, 2, + 0, 3, 1, 4, 2, 5, + 3, 4, 4, 5, + 3, 6, 4, 7, 5, 8, + 6, 7, 7, 8, + -1); + igraph_vector_int_init_int_end(&terminals_new, -1, 0, 4, 2, 8, -1); + igraph_vector_int_init(&tree_edges_new, 0); + igraph_steiner_dreyfus_wagner(&g_new_3, &terminals_new, &weights, &value_new, &tree_edges_new); + printf("Tree edges:\n"); + igraph_vector_int_print(&tree_edges_new); + printf("value: %f\n", value_new); + + igraph_destroy(&g_new); igraph_destroy(&g_new_1); igraph_destroy(&g_new_2); + igraph_destroy(&g_new_3); igraph_vector_int_destroy(&tree_edges_new); igraph_vector_int_destroy(&terminals_new); + igraph_vector_destroy(&weights); VERIFY_FINALLY_STACK(); diff --git a/tests/unit/igraph_steiner_tree_fpt.out b/tests/unit/igraph_steiner_tree_fpt.out index fd0ed2a4b9..8385bc2222 100644 --- a/tests/unit/igraph_steiner_tree_fpt.out +++ b/tests/unit/igraph_steiner_tree_fpt.out @@ -73,4 +73,9 @@ value: 6.000000 A simple square graph with few more edges outside the square Tree edges: 0 2 3 5 -value: 4.000000 \ No newline at end of file +value: 4.000000 + +A 3 vertices square graph with weights +Tree edges: +0 1 3 8 11 +value: 500.000000