GRASS Programmer's Manual
6.4.2(2012)
|
00001 00017 #include <stdio.h> 00018 #include <stdlib.h> 00019 #include <grass/gis.h> 00020 #include <grass/Vect.h> 00021 #include <grass/glocale.h> 00022 #include <grass/dgl/graph.h> 00023 #include <grass/neta.h> 00024 00025 dglInt32_t sign(dglInt32_t x) 00026 { 00027 if (x >= 0) 00028 return 1; 00029 return -1; 00030 } 00031 00047 int NetA_flow(dglGraph_s * graph, struct ilist *source_list, 00048 struct ilist *sink_list, int *flow) 00049 { 00050 int nnodes, nlines, i; 00051 dglEdgesetTraverser_s et; 00052 dglInt32_t *queue; 00053 dglInt32_t **prev; 00054 char *is_source, *is_sink; 00055 int begin, end, total_flow; 00056 00057 nnodes = dglGet_NodeCount(graph); 00058 nlines = dglGet_EdgeCount(graph) / 2; /*each line corresponds to two edges. One in each direction */ 00059 queue = (dglInt32_t *) G_calloc(nnodes + 3, sizeof(dglInt32_t)); 00060 prev = (dglInt32_t **) G_calloc(nnodes + 3, sizeof(dglInt32_t *)); 00061 is_source = (char *)G_calloc(nnodes + 3, sizeof(char)); 00062 is_sink = (char *)G_calloc(nnodes + 3, sizeof(char)); 00063 if (!queue || !prev || !is_source || !is_sink) { 00064 G_fatal_error(_("Out of memory")); 00065 return -1; 00066 } 00067 00068 for (i = 0; i < source_list->n_values; i++) 00069 is_source[source_list->value[i]] = 1; 00070 for (i = 0; i < sink_list->n_values; i++) 00071 is_sink[sink_list->value[i]] = 1; 00072 00073 for (i = 0; i <= nlines; i++) 00074 flow[i] = 0; 00075 00076 total_flow = 0; 00077 while (1) { 00078 dglInt32_t node, edge_id, min_residue; 00079 int found = -1; 00080 00081 begin = end = 0; 00082 for (i = 0; i < source_list->n_values; i++) 00083 queue[end++] = source_list->value[i]; 00084 00085 for (i = 1; i <= nnodes; i++) { 00086 prev[i] = NULL; 00087 } 00088 while (begin != end && found == -1) { 00089 dglInt32_t vertex = queue[begin++]; 00090 dglInt32_t *edge, *node = dglGetNode(graph, vertex); 00091 00092 dglEdgeset_T_Initialize(&et, graph, 00093 dglNodeGet_OutEdgeset(graph, node)); 00094 for (edge = dglEdgeset_T_First(&et); edge; 00095 edge = dglEdgeset_T_Next(&et)) { 00096 dglInt32_t cap = dglEdgeGet_Cost(graph, edge); 00097 dglInt32_t id = dglEdgeGet_Id(graph, edge); 00098 dglInt32_t to = 00099 dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge)); 00100 if (!is_source[to] && prev[to] == NULL && 00101 cap > sign(id) * flow[abs(id)]) { 00102 prev[to] = edge; 00103 if (is_sink[to]) { 00104 found = to; 00105 break; 00106 } 00107 queue[end++] = to; 00108 } 00109 } 00110 dglEdgeset_T_Release(&et); 00111 } 00112 if (found == -1) 00113 break; /*no augmenting path */ 00114 /*find minimum residual capacity along the augmenting path */ 00115 node = found; 00116 edge_id = dglEdgeGet_Id(graph, prev[node]); 00117 min_residue = 00118 dglEdgeGet_Cost(graph, 00119 prev[node]) - sign(edge_id) * flow[abs(edge_id)]; 00120 while (!is_source[node]) { 00121 dglInt32_t residue; 00122 00123 edge_id = dglEdgeGet_Id(graph, prev[node]); 00124 residue = 00125 dglEdgeGet_Cost(graph, 00126 prev[node]) - 00127 sign(edge_id) * flow[abs(edge_id)]; 00128 if (residue < min_residue) 00129 min_residue = residue; 00130 node = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[node])); 00131 } 00132 total_flow += min_residue; 00133 /*update flow along the augmenting path */ 00134 node = found; 00135 while (!is_source[node]) { 00136 edge_id = dglEdgeGet_Id(graph, prev[node]); 00137 flow[abs(edge_id)] += sign(edge_id) * min_residue; 00138 node = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[node])); 00139 } 00140 } 00141 00142 G_free(queue); 00143 G_free(prev); 00144 G_free(is_source); 00145 G_free(is_sink); 00146 00147 return total_flow; 00148 } 00149 00165 int NetA_min_cut(dglGraph_s * graph, struct ilist *source_list, 00166 struct ilist *sink_list, int *flow, struct ilist *cut) 00167 { 00168 int nnodes, i; 00169 dglEdgesetTraverser_s et; 00170 dglInt32_t *queue; 00171 char *visited; 00172 int begin, end, total_flow; 00173 00174 nnodes = dglGet_NodeCount(graph); 00175 queue = (dglInt32_t *) G_calloc(nnodes + 3, sizeof(dglInt32_t)); 00176 visited = (char *)G_calloc(nnodes + 3, sizeof(char)); 00177 if (!queue || !visited) { 00178 G_fatal_error(_("Out of memory")); 00179 return -1; 00180 } 00181 00182 total_flow = begin = end = 0; 00183 00184 for (i = 1; i <= nnodes; i++) 00185 visited[i] = 0; 00186 00187 for (i = 0; i < source_list->n_values; i++) { 00188 queue[end++] = source_list->value[i]; 00189 visited[source_list->value[i]] = 1; 00190 } 00191 00192 /* find vertices reachable from source(s) using only non-saturated edges */ 00193 while (begin != end) { 00194 dglInt32_t vertex = queue[begin++]; 00195 dglInt32_t *edge, *node = dglGetNode(graph, vertex); 00196 00197 dglEdgeset_T_Initialize(&et, graph, 00198 dglNodeGet_OutEdgeset(graph, node)); 00199 for (edge = dglEdgeset_T_First(&et); edge; 00200 edge = dglEdgeset_T_Next(&et)) { 00201 dglInt32_t cap = dglEdgeGet_Cost(graph, edge); 00202 dglInt32_t id = dglEdgeGet_Id(graph, edge); 00203 dglInt32_t to = 00204 dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge)); 00205 if (!visited[to] && cap > sign(id) * flow[abs(id)]) { 00206 visited[to] = 1; 00207 queue[end++] = to; 00208 } 00209 } 00210 dglEdgeset_T_Release(&et); 00211 } 00212 /*saturated edges from reachable vertices to non-reachable ones form a minimum cost */ 00213 Vect_reset_list(cut); 00214 for (i = 1; i <= nnodes; i++) { 00215 if (!visited[i]) 00216 continue; 00217 dglInt32_t *node, *edgeset, *edge; 00218 00219 node = dglGetNode(graph, i); 00220 edgeset = dglNodeGet_OutEdgeset(graph, node); 00221 dglEdgeset_T_Initialize(&et, graph, edgeset); 00222 for (edge = dglEdgeset_T_First(&et); edge; 00223 edge = dglEdgeset_T_Next(&et)) { 00224 dglInt32_t to, edge_id; 00225 00226 to = dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge)); 00227 edge_id = abs(dglEdgeGet_Id(graph, edge)); 00228 if (!visited[to] && flow[edge_id] != 0) { 00229 Vect_list_append(cut, edge_id); 00230 total_flow += abs(flow[abs(edge_id)]); 00231 } 00232 } 00233 dglEdgeset_T_Release(&et); 00234 } 00235 00236 G_free(visited); 00237 G_free(queue); 00238 return total_flow; 00239 } 00240 00258 int NetA_split_vertices(dglGraph_s * in, dglGraph_s * out, int *node_costs) 00259 { 00260 dglInt32_t opaqueset[16] = 00261 { 360000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; 00262 dglNodeTraverser_s nt; 00263 dglInt32_t nnodes, edge_cnt; 00264 dglInt32_t *cur_node; 00265 00266 nnodes = dglGet_NodeCount(in); 00267 dglInitialize(out, (dglByte_t) 1, (dglInt32_t) 0, (dglInt32_t) 0, 00268 opaqueset); 00269 dglNode_T_Initialize(&nt, in); 00270 edge_cnt = 0; 00271 dglInt32_t max_node_cost = 0; 00272 00273 for (cur_node = dglNode_T_First(&nt); cur_node; 00274 cur_node = dglNode_T_Next(&nt)) { 00275 dglInt32_t v = dglNodeGet_Id(in, cur_node); 00276 00277 edge_cnt++; 00278 dglInt32_t cost = 1; 00279 00280 if (node_costs) 00281 cost = node_costs[v]; 00282 if (cost > max_node_cost) 00283 max_node_cost = cost; 00284 dglAddEdge(out, 2 * v - 1, 2 * v, cost, edge_cnt); 00285 dglAddEdge(out, 2 * v, 2 * v - 1, (dglInt32_t) 0, -edge_cnt); 00286 } 00287 dglNode_T_Release(&nt); 00288 dglNode_T_Initialize(&nt, in); 00289 for (cur_node = dglNode_T_First(&nt); cur_node; 00290 cur_node = dglNode_T_Next(&nt)) { 00291 dglEdgesetTraverser_s et; 00292 dglInt32_t *edge; 00293 dglInt32_t v = dglNodeGet_Id(in, cur_node); 00294 00295 dglEdgeset_T_Initialize(&et, in, dglNodeGet_OutEdgeset(in, cur_node)); 00296 for (edge = dglEdgeset_T_First(&et); edge; 00297 edge = dglEdgeset_T_Next(&et)) { 00298 dglInt32_t to; 00299 00300 to = dglNodeGet_Id(in, dglEdgeGet_Tail(in, edge)); 00301 edge_cnt++; 00302 dglAddEdge(out, 2 * v, 2 * to - 1, max_node_cost + 1, edge_cnt); 00303 dglAddEdge(out, 2 * to - 1, 2 * v, (dglInt32_t) 0, -edge_cnt); 00304 } 00305 dglEdgeset_T_Release(&et); 00306 } 00307 dglNode_T_Release(&nt); 00308 if (dglFlatten(out) < 0) 00309 G_fatal_error(_("GngFlatten error")); 00310 return edge_cnt; 00311 }