1 #include <stdlib.h>
2 #include <string.h>
3 #include <time.h>
4 #include <grass/gis.h>
5 #include <grass/vector.h>
6 #include <grass/dbmi.h>
7 #include <grass/dgl/graph.h>
8 #include <grass/glocale.h>
9 #include "alloc.h"
10
alloc_from_centers_loop_tt(struct Map_info * Map,NODE * Nodes,CENTER * Centers,int ncenters,int tucfield)11 int alloc_from_centers_loop_tt(struct Map_info *Map, NODE *Nodes,
12 CENTER *Centers, int ncenters,
13 int tucfield)
14 {
15 int center, line, nlines, i;
16 int node1;
17 int cat;
18 struct line_cats *Cats;
19 struct line_pnts *Points;
20 double cost, n1cost, n2cost;
21 int ret;
22
23 Cats = Vect_new_cats_struct();
24 Points = Vect_new_line_struct();
25
26 nlines = Vect_get_num_lines(Map);
27
28 for (i = 2; i <= (nlines * 2 + 2); i++) {
29 Nodes[i].center = -1;/* NOTE: first two items of Nodes are not used */
30 Nodes[i].cost = -1;
31 Nodes[i].edge = 0;
32 }
33
34 for (center = 0; center < ncenters; center++) {
35 G_percent(center, ncenters, 1);
36 node1 = Centers[center].node;
37 Vect_net_get_node_cost(Map, node1, &n1cost);
38 G_debug(2, "center = %d node = %d cat = %d", center, node1,
39 Centers[center].cat);
40
41 for (line = 1; line <= nlines; line++) {
42 G_debug(5, " node1 = %d line = %d", node1, line);
43 Vect_net_get_node_cost(Map, line, &n2cost);
44 /* closed, left it as not attached */
45
46 if (Vect_read_line(Map, Points, Cats, line) < 0)
47 continue;
48 if (Vect_get_line_type(Map, line) != GV_LINE)
49 continue;
50 if (!Vect_cat_get(Cats, tucfield, &cat))
51 continue;
52
53 for (i = 0; i < 2; i++) {
54 if (i == 1)
55 cat *= -1;
56
57 ret =
58 Vect_net_ttb_shortest_path(Map, node1, 0, cat, 1,
59 tucfield, NULL,
60 &cost);
61 if (ret == -1) {
62 continue;
63 } /* node unreachable */
64
65 /* We must add center node costs (not calculated by Vect_net_shortest_path() ), but
66 * only if center and node are not identical, because at the end node cost is add later */
67 if (ret != 1)
68 cost += n1cost;
69
70 G_debug(5,
71 "Arc nodes: %d %d cost: %f (x old cent: %d old cost %f",
72 node1, line, cost, Nodes[line * 2 + i].center,
73 Nodes[line * 2 + i].cost);
74 if (Nodes[line * 2 + i].center == -1 ||
75 Nodes[line * 2 + i].cost > cost) {
76 Nodes[line * 2 + i].cost = cost;
77 Nodes[line * 2 + i].center = center;
78 }
79 }
80 }
81 }
82 G_percent(1, 1, 1);
83
84 Vect_destroy_cats_struct(Cats);
85 Vect_destroy_line_struct(Points);
86
87 return 0;
88 }
89
alloc_to_centers_loop_tt(struct Map_info * Map,NODE * Nodes,CENTER * Centers,int ncenters,int tucfield)90 int alloc_to_centers_loop_tt(struct Map_info *Map, NODE *Nodes,
91 CENTER *Centers, int ncenters,
92 int tucfield)
93 {
94 int center, line, nlines, i;
95 int node1;
96 int cat;
97 struct line_cats *Cats;
98 struct line_pnts *Points;
99 double cost, n1cost, n2cost;
100 int ret;
101
102 Cats = Vect_new_cats_struct();
103 Points = Vect_new_line_struct();
104
105 nlines = Vect_get_num_lines(Map);
106
107 for (i = 2; i <= (nlines * 2 + 2); i++) {
108 Nodes[i].center = -1;/* NOTE: first two items of Nodes are not used */
109 Nodes[i].cost = -1;
110 Nodes[i].edge = 0;
111 }
112
113 for (line = 1; line <= nlines; line++) {
114 G_debug(5, " line = %d", line);
115 Vect_net_get_node_cost(Map, line, &n2cost);
116 /* closed, left it as not attached */
117
118 if (Vect_read_line(Map, Points, Cats, line) < 0)
119 continue;
120 if (Vect_get_line_type(Map, line) != GV_LINE)
121 continue;
122 if (!Vect_cat_get(Cats, tucfield, &cat))
123 continue;
124
125 for (center = 0; center < ncenters; center++) {
126 G_percent(center, ncenters, 1);
127 node1 = Centers[center].node;
128 Vect_net_get_node_cost(Map, node1, &n1cost);
129 G_debug(2, "center = %d node = %d cat = %d", center, node1,
130 Centers[center].cat);
131
132 for (i = 0; i < 2; i++) {
133 if (i == 1)
134 cat *= -1;
135
136 ret =
137 Vect_net_ttb_shortest_path(Map, cat, 1, node1, 0,
138 tucfield, NULL,
139 &cost);
140 if (ret == -1) {
141 continue;
142 } /* node unreachable */
143
144 /* We must add center node costs (not calculated by Vect_net_shortest_path() ), but
145 * only if center and node are not identical, because at the end node cost is add later */
146 if (ret != 1)
147 cost += n1cost;
148
149 G_debug(5,
150 "Arc nodes: %d %d cost: %f (x old cent: %d old cost %f",
151 node1, line, cost, Nodes[line * 2 + i].center,
152 Nodes[line * 2 + i].cost);
153 if (Nodes[line * 2 + i].center == -1 ||
154 Nodes[line * 2 + i].cost > cost) {
155 Nodes[line * 2 + i].cost = cost;
156 Nodes[line * 2 + i].center = center;
157 }
158 }
159 }
160 }
161 G_percent(1, 1, 1);
162
163 Vect_destroy_cats_struct(Cats);
164 Vect_destroy_line_struct(Points);
165
166 return 0;
167 }
168
alloc_from_centers(dglGraph_s * graph,NODE * Nodes,CENTER * Centers,int ncenters)169 int alloc_from_centers(dglGraph_s *graph, NODE *Nodes, CENTER *Centers, int ncenters)
170 {
171 int i, nnodes;
172 dglHeap_s heap;
173 dglEdgesetTraverser_s et;
174 int have_node_costs;
175 dglInt32_t ncost;
176
177 nnodes = dglGet_NodeCount(graph);
178
179 /* initialize nodes */
180 for (i = 1; i <= nnodes; i++) {
181 Nodes[i].cost = -1;
182 Nodes[i].center = -1;
183 Nodes[i].edge = 0;
184 }
185
186 ncost = 0;
187 have_node_costs = dglGet_NodeAttrSize(graph);
188
189 dglHeapInit(&heap);
190
191 for (i = 0; i < ncenters; i++) {
192 int v = Centers[i].node;
193
194 if (Nodes[v].cost == 0)
195 continue; /* ignore duplicates */
196 Nodes[v].cost = 0; /* make sure all centers are processed first */
197 Nodes[v].center = i;
198 dglHeapData_u heap_data;
199
200 heap_data.ul = v;
201 dglHeapInsertMin(&heap, 0, ' ', heap_data);
202 }
203 while (1) {
204 dglInt32_t v, dist;
205 dglHeapNode_s heap_node;
206 dglHeapData_u heap_data;
207 dglInt32_t *edge;
208 dglInt32_t *node;
209
210 if (!dglHeapExtractMin(&heap, &heap_node))
211 break;
212 v = heap_node.value.ul;
213 dist = heap_node.key;
214 if (Nodes[v].cost < dist)
215 continue;
216
217 node = dglGetNode(graph, v);
218
219 if (have_node_costs && Nodes[v].edge) {
220 memcpy(&ncost, dglNodeGet_Attr(graph, node),
221 sizeof(ncost));
222 if (ncost > 0)
223 dist += ncost;
224 /* do not go through closed nodes */
225 if (ncost < 0)
226 continue;
227 }
228
229 dglEdgeset_T_Initialize(&et, graph,
230 dglNodeGet_OutEdgeset(graph, node));
231
232 for (edge = dglEdgeset_T_First(&et); edge;
233 edge = dglEdgeset_T_Next(&et)) {
234 dglInt32_t *to = dglEdgeGet_Tail(graph, edge);
235 dglInt32_t to_id = dglNodeGet_Id(graph, to);
236 dglInt32_t d = dglEdgeGet_Cost(graph, edge);
237
238 if (Nodes[to_id].cost < 0 || Nodes[to_id].cost > dist + d) {
239 Nodes[to_id].cost = dist + d;
240 Nodes[to_id].edge = dglEdgeGet_Id(graph, edge);
241 Nodes[to_id].center = Nodes[v].center;
242 heap_data.ul = to_id;
243 dglHeapInsertMin(&heap, dist + d, ' ', heap_data);
244 }
245 }
246
247 dglEdgeset_T_Release(&et);
248 }
249
250 dglHeapFree(&heap, NULL);
251
252 return 0;
253 }
254
alloc_to_centers(dglGraph_s * graph,NODE * Nodes,CENTER * Centers,int ncenters)255 int alloc_to_centers(dglGraph_s *graph, NODE *Nodes, CENTER *Centers, int ncenters)
256 {
257 int i, nnodes;
258 dglHeap_s heap;
259 dglEdgesetTraverser_s et;
260 int have_node_costs;
261 dglInt32_t ncost;
262
263 if (graph->Version < 2) {
264 G_warning("Directed graph must be version 2 or 3 for distances to centers");
265 return -1;
266 }
267
268 nnodes = dglGet_NodeCount(graph);
269
270 /* initialize nodes */
271 for (i = 1; i <= nnodes; i++) {
272 Nodes[i].cost = -1;
273 Nodes[i].center = -1;
274 Nodes[i].edge = 0;
275 }
276
277 ncost = 0;
278 have_node_costs = dglGet_NodeAttrSize(graph);
279
280 dglHeapInit(&heap);
281
282 for (i = 0; i < ncenters; i++) {
283 int v = Centers[i].node;
284
285 if (Nodes[v].cost == 0)
286 continue; /* ignore duplicates */
287 Nodes[v].cost = 0; /* make sure all centers are processed first */
288 Nodes[v].center = i;
289 dglHeapData_u heap_data;
290
291 heap_data.ul = v;
292 dglHeapInsertMin(&heap, 0, ' ', heap_data);
293 }
294 while (1) {
295 dglInt32_t v, dist;
296 dglHeapNode_s heap_node;
297 dglHeapData_u heap_data;
298 dglInt32_t *edge;
299 dglInt32_t *node;
300
301 if (!dglHeapExtractMin(&heap, &heap_node))
302 break;
303 v = heap_node.value.ul;
304 dist = heap_node.key;
305 if (Nodes[v].cost < dist)
306 continue;
307
308 node = dglGetNode(graph, v);
309
310 if (have_node_costs && Nodes[v].edge) {
311 memcpy(&ncost, dglNodeGet_Attr(graph, node),
312 sizeof(ncost));
313 if (ncost > 0)
314 dist += ncost;
315 /* do not go through closed nodes */
316 if (ncost < 0)
317 continue;
318 }
319
320 dglEdgeset_T_Initialize(&et, graph,
321 dglNodeGet_InEdgeset(graph, node));
322
323 for (edge = dglEdgeset_T_First(&et); edge;
324 edge = dglEdgeset_T_Next(&et)) {
325 dglInt32_t *from = dglEdgeGet_Head(graph, edge);
326 dglInt32_t from_id = dglNodeGet_Id(graph, from);
327 dglInt32_t d = dglEdgeGet_Cost(graph, edge);
328
329 if (Nodes[from_id].cost < 0 || Nodes[from_id].cost > dist + d) {
330 Nodes[from_id].cost = dist + d;
331 Nodes[from_id].edge = dglEdgeGet_Id(graph, edge);
332 Nodes[from_id].center = Nodes[v].center;
333 heap_data.ul = from_id;
334 dglHeapInsertMin(&heap, dist + d, ' ', heap_data);
335 }
336 }
337
338 dglEdgeset_T_Release(&et);
339 }
340
341 dglHeapFree(&heap, NULL);
342
343 return 0;
344 }
345
346