Merge pull request #390 from dperry17/Min-Cut-dev

Min cut dev

Author: DrTimothyAldenDavis

data/test_zero_cap.mtx

@@ -0,0 +1,13 @@
+%%MatrixMarket matrix coordinate real general
+% Example weighted directed graph adjacency (1-indexed)
+% All weights are positive; some edges explicitly have zero weight
+5 5 9
+1 2 0.500000
+1 4 0.000000
+2 3 1.250000
+2 5 2.000000
+3 1 0.000000
+3 4 3.141590
+4 2 0.750000
+5 4 0.125000
+5 5 1.000000

experimental/algorithm/LAGr_MaxFlow.c

@@ -612,6 +612,7 @@ int LAGr_MaxFlow
     // output:
     double *f,              // max flow from src node to sink node
     GrB_Matrix *flow_mtx,   // optional output flow matrix
+    GrB_Matrix *res_mtx,    // optional output for the residual matrix. Used in the min-cut
     // input:
     LAGraph_Graph G,        // graph to compute maxflow on
     GrB_Index src,          // source node
@@ -711,6 +712,10 @@ int LAGr_MaxFlow
   {
     (*flow_mtx) = NULL ;
   }
+  if (res_mtx != NULL)
+  {
+    (*res_mtx) = NULL ;
+  }
   LG_TRY(LAGraph_CheckGraph(G, msg));
   LG_ASSERT (f != NULL, GrB_NULL_POINTER) ;
   (*f) = 0;
@@ -720,6 +725,7 @@ int LAGr_MaxFlow
   LG_ASSERT_MSG(nrows == n, GrB_INVALID_VALUE, "Matrix must be square");
   LG_ASSERT_MSG(src < n && src >= 0 && sink < n && sink >= 0,
         GrB_INVALID_VALUE, "src and sink must be a value between [0, n)");
+  LG_ASSERT(G->emin_state != LAGraph_BOOLEAN_UNKNOWN, GrB_UNINITIALIZED_OBJECT) ;
   LG_ASSERT_MSG(G->emin > 0, GrB_INVALID_VALUE,
         "the edge weights (capacities) must be greater than 0");
 
@@ -1160,6 +1166,13 @@ int LAGr_MaxFlow
     GRB_TRY(GrB_select(*flow_mtx, NULL, NULL, GrB_VALUEGT_FP64, *flow_mtx, 0, NULL));
   }
 
+  if(res_mtx != NULL){
+   GRB_TRY(GrB_Matrix_new(res_mtx, GrB_FP64, n, n)); 
+   GRB_TRY(GrB_apply(*res_mtx, NULL, NULL, GetResidual, R, NULL)) ;
+    // prune zeros and negative entries from R_hat
+   GRB_TRY(GrB_select(*res_mtx, NULL, NULL, GrB_VALUEGT_FP64, *res_mtx, 0, NULL)) ; 
+  }
+
   //----------------------------------------------------------------------------
   // for test coverage only
   //----------------------------------------------------------------------------

experimental/algorithm/LAGraph_MinCut.c

@@ -0,0 +1,69 @@
+#include <LAGraph.h>
+#include "LG_internal.h"
+#include <LAGraph.h>
+
+#undef LG_FREE_ALL
+#undef LG_FREE_WORK
+
+#define LG_FREE_WORK				\
+{						\
+ GrB_free(&S_diag);				\
+ GrB_free(&S_bar_diag);				\
+ LAGraph_Delete(&G, msg);			\
+}
+
+#define LG_FREE_ALL                             \
+  { LG_FREE_WORK }
+
+int LAGraph_MinCut
+(
+    // outputs
+    GrB_Vector* S,
+    GrB_Vector* S_bar,
+    GrB_Matrix* cut_set,
+    // inputs
+    LAGraph_Graph G_origin,
+    GrB_Matrix R,
+    GrB_Index s,
+    GrB_Index t,
+    char *msg
+)
+{
+  //do a bfs from the source to the sink, stop if the frontier is empty 
+
+  LAGraph_Graph G = NULL;
+  GrB_Matrix S_diag=NULL, S_bar_diag=NULL;
+  GrB_Index n = 0;
+  GrB_Matrix_nrows(&n, R);
+
+  LG_TRY(LAGraph_CheckGraph(G_origin, msg));
+  LG_ASSERT (S != NULL, GrB_NULL_POINTER) ;
+  LG_ASSERT (S_bar != NULL, GrB_NULL_POINTER) ;
+  LG_ASSERT (cut_set != NULL, GrB_NULL_POINTER) ;
+  LG_ASSERT (s >= 0 && t >= 0, GrB_INVALID_VALUE) ;
+
+  //  printf("hit\n");
+  GrB_Matrix A = G_origin->A;
+  
+  LG_TRY(GrB_Matrix_new(cut_set, GrB_FP64, n, n));
+  LG_TRY(GrB_Vector_new(S_bar, GrB_INT64, n));
+  //S is allocated during the bfs
+  
+  LG_TRY(LAGraph_New(&G, &R, LAGraph_ADJACENCY_DIRECTED, msg));
+  LG_TRY(LAGr_BreadthFirstSearch(S, NULL, G, s, msg));
+
+
+  LG_TRY(GrB_assign(*S_bar, *S, NULL, 1, GrB_ALL, n, GrB_DESC_SC));
+  LG_TRY(GrB_assign(*S, *S, NULL, 1, GrB_ALL, n, GrB_DESC_S));
+
+  GRB_TRY(GrB_Matrix_diag(&S_diag, *S, 0));
+  GRB_TRY(GrB_Matrix_diag(&S_bar_diag, *S_bar, 0));
+
+  
+  GRB_TRY(GrB_mxm(*cut_set, NULL, NULL, GxB_PLUS_TIMES_INT64, A, S_bar_diag, NULL));
+  GRB_TRY(GrB_mxm(*cut_set, NULL, NULL, GxB_PLUS_TIMES_INT64, S_diag, *cut_set, NULL));
+
+  
+  LG_FREE_ALL;
+  return (GrB_SUCCESS);
+}

experimental/benchmark/maxflow_demo.c

@@ -58,15 +58,15 @@ int main (int argc, char ** argv){
 
   // LG_SET_BURBLE(1);
   double time = LAGraph_WallClockTime();
-  LAGRAPH_TRY(LAGr_MaxFlow(&flow, NULL, G, S, T, msg));
+  LAGRAPH_TRY(LAGr_MaxFlow(&flow, NULL, NULL, G, S, T, msg));
   time = LAGraph_WallClockTime() - time;
   printf("Time for LAGraph_MaxFlow: %g sec\n", time);
   printf("Max Flow is: %lf\n", flow);
 
   printf("Starting max flow from %" PRIu64 " to %" PRIu64
     ", with flow_matrix returned\n", S, T);
   time = LAGraph_WallClockTime();
-  LAGRAPH_TRY(LAGr_MaxFlow(&flow, &flow_matrix, G, S, T, msg));
+  LAGRAPH_TRY(LAGr_MaxFlow(&flow, &flow_matrix, NULL, G, S, T, msg));
   time = LAGraph_WallClockTime() - time;
   printf("Time for LAGraph_MaxFlow with flow matrix: %g sec\n", time);
   printf("Max Flow is: %lf\n", flow);

experimental/test/test_MaxFlow.c

@@ -30,7 +30,7 @@ GrB_Matrix A = NULL;
 #ifdef GRAPHBLAS_HAS_CUDA
 #define NTESTS 4
 #else
-#define NTESTS 7
+#define NTESTS 8
 #endif
 char filename[LEN + 1];
 
@@ -47,6 +47,7 @@ test_info tests[] = {
   {"matrix_random_flow.mtx", 0,9, 22, LAGraph_ADJACENCY_DIRECTED},
   {"rand.mtx", 0, 19, 37, LAGraph_ADJACENCY_DIRECTED},
   {"mcl.mtx", 0, 9, 0, LAGraph_ADJACENCY_DIRECTED},
+  {"test_zero_cap.mtx", 0, 4, 0.5, LAGraph_ADJACENCY_DIRECTED},
 #ifndef GRAPHBLAS_HAS_CUDA
 // FIXME: the CUDA cases are currently very slow for these matrices,
 // when the GPU is hacked to always be used regardless of problem size:
@@ -84,14 +85,14 @@ void test_MaxFlow(void) {
     // test with JIT
     OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, GxB_JIT_ON));
     double flow = 0;
-    OK(LAGr_MaxFlow(&flow, NULL, G, tests[test].S, tests[test].T, msg));
+    OK(LAGr_MaxFlow(&flow, NULL, NULL, G, tests[test].S, tests[test].T, msg));
     printf("%s\n", msg);
     printf("flow is: %lf\n", flow);
     TEST_CHECK(flow == tests[test].F);
 
     // test without JIT
     OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, GxB_JIT_OFF));
-    OK(LAGr_MaxFlow(&flow, NULL, G, tests[test].S, tests[test].T, msg));
+    OK(LAGr_MaxFlow(&flow, NULL, NULL, G, tests[test].S, tests[test].T, msg));
     TEST_CHECK(flow == tests[test].F);
     OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, GxB_JIT_ON));
 
@@ -107,7 +108,7 @@ void test_MaxFlow(void) {
             {
               printf("src: %d, dest: %d\n", (int) src, (int) dest);
                 if (src == dest) continue ;
-                OK(LAGr_MaxFlow(&flow, NULL, G, src, dest, msg));
+                OK(LAGr_MaxFlow(&flow, NULL, NULL, G, src, dest, msg));
             }
         }
     }
@@ -152,7 +153,7 @@ void test_MaxFlowMtx(void) {
     // test with JIT
     OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, GxB_JIT_ON));
     double flow = 0;
-    OK(LAGr_MaxFlow(&flow, &flow_mtx, G, tests[test].S, tests[test].T, msg));
+    OK(LAGr_MaxFlow(&flow, &flow_mtx, NULL, G, tests[test].S, tests[test].T, msg));
     TEST_CHECK (flow_mtx != NULL) ;
     GxB_print (flow_mtx, 2) ;
     int status = LG_check_flow(flow_mtx, msg);
@@ -165,7 +166,7 @@ void test_MaxFlowMtx(void) {
 
     // test without JIT
     OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, GxB_JIT_OFF));
-    OK(LAGr_MaxFlow(&flow, &flow_mtx, G, tests[test].S, tests[test].T, msg));
+    OK(LAGr_MaxFlow(&flow, &flow_mtx, NULL, G, tests[test].S, tests[test].T, msg));
     TEST_CHECK (flow_mtx != NULL) ;
     status = LG_check_flow(flow_mtx, msg);
     TEST_CHECK (status == GrB_SUCCESS) ;

experimental/test/test_MinCut.c

@@ -0,0 +1,93 @@
+#include <stdio.h>
+#include <acutest.h>
+#include <LAGraphX.h>
+#include <LAGraph_test.h>
+#include <LG_Xtest.h>
+
+
+char msg[LAGRAPH_MSG_LEN];
+LAGraph_Graph G = NULL;
+GrB_Matrix A = NULL;
+#define LEN 512
+#define NTESTS 7
+char filename[LEN + 1];
+
+
+typedef struct {
+  char* filename;
+  GrB_Index s;
+  GrB_Index t;
+  LAGraph_Kind kind;
+} test_case;
+
+test_case tests[] = {
+  {"wiki.mtx", 0, 5, LAGraph_ADJACENCY_DIRECTED},
+  {"matrix_random_flow.mtx", 0,9, LAGraph_ADJACENCY_DIRECTED},
+  {"rand.mtx", 0, 19, LAGraph_ADJACENCY_DIRECTED},
+  {"mcl.mtx", 0, 9, LAGraph_ADJACENCY_DIRECTED},
+  {"cycle_flow.mtx", 0, 89,  LAGraph_ADJACENCY_DIRECTED},
+  {"random_weighted_general2.mtx", 0, 299, LAGraph_ADJACENCY_UNDIRECTED},
+  {"random_weighted_general1.mtx", 0, 499, LAGraph_ADJACENCY_UNDIRECTED}
+};
+
+
+void test_MinCut() {
+
+  LAGraph_Init(msg);
+#if LG_SUITESPARSE_GRAPHBLAS_V10
+  OK(LG_SET_BURBLE(0));
+  for(uint8_t test = 0; test < NTESTS; test++){
+    GrB_Matrix A=NULL, R=NULL, cut_set=NULL;
+    GrB_Vector S=NULL, S_bar=NULL;
+    double min_cut = 0;
+    GrB_Index n = 0;
+    printf ("\nMatrix: %s\n", tests[test].filename);
+    TEST_CASE(tests[test].filename);
+    
+    snprintf(filename, LEN, LG_DATA_DIR "%s", tests[test].filename);
+    FILE* f = fopen(filename, "r");
+    TEST_CHECK(f != NULL);
+    
+    OK(LAGraph_MMRead(&A, f, msg));
+    
+    OK(fclose(f));
+    LAGraph_Kind kind = tests [test].kind ;
+    OK(LAGraph_New(&G, &A, kind, msg));
+    if (kind == LAGraph_ADJACENCY_DIRECTED)
+    {
+        OK(LAGraph_Cached_AT(G, msg));
+    }
+
+    OK(LAGraph_Cached_EMin(G, msg));
+
+    // test with JIT
+    OK(GxB_Global_Option_set(GxB_JIT_C_CONTROL, GxB_JIT_ON));
+    double flow = 0;
+    OK(LAGr_MaxFlow(&flow, NULL, &R, G, tests[test].s, tests[test].t, msg));
+    printf("%s\n", msg);
+    printf("flow is: %lf\n", flow);
+
+    OK(LAGraph_MinCut(&S, &S_bar, &cut_set, G, R, tests[test].s, tests[test].t, msg));
+    printf("%s\n", msg);
+
+    OK(GrB_reduce(&min_cut, NULL, GrB_PLUS_MONOID_FP64, cut_set, NULL));
+
+    TEST_CHECK(flow == min_cut);
+    printf("The min cut: %lf\n", min_cut);
+
+    GrB_free(&A);
+    GrB_free(&R);
+    GrB_free(&S);
+    GrB_free(&S_bar);
+    GrB_free(&cut_set);
+    LAGraph_Delete(&G, msg);
+  }
+
+#endif
+  
+
+  LAGraph_Finalize(msg);
+  
+}
+
+TEST_LIST = {{"MinCut", test_MinCut}, {NULL, NULL}};

include/LAGraphX.h

@@ -1507,6 +1507,7 @@ int LAGr_MaxFlow(
     //outputs
     double* f,
     GrB_Matrix* flow_mtx,
+    GrB_Matrix* res_mtx,
     //inputs
     LAGraph_Graph G,
     GrB_Index src, //source node index
@@ -1515,6 +1516,19 @@ int LAGr_MaxFlow(
     char* msg
 );
 
+LAGRAPHX_PUBLIC
+int LAGraph_MinCut(
+    //outputs
+    GrB_Vector* S,
+    GrB_Vector* S_bar,
+    GrB_Matrix* cut_set,
+    // inputs
+    LAGraph_Graph G_origin, //original graph with capacities
+    GrB_Matrix R, //residual graph
+    GrB_Index s, //source node index
+    GrB_Index t, //sink node index
+    char *msg
+);
 
 #if defined ( __cplusplus )
 }