Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_SwapEdges: randomly swaps edges in a graph
3 : //------------------------------------------------------------------------------
4 :
5 : // LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved.
6 : // SPDX-License-Identifier: BSD-2-Clause
7 : //
8 : // For additional details (including references to third party source code and
9 : // other files) see the LICENSE file or contact permission@sei.cmu.edu. See
10 : // Contributors.txt for a full list of contributors. Created, in part, with
11 : // funding and support from the U.S. Government (see Acknowledgments.txt file).
12 : // DM22-0790
13 :
14 : // Contributed by Gabriel Gomez, Texas A&M University
15 :
16 : //------------------------------------------------------------------------------
17 :
18 : // References:
19 :
20 : // R. Milo, N. Kashtan, S. Itzkovitz, M. E. J. Newman, and U. Alon, “On the
21 : // uniform generation of random graphs with prescribed degree sequences,” 2004.
22 :
23 : #define FREE_LOOP \
24 : { \
25 : GrB_free (&M) ; \
26 : GrB_free(&dup_swaps_v); \
27 : GrB_free(&new_hashed_edges); \
28 : GrB_free(&edge_perm) ; \
29 : }
30 :
31 : #define LG_FREE_WORK \
32 : { \
33 : /* free any workspace used here */ \
34 : GrB_free (&A_tril) ; \
35 : GrB_free (&Ai) ; \
36 : GrB_free (&Aj) ; \
37 : GrB_free (&random_v) ; \
38 : GrB_free (&r_permute) ; \
39 : GrB_free (&ramp_v) ; \
40 : GrB_free (&E_temp) ; \
41 : GrB_free (&r_60) ; \
42 : GrB_free (&exists) ; \
43 : GrB_free (&E_vec) ; \
44 : GrB_free (&swap_pair) ; \
45 : GrB_free (&hash_seed_e) ; \
46 : GrB_free (&add_term_biop) ; \
47 : GrB_free (&add_term_monoid) ; \
48 : GrB_free (&plus_term_one) ; \
49 : GrB_free (&second_edge) ; \
50 : GrB_free (&second_bool_edge) ; \
51 : GrB_free (&second_edge_monoid) ; \
52 : GrB_free (&second_second_edge) ; \
53 : GrB_free (&lg_shiftland) ; \
54 : GrB_free (&lg_edge) ; \
55 : GrB_free (&lg_swap) ; \
56 : GrB_free (&hashed_edges) ; \
57 : GrB_free (&con) ; \
58 : GrB_free (&one8) ; \
59 : GrB_free (&x) ; \
60 : GrB_free (&temp_i) ; \
61 : LAGraph_Free ((void ** )&indices, NULL) ; \
62 : LAGraph_Free ((void **) &dup_swaps, NULL); \
63 : FREE_LOOP ; \
64 : }
65 :
66 : #define LG_FREE_ALL \
67 : { \
68 : /* free any workspace used here */ \
69 : LG_FREE_WORK ; \
70 : /* free all the output variable(s) */ \
71 : GrB_free (&A_new) ; \
72 : LAGraph_Delete(G_new, NULL) ; \
73 : /* take any other corrective action */ \
74 : }
75 :
76 : #include "LG_internal.h"
77 : #include "LAGraphX.h"
78 14091831 : void shift_and
79 : (uint16_t *z, const uint16_t *x)
80 : {
81 14091831 : (*z) = (*x) & ((*x) << 8);
82 14091831 : (*z) |= (*z) >> 8;
83 14091831 : }
84 : #define SHIFT_AND \
85 : "void shift_and \n"\
86 : " (uint16_t *z, const uint16_t *x) \n"\
87 : " { \n"\
88 : " (*z) = (*x) & ((*x) << 8); \n"\
89 : " (*z) |= (*z) >> 8; \n"\
90 : " }"
91 :
92 : typedef struct {
93 : uint64_t a;
94 : uint64_t b;
95 : } edge_type64;
96 : #define EDGE_TYPE64 \
97 : "typedef struct { uint64_t a; uint64_t b; } edge_type64;"
98 :
99 : typedef struct {
100 : uint64_t a;
101 : uint64_t b;
102 : uint64_t c;
103 : uint64_t d;
104 : } swap_type64;
105 : #define SWAP_TYPE64 \
106 : "typedef struct { \n"\
107 : " uint64_t a; uint64_t b; uint64_t c; uint64_t d; \n"\
108 : "} swap_type64;"
109 :
110 : typedef struct {
111 : uint32_t a;
112 : uint32_t b;
113 : } edge_type32;
114 : #define EDGE_TYPE32 \
115 : "typedef struct { uint32_t a; uint32_t b; } edge_type32;"
116 :
117 : typedef struct {
118 : uint32_t a;
119 : uint32_t b;
120 : uint32_t c;
121 : uint32_t d;
122 : } swap_type32;
123 : #define SWAP_TYPE32 \
124 : "typedef struct { \n"\
125 : " uint32_t a; uint32_t b; uint32_t c; uint32_t d; \n"\
126 : "}swap_type32;"
127 :
128 14506920 : void swap_bc64
129 : (swap_type64 *z, const swap_type64 *x, GrB_Index I, GrB_Index J, const bool *y)
130 : {
131 14506920 : memcpy(z, x, sizeof(*z)) ; //unnessesary when aliassed but done for safety.
132 14506920 : if(z->a == z->c || z->b == z->c || z->a == z->d || z->b == z->d ) return;
133 14432826 : if(I & 1)
134 : {
135 7209524 : uint64_t temp = z->d;
136 7209524 : z->d = z->b;
137 7209524 : z->b = temp;
138 : }
139 : else
140 : {
141 7223302 : uint64_t temp = z->c;
142 7223302 : z->c = z->b;
143 7223302 : z->b = temp;
144 : }
145 : }
146 276453 : void swap_bc32
147 : (swap_type32 *z, const swap_type32 *x, GrB_Index I, GrB_Index J, const bool *y)
148 : {
149 276453 : memcpy(z, x, sizeof(*z)) ; //unnessesary when aliassed but done for safety.
150 276453 : if(z->a == z->c || z->b == z->c || z->a == z->d || z->b == z->d ) return;
151 259407 : if(I & 1)
152 : {
153 131500 : uint32_t temp = z->d;
154 131500 : z->d = z->b;
155 131500 : z->b = temp;
156 : }
157 : else
158 : {
159 127907 : uint32_t temp = z->c;
160 127907 : z->c = z->b;
161 127907 : z->b = temp;
162 : }
163 : }
164 : #define SWAP_BC64 \
165 : "void swap_bc64 \n"\
166 : "(swap_type64 *z, const swap_type64 *x, GrB_Index I, GrB_Index J, const bool *y)\n"\
167 : "{ \n"\
168 : " memcpy(z, x, sizeof(*z)) ; //unnessesary when aliassed but done for safety.\n"\
169 : " if(z->a == z->c || z->b == z->c || z->a == z->d || z->b == z->d ) return; \n"\
170 : " if(I & 1) \n"\
171 : " { \n"\
172 : " uint64_t temp = z->d; \n"\
173 : " z->d = z->b; \n"\
174 : " z->b = temp; \n"\
175 : " } \n"\
176 : " else \n"\
177 : " { \n"\
178 : " uint64_t temp = z->c; \n"\
179 : " z->c = z->b; \n"\
180 : " z->b = temp; \n"\
181 : " } \n"\
182 : "}"
183 : #define SWAP_BC32 \
184 : "void swap_bc32 \n"\
185 : "(swap_type32 *z, const swap_type32 *x, GrB_Index I, GrB_Index J, const bool *y)\n"\
186 : "{ \n"\
187 : " memcpy(z, x, sizeof(*z)) ; //unnessesary when aliassed but done for safety.\n"\
188 : " if(z->a == z->c || z->b == z->c || z->a == z->d || z->b == z->d ) return; \n"\
189 : " if(I & 1) \n"\
190 : " { \n"\
191 : " uint32_t temp = z->d; \n"\
192 : " z->d = z->b; \n"\
193 : " z->b = temp; \n"\
194 : " } \n"\
195 : " else \n"\
196 : " { \n"\
197 : " uint32_t temp = z->c; \n"\
198 : " z->c = z->b; \n"\
199 : " z->b = temp; \n"\
200 : " } \n"\
201 : "}"
202 :
203 : // using xorshift, from https://en.wikipedia.org/wiki/Xorshift
204 : // with a state of uint64_t, or xorshift64star.
205 71134400 : void hash_edge64
206 : (uint64_t *z, const edge_type64 *x, const uint64_t *mask)
207 : {
208 71134400 : (*z) = x->a ^ x->b;
209 71134400 : (*z) ^= (*z) << 13;
210 71134400 : (*z) ^= (*z) >> 7;
211 71134400 : (*z) ^= (x->a < x->b)? x->a: x->b;
212 71134400 : (*z) ^= (*z) << 17;
213 71134400 : (*z) &= (*mask);
214 71134400 : }
215 1376528 : void hash_edge32
216 : (uint64_t *z, const edge_type32 *x, const uint64_t *mask)
217 : {
218 1376528 : (*z) = x->a ^ x->b;
219 1376528 : (*z) ^= (*z) << 13;
220 1376528 : (*z) ^= (*z) >> 7;
221 1376528 : (*z) ^= (uint64_t)((x->a < x->b)? x->a: x->b);
222 1376528 : (*z) ^= (*z) << 17;
223 1376528 : (*z) &= (*mask);
224 1376528 : }
225 : #define HASH_EDGE64 \
226 : "void hash_edge64 \n"\
227 : "(uint64_t *z, const edge_type64 *x, const uint64_t *mask) \n"\
228 : "{ \n"\
229 : " (*z) = x->a ^ x->b; \n"\
230 : " (*z) ^= (*z) << 13; \n"\
231 : " (*z) ^= (*z) >> 7; \n"\
232 : " (*z) ^= (uint64_t)((x->a < x->b)? x->a: x->b); \n"\
233 : " (*z) ^= (*z) << 17; \n"\
234 : " (*z) &= (*mask); \n"\
235 : "}"
236 : #define HASH_EDGE32 \
237 : "void hash_edge32 \n"\
238 : "(uint64_t *z, const edge_type32 *x, const uint64_t *mask) \n"\
239 : "{ \n"\
240 : " (*z) = x->a ^ x->b; \n"\
241 : " (*z) ^= (*z) << 13; \n"\
242 : " (*z) ^= (*z) >> 7; \n"\
243 : " (*z) ^= (uint64_t)((x->a < x->b)? x->a: x->b); \n"\
244 : " (*z) ^= (*z) << 17; \n"\
245 : " (*z) &= (*mask); \n"\
246 : "}"
247 :
248 28599930 : void add_term
249 : (int8_t *z, const int8_t *x, const int8_t *y)
250 : {
251 28599930 : (*z) = (*x) | (*y) + ((int8_t)1 & (*x) & (*y)) ;
252 28599930 : }
253 : #define ADD_TERM \
254 : "void add_term \n"\
255 : "(int8_t *z, const int8_t *x, const int8_t *y) \n"\
256 : "{ \n"\
257 : " (*z) = (*x) | (*y) + ((int8_t)1 & (*x) & (*y)) ; \n"\
258 : "}"
259 :
260 12637456 : void edge2nd64_bool
261 : (edge_type64 *z, const bool *x, const edge_type64 *y)
262 : {
263 12637456 : z->a = y->a;
264 12637456 : z->b = y->b;
265 12637456 : }
266 219630 : void edge2nd32_bool
267 : (edge_type32 *z, const bool *x, const edge_type32 *y)
268 : {
269 219630 : z->a = y->a;
270 219630 : z->b = y->b;
271 219630 : }
272 12637456 : void edge2nd64_edge
273 : (edge_type64 *z, const edge_type64 *x, const edge_type64 *y)
274 : {
275 12637456 : z->a = y->a;
276 12637456 : z->b = y->b;
277 12637456 : }
278 219630 : void edge2nd32_edge
279 : (edge_type32 *z, const edge_type32 *x, const edge_type32 *y)
280 : {
281 219630 : z->a = y->a;
282 219630 : z->b = y->b;
283 219630 : }
284 : #define EDGE2ND32_BOOL \
285 : "void edge2nd32_bool \n"\
286 : "(edge_type32 *z, const bool *x, const edge_type32 *y) \n"\
287 : "{ \n"\
288 : " //if(y->a == 0 && y->b == 0) return; \n"\
289 : " z->a = y->a; \n"\
290 : " z->b = y->b; \n"\
291 : "}"
292 : #define EDGE2ND64_BOOL \
293 : "void edge2nd64_bool \n"\
294 : "(edge_type64 *z, const bool *x, const edge_type64 *y) \n"\
295 : "{ \n"\
296 : " //if(y->a == 0 && y->b == 0) return; \n"\
297 : " z->a = y->a; \n"\
298 : " z->b = y->b; \n"\
299 : "}"
300 : #define EDGE2ND32_EDGE \
301 : "void edge2nd32_edge \n"\
302 : "(edge_type32 *z, const edge_type32 *x, const edge_type32 *y) \n"\
303 : "{ \n"\
304 : " //if(y->a == 0 && y->b == 0) return; \n"\
305 : " z->a = y->a; \n"\
306 : " z->b = y->b; \n"\
307 : "}"
308 : #define EDGE2ND64_EDGE \
309 : "void edge2nd64_edge \n"\
310 : "(edge_type64 *z, const edge_type64 *x, const edge_type64 *y) \n"\
311 : "{ \n"\
312 : " //if(y->a == 0 && y->b == 0) return; \n"\
313 : " z->a = y->a; \n"\
314 : " z->b = y->b; \n"\
315 : "}"
316 :
317 5245 : int LAGr_SwapEdges
318 : (
319 : // output
320 : LAGraph_Graph *G_new, // A new graph with the same degree for each node
321 : uint64_t *pSwaps, // Actual number of Swaps proformed
322 : // input: not modified
323 : const LAGraph_Graph G, // Graph to be randomized.
324 : double loopTry, // Percent of edges to involve per loop [0,1]
325 : double loopMin, // Minimum Swaps percent per loop [0,1)
326 : uint64_t totSwaps, // Desired Swaps
327 : uint64_t seed, // Random Seed
328 : char *msg
329 : )
330 : {
331 : #if LG_SUITESPARSE_GRAPHBLAS_V10
332 : //--------------------------------------------------------------------------
333 : // Declorations
334 : //--------------------------------------------------------------------------
335 5245 : GrB_Matrix A = NULL, A_new = NULL; // n x n Adjacency Matrix
336 5245 : GrB_Vector Ai = NULL, Aj = NULL;
337 : // e x 1 vector, each entry is an edge.
338 5245 : GrB_Vector E_vec = NULL, E_temp = NULL;
339 :
340 : // swaps x 4
341 : // Each row contains 4 entries corresponding to the verticies
342 : // that are involved in the swap.
343 5245 : GrB_Vector M = NULL;
344 5245 : GxB_Container con = NULL;
345 :
346 : // n = |V| e = |E|
347 5245 : GrB_Index n = 0, e = 0;
348 :
349 : // n x n
350 : // Lower triangle of adjacency matrix
351 5245 : GrB_Matrix A_tril = NULL ;
352 :
353 : // e x 1 random vectors
354 5245 : GrB_Vector random_v = NULL, r_permute = NULL;
355 :
356 : // indicies for A
357 5245 : void *indices = NULL;
358 :
359 5245 : GrB_Vector ramp_v = NULL;
360 :
361 : // edge permutation
362 5245 : GrB_Vector edge_perm = NULL;
363 5245 : bool iso = false;
364 :
365 : // Number of values kept in each phase
366 5245 : GrB_Index n_keep = 0, M_nvals = 0, dup_arr_size = 0;
367 :
368 : // swaps x 2 matrix which holds the hashes of each planned edge.
369 5245 : GrB_Vector new_hashed_edges = NULL;
370 :
371 : // e holds hashes of old edges
372 5245 : GrB_Vector hashed_edges = NULL;
373 :
374 : // 2^60 holds the buckets in which hashes collided.
375 5245 : GrB_Vector exists = NULL;
376 :
377 5245 : GrB_UnaryOp lg_shiftland = NULL;
378 :
379 : // b1 <---> a2 or b1 <---> b2
380 5245 : GrB_IndexUnaryOp swap_pair = NULL;
381 :
382 : // z = h_y(x)
383 5245 : GrB_BinaryOp hash_seed_e = NULL;
384 :
385 : // This monoid has only been designed for inputs in {0,1,2}, other behavior
386 : // is undefined.
387 : // (0,x) -> x, (1,1) -> 2, (2,x) -> 2 (and commutative)
388 : // Aka z = min(2, x + y)
389 5245 : GrB_BinaryOp add_term_biop = NULL;
390 5245 : GrB_Monoid add_term_monoid = NULL;
391 5245 : GrB_Semiring plus_term_one = NULL;
392 :
393 : // z = y
394 5245 : GrB_BinaryOp second_edge = NULL;
395 5245 : GrB_BinaryOp second_bool_edge = NULL;
396 5245 : GrB_Monoid second_edge_monoid = NULL;
397 5245 : GrB_Semiring second_second_edge = NULL;
398 :
399 : // Toople types
400 5245 : GrB_Type lg_edge = NULL, lg_swap = NULL;
401 : // Unload types
402 5245 : GrB_Type M_type = NULL, E_type = NULL, Ai_type = NULL, Aj_type = NULL;
403 5245 : int M_hand = 0, E_hand = 0;
404 :
405 5245 : int16_t *dup_swaps = NULL;
406 5245 : GrB_Vector dup_swaps_v = NULL;
407 : // BOOL swaps * 2 vector that holds false if an edge in the swap did not "work"
408 :
409 5245 : GrB_Vector sort_h = NULL;
410 5245 : GrB_Vector r_60 = NULL;
411 : // Only Used for coverage typecasting
412 5245 : GrB_Vector temp_i = NULL;
413 : // count swaps
414 5245 : GrB_Index num_swaps = 0, num_attempts = 0;
415 :
416 : // Constants ---------------------------------------------------------------
417 5245 : GrB_Vector x = NULL;
418 :
419 5245 : GrB_Scalar one8 = NULL;
420 :
421 5245 : GrB_Index ind_size = 0;
422 :
423 : //--------------------------------------------------------------------------
424 : // Check inputs
425 : //--------------------------------------------------------------------------
426 5245 : LG_ASSERT_MSG (
427 : G->kind == LAGraph_ADJACENCY_UNDIRECTED,
428 : LAGRAPH_INVALID_GRAPH,
429 : "G must be undirected"
430 : ) ;
431 : // char type[LAGRAPH_MAX_NAME_LEN];
432 5245 : LG_ASSERT_MSG (G->nself_edges == 0, LAGRAPH_NO_SELF_EDGES_ALLOWED,
433 : "G->nself_edges must be zero") ;
434 5245 : LG_ASSERT (G_new != NULL, GrB_NULL_POINTER) ;
435 5245 : LG_ASSERT (pSwaps != NULL, GrB_NULL_POINTER) ;
436 5245 : *G_new = NULL ;
437 :
438 : //--------------------------------------------------------------------------
439 : // Initializations
440 : //--------------------------------------------------------------------------
441 5245 : A = G->A ;
442 :
443 : //--------------------------------------------------------------------------
444 : // Extract edges from A
445 : //--------------------------------------------------------------------------
446 :
447 5245 : GRB_TRY (GrB_Matrix_nrows (&n, A)) ;
448 5245 : GRB_TRY (GrB_Matrix_new (&A_tril, GrB_BOOL, n, n)) ;
449 5230 : GRB_TRY (GrB_Vector_new(&Ai, GrB_BOOL, 0)) ;
450 5220 : GRB_TRY (GrB_Vector_new(&Aj, GrB_BOOL, 0)) ;
451 :
452 : // Extract lower triangular edges
453 5210 : GRB_TRY (GrB_select (A_tril, NULL, NULL, GrB_TRIL, A, 0, NULL)) ;
454 5180 : GRB_TRY (GrB_Matrix_nvals(&e, A_tril)) ;
455 5180 : GRB_TRY (GxB_Matrix_extractTuples_Vector(Ai, Aj, NULL, A_tril, NULL)) ;
456 : #ifdef COVERAGE
457 5170 : if(n > 100)
458 : {
459 : // Make Ai and Aj 64 bit
460 2977 : GRB_TRY (GrB_Vector_new(&temp_i, GrB_INT64, e)) ;
461 2971 : GRB_TRY (GrB_assign(temp_i, NULL, NULL, Ai, GrB_ALL, 0, NULL)) ;
462 2968 : GRB_TRY (GrB_free(&Ai)) ;
463 2968 : Ai = temp_i;
464 2968 : temp_i = NULL;
465 2968 : GRB_TRY (GrB_Vector_new(&temp_i, GrB_INT64, e)) ;
466 2962 : GRB_TRY (GrB_assign(temp_i, NULL, NULL, Aj, GrB_ALL, 0, NULL)) ;
467 2959 : GRB_TRY (GrB_free(&Aj)) ;
468 2959 : Aj = temp_i;
469 2959 : temp_i = NULL;
470 : }
471 : #endif
472 5152 : int codei = 0, codej = 0;
473 5152 : GRB_TRY (GxB_Vector_type(&Ai_type, Ai)) ;
474 5152 : GrB_get(Ai, &codei, GrB_EL_TYPE_CODE);
475 5152 : GrB_get(Aj, &codej, GrB_EL_TYPE_CODE);
476 :
477 5152 : LG_ASSERT_MSG (
478 : codei == codej, GrB_INVALID_VALUE,
479 : "extractTuples_Vector returned different types for Ai and Aj"
480 : ) ;
481 : //--------------------------------------------------------------------------
482 : // Initialize all operators and types
483 : //--------------------------------------------------------------------------
484 5152 : if(codei == GrB_UINT32_CODE) // Use uint32 if possible
485 : {
486 2193 : GRB_TRY (GxB_Type_new(
487 : &lg_edge, sizeof(edge_type32), "edge_type32", EDGE_TYPE32)) ;
488 2189 : GRB_TRY (GxB_Type_new(
489 : &lg_swap, sizeof(swap_type32), "swap_type32", SWAP_TYPE32)) ;
490 2185 : GRB_TRY(GxB_BinaryOp_new(
491 : &hash_seed_e, (GxB_binary_function) (&hash_edge32),
492 : GrB_UINT64, lg_edge, GrB_UINT64, "hash_edge32", HASH_EDGE32
493 : )) ;
494 2181 : GRB_TRY (GxB_IndexUnaryOp_new (
495 : &swap_pair, (GxB_index_unary_function) (&swap_bc32),
496 : lg_swap, lg_swap, GrB_BOOL, "swap_bc32", SWAP_BC32
497 : )) ;
498 2177 : GRB_TRY(GxB_BinaryOp_new(
499 : &second_edge, (GxB_binary_function) (&edge2nd32_edge),
500 : lg_edge, lg_edge, lg_edge, "edge2nd32_edge", EDGE2ND32_EDGE
501 : )) ;
502 2173 : GRB_TRY(GxB_BinaryOp_new(
503 : &second_bool_edge, (GxB_binary_function) (&edge2nd32_bool),
504 : lg_edge, GrB_BOOL, lg_edge, "edge2nd32_bool", EDGE2ND32_BOOL
505 : )) ;
506 : }
507 : else //uint64 types
508 : {
509 2959 : GRB_TRY (GxB_Type_new(
510 : &lg_edge, sizeof(edge_type64), "edge_type64", EDGE_TYPE64)) ;
511 2953 : GRB_TRY (GxB_Type_new(
512 : &lg_swap, sizeof(swap_type64), "swap_type64", SWAP_TYPE64)) ;
513 2947 : GRB_TRY(GxB_BinaryOp_new(
514 : &hash_seed_e, (GxB_binary_function) (&hash_edge64),
515 : GrB_UINT64, lg_edge, GrB_UINT64, "hash_edge64", HASH_EDGE64
516 : )) ;
517 2941 : GRB_TRY (GxB_IndexUnaryOp_new (
518 : &swap_pair, (GxB_index_unary_function) (&swap_bc64),
519 : lg_swap, lg_swap, GrB_BOOL, "swap_bc64", SWAP_BC64
520 : )) ;
521 2935 : GRB_TRY(GxB_BinaryOp_new(
522 : &second_edge, (GxB_binary_function) (&edge2nd64_edge),
523 : lg_edge, lg_edge, lg_edge, "edge2nd64_edge", EDGE2ND64_EDGE
524 : )) ;
525 2929 : GRB_TRY(GxB_BinaryOp_new(
526 : &second_bool_edge, (GxB_binary_function) (&edge2nd64_bool),
527 : lg_edge, GrB_BOOL, lg_edge, "edge2nd64_bool", EDGE2ND64_BOOL
528 : )) ;
529 : }
530 :
531 5092 : GRB_TRY (GxB_UnaryOp_new (
532 : &lg_shiftland, (GxB_unary_function) (&shift_and),
533 : GrB_UINT16, GrB_UINT16, "shift_and", SHIFT_AND
534 : )) ;
535 5082 : GRB_TRY(GxB_BinaryOp_new(
536 : &add_term_biop, (GxB_binary_function) (&add_term),
537 : GrB_INT8, GrB_INT8, GrB_INT8, "add_term", ADD_TERM
538 : )) ;
539 :
540 5072 : GRB_TRY (GxB_Monoid_terminal_new_INT8(
541 : &add_term_monoid, add_term_biop, (int8_t) 0, (int8_t) 2
542 : )) ;
543 :
544 5057 : edge_type64 iden_second = {0,0};
545 5057 : GRB_TRY (GrB_Monoid_new_UDT(
546 : &second_edge_monoid, second_edge, (void *) &iden_second
547 : )) ;
548 :
549 5047 : GRB_TRY(GrB_Semiring_new(
550 : &plus_term_one, add_term_monoid, GrB_ONEB_INT8
551 : )) ;
552 5037 : GRB_TRY(GrB_Semiring_new(
553 : &second_second_edge, second_edge_monoid, second_bool_edge
554 : )) ;
555 :
556 : //--------------------------------------------------------------------------
557 : // Make E Vector
558 : //--------------------------------------------------------------------------
559 :
560 5027 : GRB_TRY (GrB_Vector_new(&E_vec, Ai_type, 2 * e)) ;
561 :
562 : // Filling out E_vec helps assign be much quicker.
563 5017 : GRB_TRY (GrB_assign(
564 : E_vec, NULL, NULL, 0, GrB_ALL, 0, NULL));
565 5012 : GrB_Index stride[] = {(GrB_Index) 0, e * 2 - 1, (GrB_Index) 2} ;
566 :
567 : // Shuffle i and j into E_vec.
568 5012 : GRB_TRY (GrB_Vector_assign(
569 : E_vec, NULL, NULL, Aj, stride, GxB_STRIDE, NULL)) ;
570 4997 : stride[GxB_BEGIN] = 1;
571 4997 : GRB_TRY (GrB_Vector_assign(
572 : E_vec, NULL, NULL, Ai, stride, GxB_STRIDE, NULL)) ;
573 :
574 : // Reinterpret E_vec as a vector of tuples.
575 4987 : GRB_TRY (GxB_Vector_unload(
576 : E_vec, &indices, &E_type, &e, &ind_size, &E_hand, NULL));
577 4987 : e /= 2;
578 4987 : GRB_TRY (GxB_Vector_load(
579 : E_vec, &indices, lg_edge, e, ind_size, E_hand, NULL));
580 :
581 : // Find Hash Size
582 4987 : int shift_e = 63 - (int) floor (log2 ((double) e)) ;
583 4987 : uint64_t ehash_size = (1ull << (67-shift_e)) ;
584 : // printf("Hash Size: %ld\n", ehash_size);
585 :
586 : //--------------------------------------------------------------------------
587 : // Initialize the rest of the vectors
588 : //--------------------------------------------------------------------------
589 :
590 : // Hash values and buckets
591 4987 : GRB_TRY (GrB_Vector_new(&exists, GrB_INT8, ehash_size)) ;
592 4977 : GRB_TRY (GrB_Vector_new(&hashed_edges, GrB_UINT64, e)) ;
593 4967 : GRB_TRY(GrB_set (exists, GxB_BITMAP | GxB_FULL, GxB_SPARSITY_CONTROL)) ;
594 :
595 : // Ramp
596 4957 : GRB_TRY (GrB_Vector_new(&ramp_v, GrB_UINT64, e + 1)) ;
597 4947 : GRB_TRY (GrB_Vector_assign_UINT64 (ramp_v, NULL, NULL, 0, GrB_ALL, 0, NULL)) ;
598 4942 : GRB_TRY (GrB_Vector_apply_IndexOp_UINT64 (ramp_v, NULL, NULL,
599 : GrB_ROWINDEX_INT64, ramp_v, 0, NULL)) ;
600 : GrB_Index ramp_size;
601 :
602 : // Constants
603 4932 : GRB_TRY (GrB_Scalar_new (&one8, GrB_UINT8)) ;
604 4922 : GRB_TRY (GrB_Scalar_setElement_UINT8 (one8, 1)) ;
605 4887 : GRB_TRY (GrB_Vector_new (&x, GrB_BOOL, e)) ;
606 :
607 : // Random Vector
608 4877 : GRB_TRY (GrB_Vector_new(&random_v, GrB_UINT64, e)) ;
609 4867 : GRB_TRY (GrB_Vector_new(&r_60, GrB_UINT64, e)) ;
610 4857 : GRB_TRY (GrB_Vector_new(&r_permute, GrB_UINT64, 1ull << (64-shift_e))) ;
611 4847 : GRB_TRY(GrB_set (r_permute, GxB_BITMAP, GxB_SPARSITY_CONTROL)) ;
612 4837 : GRB_TRY (GrB_Vector_assign_UINT64 (
613 : random_v, NULL, NULL, 0, GrB_ALL, e, NULL)) ;
614 4832 : LG_TRY(
615 : LAGraph_Random_Seed(random_v, seed, msg)) ;
616 :
617 : // printf("Entering loop, Good Luck:\n") ;
618 23346 : while(num_swaps < totSwaps)
619 : {
620 : GrB_Index perm_size, arr_size, junk_size;
621 : // r_60 has the random vector shifted by some amount.
622 27880 : GRB_TRY (GrB_Vector_apply_BinaryOp2nd_UINT64(
623 : r_60, NULL, NULL, GxB_BSHIFT_UINT64, random_v, -(shift_e), NULL
624 : )) ;
625 23187 : GRB_TRY (GrB_Vector_clear(x)) ;
626 23159 : GRB_TRY (GrB_Vector_resize(x, e)) ;
627 23113 : GRB_TRY (GrB_Vector_assign_BOOL(
628 : x, NULL, NULL, true, GrB_ALL, 0, NULL)) ;
629 23085 : LG_TRY (LAGraph_FastAssign_Semiring(
630 : r_permute, NULL, NULL, r_60, x, ramp_v, GxB_ANY_FIRSTJ_INT64,
631 : NULL, msg)) ;
632 :
633 22441 : GrB_Index edges_permed = 0;
634 22441 : GRB_TRY (GrB_Vector_nvals(&edges_permed, r_permute)) ;
635 22441 : GRB_TRY (GrB_Vector_new(&edge_perm, GrB_BOOL, edges_permed)) ;
636 22385 : GRB_TRY (GrB_Vector_extractTuples(NULL, edge_perm, r_permute, NULL)) ;
637 22329 : n_keep = LAGRAPH_MIN((int)(e * loopTry * 0.5), edges_permed / 2) ;
638 :
639 : // Chose only the edges we need from edge_perm
640 22329 : GRB_TRY (GrB_Vector_resize(edge_perm, n_keep * 2)) ;
641 :
642 : // Get the desired edges from the E_vec array.
643 22329 : GRB_TRY (GrB_Vector_new(&M, lg_edge, n_keep * 2)) ;
644 22273 : GRB_TRY (GxB_Vector_extract_Vector(
645 : M, NULL, NULL, E_vec, edge_perm, NULL
646 : )) ;
647 :
648 : // Make the swaps via the swap_pair unary op.
649 22217 : GRB_TRY (GxB_Vector_unload(
650 : M, (void **) &indices, &M_type, &M_nvals, &ind_size, &M_hand,
651 : NULL
652 : )) ;
653 22217 : GRB_TRY (GxB_Vector_load(
654 : M, (void **) &indices, lg_swap, M_nvals / 2, ind_size, M_hand, NULL
655 : )) ;
656 22217 : GRB_TRY (GrB_Vector_apply_IndexOp_BOOL(
657 : M, NULL, NULL, swap_pair, M, false, NULL)) ;
658 22217 : GRB_TRY (GxB_Vector_unload(
659 : M, (void **) &indices, &M_type, &M_nvals, &ind_size, &M_hand,
660 : NULL
661 : )) ;
662 22217 : GRB_TRY (GxB_Vector_load(
663 : M, (void **) &indices, lg_edge, M_nvals * 2, ind_size, M_hand,
664 : NULL
665 : )) ;
666 :
667 : // Hash Edges ----------------------------------------------------------
668 22217 : GRB_TRY (GrB_Vector_new(
669 : &new_hashed_edges, GrB_UINT64, n_keep * 2)) ;
670 :
671 22161 : GRB_TRY (GrB_Vector_apply_BinaryOp2nd_UINT64(
672 : new_hashed_edges, NULL, NULL, hash_seed_e, M,
673 : ehash_size - 1ll, NULL
674 : )) ;
675 22133 : GRB_TRY (GrB_Vector_apply_BinaryOp2nd_UINT64(
676 : hashed_edges, NULL, NULL, hash_seed_e, E_vec,
677 : ehash_size - 1ll, NULL
678 : )) ;
679 :
680 : //----------------------------------------------------------------------
681 : // Build Hash Buckets
682 : //----------------------------------------------------------------------
683 :
684 22105 : GRB_TRY (GrB_Vector_new(&dup_swaps_v, GrB_INT8, n_keep * 2)) ;
685 22049 : GRB_TRY (GrB_set(dup_swaps_v, GxB_BITMAP, GxB_SPARSITY_CONTROL)) ;
686 :
687 21993 : GRB_TRY (GrB_Vector_clear(x)) ;
688 21965 : GRB_TRY (GrB_Vector_resize(x, e)) ;
689 21965 : GRB_TRY (GrB_Vector_assign_BOOL(
690 : x, NULL, NULL, true, GrB_ALL, 0, NULL)) ;
691 : // place a one in any bucket that coresponds to an edge currently in E
692 21937 : LG_TRY (LAGraph_FastAssign_Semiring(
693 : exists, NULL, NULL, hashed_edges, x, ramp_v, GxB_ANY_PAIR_UINT8,
694 : NULL, msg
695 : )) ;
696 :
697 21293 : GRB_TRY (GrB_Vector_clear(x)) ;
698 21265 : GRB_TRY (GrB_Vector_resize(x, n_keep * 2)) ;
699 21013 : GRB_TRY (GrB_Vector_assign_BOOL(
700 : x, NULL, NULL, true, GrB_ALL, 0, NULL)) ;
701 :
702 : // exists cannot possibly be full at this point.
703 : // Fill out exists in O(1) time.
704 20985 : GRB_TRY (GxB_Container_new(&con)) ;
705 20817 : GRB_TRY (GxB_unload_Vector_into_Container(exists, con, NULL)) ;
706 : // Sanity check
707 20817 : LG_ASSERT (con->format == GxB_BITMAP, GrB_INVALID_VALUE) ;
708 20817 : GrB_free(&exists);
709 20817 : exists = con->b;
710 20817 : con->b = NULL;
711 20817 : GRB_TRY (GrB_free(&con)) ;
712 : // exist has to be full at this point
713 :
714 :
715 : // "Count" all of the edges that fit into each bucket. Stop counting at
716 : // 2 since we will have to throw that whole bucket away anyway.
717 20817 : LG_TRY (LAGraph_FastAssign_Semiring(
718 : exists, NULL, add_term_biop, new_hashed_edges, x, ramp_v,
719 : plus_term_one, NULL, msg
720 : )) ;
721 20201 : GRB_TRY(GrB_set (exists, GxB_BITMAP | GxB_FULL, GxB_SPARSITY_CONTROL)) ;
722 : // Select buckets with only one corresponding value
723 20201 : GRB_TRY (GrB_Vector_select_INT8(
724 : exists, NULL, NULL, GrB_VALUEEQ_UINT8, exists, (int8_t) 1,
725 : NULL
726 : )) ;
727 :
728 : // Find each hashed edge's bucket, dup_swaps_v is 1 if exists[edge] = 1
729 20089 : LG_TRY (LAGraph_FastAssign_Semiring(
730 : dup_swaps_v, NULL, NULL, new_hashed_edges, exists, ramp_v,
731 : GxB_ANY_PAIR_INT8, GrB_DESC_T0, msg
732 : )) ;
733 : // GRB_TRY (GxB_Vector_extract_Vector(
734 : // dup_swaps_v, NULL, NULL, exists, new_hashed_edges, NULL)) ;
735 :
736 : // Fill out dup_swaps_v in O(1) time.
737 19529 : GRB_TRY (GxB_Container_new(&con)) ;
738 19361 : GRB_TRY (GxB_unload_Vector_into_Container(dup_swaps_v, con, NULL)) ;
739 19361 : n_keep = con->nvals;
740 19361 : GRB_TRY (GrB_free(&dup_swaps_v));
741 19361 : dup_swaps_v = con->b;
742 19361 : con->b = NULL;
743 19361 : GRB_TRY (GrB_free(&con)) ;
744 19361 : GRB_TRY (GxB_Vector_unload(
745 : dup_swaps_v, (void **) &dup_swaps, &M_type, &M_nvals, &dup_arr_size,
746 : &M_hand, NULL)) ;
747 19361 : GRB_TRY (GxB_Vector_load(
748 : dup_swaps_v, (void **) &dup_swaps, GrB_INT16, M_nvals / 2,
749 : dup_arr_size, M_hand, NULL)) ;
750 19361 : GRB_TRY (GrB_apply(
751 : dup_swaps_v, NULL, NULL, lg_shiftland, dup_swaps_v, NULL)) ;
752 19305 : GRB_TRY (GxB_Vector_unload(
753 : dup_swaps_v, (void **) &dup_swaps, &M_type, &M_nvals, &dup_arr_size,
754 : &M_hand, NULL)) ;
755 19305 : GRB_TRY (GxB_Vector_load(
756 : dup_swaps_v, (void **) &dup_swaps, GrB_INT8, M_nvals * 2,
757 : dup_arr_size, M_hand, NULL)) ;
758 :
759 19305 : GRB_TRY (GrB_Vector_clear(exists)) ;
760 : // ---------------------------------------------------------------------
761 : // Place Good Swaps back into E_vec
762 : // ---------------------------------------------------------------------
763 :
764 19305 : GRB_TRY (GxB_Container_new(&con)) ;
765 19137 : GRB_TRY (GxB_unload_Vector_into_Container(M, con, NULL)) ;
766 19137 : GRB_TRY (GrB_free(&(con->b))) ;
767 19137 : con->b = dup_swaps_v;
768 : // n_keep = sum (dup_swaps_v), to count the # of 1's that now appear in
769 : // dup_swaps_v, which will become nvals(M) after loading M from the
770 : // container con.
771 19137 : GRB_TRY (GrB_reduce (&n_keep, NULL, GrB_PLUS_MONOID_UINT64,
772 : dup_swaps_v, NULL)) ;
773 19137 : con->nvals = n_keep;
774 19137 : con->format = GxB_BITMAP;
775 19137 : dup_swaps_v = NULL;
776 19137 : GRB_TRY (GxB_load_Vector_from_Container(M, con, NULL)) ;
777 : // GRB_TRY (GxB_print (M, 1)) ; // for debugging; must be OK here
778 19137 : GRB_TRY (GrB_free(&con)) ;
779 19137 : GRB_TRY (LAGraph_FastAssign_Semiring(
780 : E_vec, NULL, second_edge, edge_perm, M, ramp_v,
781 : second_second_edge, NULL, msg)) ;
782 :
783 18521 : n_keep /= 2;
784 :
785 18521 : FREE_LOOP ; // Free Matricies that have to be rebuilt
786 :
787 18521 : num_swaps += n_keep ;
788 18521 : LG_TRY (LAGraph_Random_Next(random_v, msg)) ;
789 : // printf("Made %ld swaps this loop. "
790 : // "[%.3f%% of Planned, %.3f%% of edges swapped]\n"
791 : // "Completed %ld swaps total. [%.3f%% of Planned]\n",
792 : // n_keep, n_keep * 100.0 / totSwaps, n_keep * 200.0 / e, num_swaps,
793 : // num_swaps * 100.0 / totSwaps) ;
794 18521 : if(n_keep < (int) (loopMin * e / 2) + 1)
795 : {
796 : // printf("Too Few Swaps occured! Exiting.\n");
797 1 : break;
798 : }
799 : }
800 133 : GRB_TRY (GxB_Vector_unload(
801 : E_vec, (void **) &indices, &lg_edge, &e, &ind_size, &E_hand, NULL));
802 133 : GRB_TRY (GxB_Vector_load(
803 : E_vec, (void **) &indices, E_type, e * 2, ind_size, E_hand, NULL));
804 133 : GRB_TRY (GrB_Vector_extract(
805 : Aj, NULL, NULL, E_vec, stride, GxB_STRIDE, NULL));
806 123 : stride[GxB_BEGIN] = 0;
807 123 : GRB_TRY (GrB_Vector_extract(
808 : Ai, NULL, NULL, E_vec, stride, GxB_STRIDE, NULL));
809 : // Build Output Matrix
810 113 : GRB_TRY (GrB_Matrix_new(&A_new, GrB_BOOL, n, n)) ;
811 98 : GRB_TRY (GxB_Matrix_build_Scalar_Vector(A_new, Ai, Aj, one8, NULL)) ;
812 63 : GRB_TRY (GrB_eWiseAdd(
813 : A_new, NULL, NULL, GrB_LOR_MONOID_BOOL, A_new,A_new, GrB_DESC_T0
814 : )) ;
815 23 : LAGRAPH_TRY (LAGraph_New (
816 : G_new, &A_new, LAGraph_ADJACENCY_UNDIRECTED, msg)) ;
817 18 : LG_FREE_WORK ;
818 18 : (*pSwaps) = num_swaps ;
819 18 : return (num_swaps >= totSwaps)? GrB_SUCCESS : LAGRAPH_INSUFFICIENT_SWAPS ;
820 : #else
821 : // printf("LAGr_SwapEdges Needs GB v10\n") ;
822 : return (GrB_NOT_IMPLEMENTED) ;
823 : #endif
824 : }
|
