1# 2# V 1.1 - Bug fixed 3# 4# By Steve Linton 5# 6# Bad documentation by Chris Jefferson 7# 8# Finds the minimal image of a Set set under a group G. 9# 10# Usage: NewSmallestImage(G, set, stab, x -> x); 11# 12# (ignore the last argument!) 13# 14# Where stab should be a subgroup of 15# Stabilizer(G,set); 16# 17# If in doubt, the best way to invoke this algorithm 18# is: 19# NewSmallestImage(G, set, Stabilizer(G, set, OnSets), x -> x); 20# 21# Returns a pair [image, stabilizer], where stabilizer is a subgroup of Stabilizer(G, set), possibly larger than the one given into the function. 22# 23# Note: The return type of this is NOT a set, but provides the pointwise mapping of the input set. 24# This means the permutation which actually provides the smallest image can be found cheaply as follows: 25# 26# res := NewSmallestImage(G, set, Group(()), x -> x); 27# perm := RepresentativeAction(G, set, res[1], OnTuples); 28 29 30 31# 32# Search node data: 33# 34# selected -- indices in set of points being mapped to minimal image at this node 35# image -- sequence-wise image of set under element represented by this node 36# substab -- Stab_K(selected) sequence stabilizer 37# children -- nodes corresponding to extensions of selected 38# parent -- node corr to all but last element of selected.] 39# childno 40# next -- across row 41# prev -- likewise 42# deleted 43# 44 45# At each level 46 47# Find the next pt of minimum image and all the corresponding nodes 48 49# If at any time in this process we notice a potential prune, we have a 50# generator of K -- construct it, close K with it and delete all tree 51# branches that are now non-canonical -- propagate new K down through 52# tree. Continue with surviving nodes at current level. 53 54_IMAGES_NSI_HASH_LIMIT :=100; 55 56 57_IMAGES_TIME_CLASSES := []; 58 59_IMAGES_DeclareTimeClass := function(name) 60 BindGlobal(name, Length(_IMAGES_TIME_CLASSES)+1); 61 Add(_IMAGES_TIME_CLASSES,MakeImmutable(name)); 62end; 63 64 65_IMAGES_DeclareTimeClass("pass1"); 66_IMAGES_DeclareTimeClass("pass2"); 67_IMAGES_DeclareTimeClass("pass3"); 68_IMAGES_DeclareTimeClass("shortcut"); 69_IMAGES_DeclareTimeClass("changeStabChain"); 70_IMAGES_DeclareTimeClass("orbit"); 71_IMAGES_DeclareTimeClass("skippedorbit"); 72_IMAGES_DeclareTimeClass("getcands"); 73_IMAGES_DeclareTimeClass("improve"); 74_IMAGES_DeclareTimeClass("check1"); 75_IMAGES_DeclareTimeClass("check2"); 76_IMAGES_DeclareTimeClass("prune"); 77_IMAGES_DeclareTimeClass("ShallowNode"); 78_IMAGES_DeclareTimeClass("DeepNode"); 79_IMAGES_DeclareTimeClass("FilterOrbCount"); 80 81_IMAGES_TIME_CLASSES := MakeImmutable(_IMAGES_TIME_CLASSES); 82 83_IMAGES_nsi_stats := ListWithIdenticalEntries(Length(_IMAGES_TIME_CLASSES),0); 84 85_IMAGES_DO_TIMING := true; 86if IsBound( MakeThreadLocal ) then 87 MakeThreadLocal("_IMAGES_DO_TIMING"); 88fi; 89 90if _IMAGES_DO_TIMING then 91 _IMAGES_StartTimer := function(cat) 92 _IMAGES_nsi_stats[cat] := _IMAGES_nsi_stats[cat] - Runtime(); 93 end; 94 95 _IMAGES_StopTimer := function(cat) 96 _IMAGES_nsi_stats[cat] := _IMAGES_nsi_stats[cat] + Runtime(); 97 end; 98 99 _IMAGES_IncCount := function(cat) 100 _IMAGES_nsi_stats[cat] := _IMAGES_nsi_stats[cat] + 1; 101 end; 102 103 _IMAGES_ResetStats := function() 104 _IMAGES_nsi_stats := ListWithIdenticalEntries(Length(_IMAGES_TIME_CLASSES),0); 105 end; 106 107 _IMAGES_ResetStats(); 108 if IsBound( MakeThreadLocal ) then 109 MakeThreadLocal("_IMAGES_nsi_stats"); 110 fi; 111 112 _IMAGES_GetStats := function() 113 local r, c; 114 r := rec(); 115 for c in _IMAGES_TIME_CLASSES do 116 r.(c) := _IMAGES_nsi_stats[ValueGlobal(c)]; 117 od; 118 return r; 119 end; 120 121else 122 _IMAGES_StartTimer := function(cat) 123 return; 124 end; 125 126 _IMAGES_StopTimer := function(cat) 127 return; 128 end; 129 130 _IMAGES_IncCount := function(cat) 131 return; 132 end; 133 134 _IMAGES_ResetStats := function() 135 return; 136 end; 137 138 _IMAGES_GetStats := function() 139 return fail; 140 end; 141 142fi; 143 144_IMAGES_Get_Hash := function(m) 145 local jenkins_hash; 146 if IsBoundGlobal("JENKINS_HASH") then 147 jenkins_hash := ValueGlobal("JENKINS_HASH"); 148 return s->jenkins_hash(s,GAPInfo.BytesPerVariable*m+GAPInfo.BytesPerVariable); 149 else 150 return s->HashKeyBag(s,57,0,GAPInfo.BytesPerVariable*m+GAPInfo.BytesPerVariable); 151 fi; 152end; 153 154 155# GAP dictionaries don't (currently) provide a way of getting the values 156# stored in them, so here we cache them seperately 157_countingDict := function(dictexample) 158 local data; 159 data := rec( 160 d := NewDictionary(dictexample, true), 161 l := [] 162 ); 163 164 return rec( 165 add := function(list) 166 local val; 167 val := LookupDictionary(data.d, list); 168 if val = fail then 169 val := 0; 170 Add(data.l, list); 171 fi; 172 val := val + 1; 173 AddDictionary(data.d, list, val); 174 end, 175 176 findElement := function(comp) 177 local smallval, smalllist, val, i; 178 smalllist := data.l[1]; 179 smallval := LookupDictionary(data.d, smalllist); 180 for i in data.l do 181 val := LookupDictionary(data.d, i); 182 if comp(val, smallval) or (val = smallval and i < smalllist) then 183 smallval := val; 184 smalllist := i; 185 fi; 186 od; 187 return smalllist; 188 end, 189 190 dump := function() return data; end 191 ); 192end; 193 194 195if not IsBound(InfoNSI) then 196 DeclareInfoClass("InfoNSI"); 197fi; 198 199_IMAGES_RATIO := function(selector) 200 return function(orbmins, orbitCounts, orbsizes) 201 local index, result, i, ret; 202 index := 1; 203 result := [selector(1, orbmins, orbitCounts, orbsizes), orbmins[1]]; 204 for i in [2..Length(orbmins)] do 205 ret := [selector(i, orbmins, orbitCounts, orbsizes), orbmins[i]]; 206 if (orbitCounts[index] = 0) or (ret < result and orbitCounts[i] <> 0) then 207 index := i; 208 result := ret; 209 fi; 210 od; 211 return index; 212 end; 213end; 214 215_IMAGES_RARE_RATIO_ORBIT := _IMAGES_RATIO( 216 function(i, orbmins, orbitCounts, orbsizes) 217 return (Log2(Float(orbitCounts[i])))/orbsizes[i]; 218 end 219); 220 221_IMAGES_COMMON_RATIO_ORBIT := _IMAGES_RATIO( 222 function(i, orbmins, orbitCounts, orbsizes) 223 return -(Log2(Float(orbitCounts[i])))/orbsizes[i]; 224 end 225); 226 227_IMAGES_RARE_RATIO_ORBIT_FIX := _IMAGES_RATIO( 228 function(i, orbmins, orbitCounts, orbsizes) 229 if(orbsizes[i]) = 1 then return Float(-(2^32)+orbitCounts[i]); fi; 230 return (Log2(Float(orbitCounts[i])))/orbsizes[i]; 231 end 232); 233 234_IMAGES_COMMON_RATIO_ORBIT_FIX := _IMAGES_RATIO( 235 function(i, orbmins, orbitCounts, orbsizes) 236 if(orbsizes[i]) = 1 then return Float(-(2^32)-orbitCounts[i]); fi; 237 return -(Log2(Float(orbitCounts[i])))/orbsizes[i]; 238 end 239); 240 241_IMAGES_RARE_ORBIT := _IMAGES_RATIO( 242 function(i, orbmins, orbitCounts, orbsizes) 243 return orbitCounts[i]; 244 end 245); 246 247_IMAGES_COMMON_ORBIT := _IMAGES_RATIO( 248 function(i, orbmins, orbitCounts, orbsizes) 249 return -orbitCounts[i]; 250 end 251); 252 253 254_NewSmallestImage := function(g,set,k,skip_func, early_exit, disableStabilizerCheck_in, config_option) 255 local leftmost_node, next_node, delete_node, delete_nodes, 256 clean_subtree, handle_new_stabilizer_element, 257 simpleOrbitReps, make_orbit, n, s, l, m, hash, 258 lastupb, root, depth, gens, orbnums, orbmins, 259 orbsizes, upb, imsets, imsetnodes, node, cands, y, 260 x, num, rep, node2, prevnode, nodect, changed, 261 newnode, image, dict, seen, he, bestim, bestnode, 262 imset, p, 263 config, 264 globalOrbitCounts, globalBestOrbit, minOrbitMset, orbitMset, 265 savedArgs, 266 countOrbDict, 267 bestOrbitMset 268 ; 269 270 if IsString(config_option) then 271 config_option := ValueGlobal(config_option); 272 fi; 273 274 if config_option.branch = "static" then 275 savedArgs := rec( config_option := config_option, g := g, k := k, set := set ); 276 if config_option.order = "MinOrbit" then 277 savedArgs.perm := MinOrbitPerm(g); 278 elif config_option.order = "MaxOrbit" then 279 savedArgs.perm := MaxOrbitPerm(g); 280 else 281 ErrorNoReturn("Invalid 'order' when branch = 'static' in CanonicalImage"); 282 fi; 283 savedArgs.perminv := savedArgs.perm^-1; 284 g := g^savedArgs.perm; 285 k := k^savedArgs.perm; 286 set := OnTuples(set, savedArgs.perm); 287 config_option := rec(branch := "minimum"); 288 else 289 savedArgs := rec(perminv := ()); 290 fi; 291 292 if config_option.branch = "minimum" then 293 config := rec( 294 skipNewOrbit := function() return (upb <= lastupb + 1); end, 295 getQuality := pt -> orbmins[pt], 296 getBasePoint := IdFunc, 297 initial_lastupb := 0, 298 initial_upb := infinity, 299 countRareOrbits := false, 300 tryImproveStabilizer := true, 301 preFilterByOrbMset := false, 302 ); 303 elif config_option.branch = "dynamic" then 304 config := rec(skipNewOrbit := ReturnFalse, 305 preFilterByOrbMset := false); 306 if config_option.order in ["MinOrbit", "MaxOrbit", "SingleMaxOrbit"] then 307 config.getBasePoint := pt->pt[2]; 308 config.initial_lastupb := [-infinity, -infinity]; 309 config.initial_upb := [infinity, infinity]; 310 config.countRareOrbits := false; 311 config.tryImproveStabilizer := true; 312 313 if config_option.order = "MinOrbit" then 314 config.getQuality := pt -> [orbsizes[pt], orbmins[pt]]; 315 elif config_option.order = "MaxOrbit" then 316 config.getQuality := pt -> [-orbsizes[pt], orbmins[pt]]; 317 elif config_option.order = "SingleMaxOrbit" then 318 config.getQuality := function(pt) 319 if orbsizes[pt] = 1 then 320 return [-(2^64), orbmins[pt]]; 321 else 322 return [-orbsizes[pt], orbmins[pt]]; 323 fi; 324 end; 325 else 326 ErrorNoReturn("?"); 327 fi; 328 elif config_option.order in ["RareOrbit", "CommonOrbit", "RareRatioOrbit", "CommonRatioOrbit", 329 "RareRatioOrbitFix", "CommonRatioOrbitFix"] then 330 config.getBasePoint := IdFunc; 331 config.initial_lastupb := 0; 332 config.initial_upb := infinity; 333 config.countRareOrbits := true; 334 config.tryImproveStabilizer := false; 335 config.getQuality := pt -> orbmins[pt]; 336 if config_option.order = "RareOrbit" then 337 config.calculateBestOrbit := _IMAGES_RARE_ORBIT; 338 elif config_option.order = "CommonOrbit" then 339 config.calculateBestOrbit := _IMAGES_COMMON_ORBIT; 340 elif config_option.order = "RareRatioOrbit" then 341 config.calculateBestOrbit := _IMAGES_RARE_RATIO_ORBIT; 342 elif config_option.order = "RareRatioOrbitFix" then 343 config.calculateBestOrbit := _IMAGES_RARE_RATIO_ORBIT_FIX; 344 elif config_option.order = "CommonRatioOrbit" then 345 config.calculateBestOrbit := _IMAGES_COMMON_RATIO_ORBIT; 346 elif config_option.order = "CommonRatioOrbitFix" then 347 config.calculateBestOrbit := _IMAGES_COMMON_RATIO_ORBIT_FIX; 348 else 349 ErrorNoReturn("?"); 350 fi; 351 else 352 ErrorNoReturn("Invalid ordering: ", config_option.order); 353 fi; 354 355 if IsBound(config_option.orbfilt) then 356 if config_option.orbfilt = "Min" then 357 # This space intensionally blank 358 elif config_option.orbfilt = "Rare" then 359 config.findBestOrbMset := function(x,y) return x < y; end; 360 elif config_option.orbfilt = "Common" then 361 config.findBestOrbMset := function(x,y) return x > y; end; 362 else 363 Error("Invalid 'orbfilt' option"); 364 fi; 365 config.preFilterByOrbMset := true; 366 fi; 367 else 368 ErrorNoReturn("'branch' must be minimum, static or dynamic"); 369 fi; 370 371 if disableStabilizerCheck_in = true then 372 config.tryImproveStabilizer := false; 373 fi; 374 375 ## Node exploration functions 376 leftmost_node := function(depth) 377 local n, i; 378 n := root; 379 while Length(n.selected) < depth -1 do 380 n := n.children[1]; 381 od; 382 return n; 383 end; 384 385 next_node := function(node) 386 local n; 387 n := node; 388 repeat 389 n := n.next; 390 until n = fail or not n.deleted; 391 return n; 392 end; 393 394 # Delete a node, and recursively deleting all it's children. 395 delete_node := function(node) 396 local i; 397 if node.deleted then 398 return; 399 fi; 400 Info(InfoNSI,3,"Deleting ",node.selected); 401 if node.prev <> fail then 402 node.prev.next := node.next; 403 fi; 404 if node.next <> fail then 405 node.next.prev := node.prev; 406 fi; 407 node.deleted := true; 408 if node.parent <> fail then 409 Remove(node.parent.children, node.childno); 410 if Length(node.parent.children) = 0 then 411 delete_node(node.parent); 412 else 413 for i in [node.childno..Length(node.parent.children)] do 414 node.parent.children[i].childno := i; 415 od; 416 fi; 417 fi; 418 if IsBound(node.children) then 419 delete_nodes(ShallowCopy(node.children)); 420 fi; 421 end; 422 delete_nodes := function(nodes) 423 local node; 424 for node in nodes do 425 delete_node(node); 426 od; 427 end; 428 429 # Filter nodes by stabilizer group, 430 # Updates the stabilizer group of the node, 431 clean_subtree := function(node) 432 local bad, seen, c, x, q, gens, olen, pt, gen, im; 433 Info(InfoNSI,3,"Cleaning at ",node.selected); 434 if not IsBound(node.children) then 435 return; 436 fi; 437 bad := []; 438 439 seen := BlistList([1..m],[]); 440 for c in node.children do 441 if IsBound(c.selectedbaselength) then 442 x := c.selected[c.selectedbaselength]; 443 else 444 x := c.selected[Length(c.selected)]; 445 fi; 446 if seen[x] then 447 Info(InfoNSI, 5, "Removing ", c, " because ", x); 448 Add(bad,c); 449 else 450 q := [x]; 451 gens := GeneratorsOfGroup(node.substab); 452 olen := 1; 453 Info(InfoNSI, 5, "Keeping ", c, " because ", x); 454 seen[x] := true; 455 for pt in q do 456 for gen in gens do 457 im := pt^gen; 458 if not seen[im] then 459 seen[im] := true; 460 Add(q,im); 461 olen := olen+1; 462 fi; 463 od; 464 od; 465 if olen < Size(node.substab)/Size(c.substab) then 466 c.substab := Stabilizer(node.substab,x); 467 clean_subtree(c); 468 fi; 469 fi; 470 od; 471 delete_nodes(bad); 472 end; 473 474 # Add a new stabilizer element, mapping node1 to node2, and then call 475 # clean_subtree to remove any new subtrees. 476 handle_new_stabilizer_element := function(node1,node2) 477 local perm1, i; 478 # so node1 and node2 represnet group elements that map set to the same 479 # place in two different ways 480 perm1 := PermListList(node1.image, node2.image); 481 Info(InfoNSI, 2, "Can map ",node1.image, " to ", node2.image, " : ", perm1); 482 Assert(1, not perm1 in l); 483 l := ClosureGroup(l,perm1); 484 Info(InfoNSI,2,"Found new stabilizer element. Stab now ",Size(l)); 485 root.substab := l; 486 clean_subtree(root); 487 end; 488 489 # Given a group 'gp' and a set 'set', find orbit representatives 490 # of 'set' in 'gp' simply. 491 simpleOrbitReps := function(gp,set) 492 local m, n, b, seed, reps, gens, q, pt, gen, im; 493 m := Length(set); 494 n := set[m]; 495 b := BlistList([1..n],set); 496 seed := set[1]; 497 reps := []; 498 gens := GeneratorsOfGroup(gp); 499 while seed <> fail and seed <= n do 500 b[seed] := false; 501 q := [seed]; 502 Add(reps,seed); 503 for pt in q do 504 for gen in gens do 505 im := pt^gen; 506 if b[im] then 507 b[im] := false; 508 Add(q,im); 509 fi; 510 od; 511 od; 512 seed := Position(b,true,seed); 513 od; 514 return reps; 515 end; 516 517 # Make orbit of x, updating orbnums, orbmins and orbsizes as approriate. 518 make_orbit := function(x) 519 local q, rep, num, pt, gen, img; 520 if orbnums[x] <> -1 then 521 return orbnums[x]; 522 fi; 523 q := [x]; 524 rep := x; 525 num := Length(orbmins)+1; 526 orbnums[x] := num; 527 for pt in q do 528 for gen in gens do 529 img := pt^gen; 530 if orbnums[img] = -1 then 531 orbnums[img] := num; 532 Add(q,img); 533 if img < rep then 534 rep := img; 535 fi; 536 fi; 537 od; 538 od; 539 Add(orbmins,rep); 540 Add(orbsizes,Length(q)); 541 return num; 542 end; 543 544 if set = [] then 545 return [ [], k^(savedArgs.perminv)]; 546 fi; 547 548 n := Maximum(LargestMovedPoint(g), Maximum(set)); 549 s := StabChainMutable(g); 550 l := Action(k,set); 551 m := Length(set); 552 hash := _IMAGES_Get_Hash(m); 553 lastupb := config.initial_lastupb; 554 root := rec(selected := [], 555 image := set, 556 imset := Immutable(Set(set)), 557 substab := l, 558 deleted := false, 559 next := fail, 560 prev := fail, 561 parent := fail); 562 for depth in [1..m] do 563 Info(InfoNSI, 3, "Stabilizer is :", s.generators); 564 gens := s.generators; 565 orbnums := ListWithIdenticalEntries(n,-1); 566 orbmins := []; 567 orbsizes := []; 568 upb := config.initial_upb; 569 imsets := []; 570 imsetnodes := []; 571 # 572 # At this point, all bottom nodes are blue 573 # first pass creates appropriate set of virtual red nodes 574 # 575 _IMAGES_StartTimer(pass1); 576 577 if IsBound(config.findBestOrbMset) then 578 countOrbDict := _countingDict([1,2,3]); 579 node := leftmost_node(depth); 580 while node <> fail do 581 _IMAGES_StartTimer(getcands); 582 cands := Difference([1..m],skip_func(node.selected)); 583 584 _IMAGES_StopTimer(getcands); 585 orbitMset := []; 586 for y in cands do 587 _IMAGES_IncCount(check1); 588 x := node.image[y]; 589 num := make_orbit(x); 590 Add(orbitMset, orbmins[num]); 591 od; 592 Sort(orbitMset); 593 countOrbDict.add(orbitMset); 594 node := next_node(node); 595 od; 596 597 bestOrbitMset := countOrbDict.findElement(config.findBestOrbMset); 598 Unbind(countOrbDict); # Free memory 599 fi; 600 601 if config.preFilterByOrbMset then 602 minOrbitMset := [infinity]; 603 node := leftmost_node(depth); 604 while node <> fail do 605 _IMAGES_StartTimer(getcands); 606 cands := Difference([1..m],skip_func(node.selected)); 607 608 _IMAGES_StopTimer(getcands); 609 orbitMset := []; 610 for y in cands do 611 _IMAGES_IncCount(check1); 612 x := node.image[y]; 613 num := make_orbit(x); 614 Add(orbitMset, orbmins[num]); 615 od; 616 Sort(orbitMset); 617 Info(InfoNSI, 5, "Considering: ", orbitMset, "::",node.selected); 618 if IsBound(bestOrbitMset) then 619 if orbitMset <> bestOrbitMset then 620 delete_node(node); 621 fi; 622 else 623 if orbitMset < minOrbitMset then 624 Info(InfoNSI, 4, "New min: ", orbitMset); 625 minOrbitMset := orbitMset; 626 node2 := node.prev; 627 while node2 <> fail do 628 Info(InfoNSI, 4, "Clean up old big set"); 629 _IMAGES_IncCount(FilterOrbCount); 630 delete_node(node2); 631 node2 := node2.prev; 632 od; 633 elif orbitMset > minOrbitMset then 634 Info(InfoNSI, 4, "Too big!"); 635 delete_node(node); 636 fi; 637 fi; 638 639 node := next_node(node); 640 od; 641 fi; 642 643 if config.countRareOrbits then 644 globalOrbitCounts := ListWithIdenticalEntries(Length(orbmins), 0) ; 645 node := leftmost_node(depth); 646 while node <> fail do 647 _IMAGES_StartTimer(getcands); 648 cands := Difference([1..m],skip_func(node.selected)); 649 if Length(cands) > 1 and not IsTrivial(node.substab) then 650 cands := simpleOrbitReps(node.substab,cands); 651 fi; 652 # 653 # These index the children of node that will 654 # not be immediately deleted under rule C 655 # 656 _IMAGES_StopTimer(getcands); 657 for y in cands do 658 _IMAGES_IncCount(check1); 659 x := node.image[y]; 660 num := make_orbit(x); 661 662 if IsBound(globalOrbitCounts[num]) then 663 globalOrbitCounts[num] := globalOrbitCounts[num] + 1; 664 else 665 globalOrbitCounts[num] := 1; 666 fi; 667 od; 668 node := next_node(node); 669 od; 670 671 globalBestOrbit := config.calculateBestOrbit(orbmins, globalOrbitCounts, orbsizes); 672 upb := orbmins[globalBestOrbit]; 673 fi; 674 675 676 node := leftmost_node(depth); 677 while node <> fail do 678 679 _IMAGES_StartTimer(getcands); 680 cands := Difference([1..m],skip_func(node.selected)); 681 if Length(cands) > 1 and not IsTrivial(node.substab) then 682 cands := simpleOrbitReps(node.substab,cands); 683 fi; 684 # 685 # These index the children of node that will 686 # not be immediately deleted under rule C 687 # 688 _IMAGES_StopTimer(getcands); 689 node.validkids := []; 690 for y in cands do 691 _IMAGES_IncCount(check1); 692 x := node.image[y]; 693 694 num := orbnums[x]; 695 if num = -1 then 696 # 697 # Need a new orbit. Also require the smallest point 698 # as the rep. 699 # 700 # 701 # If there is no prospect of the new orbit being 702 # better than the current best then go on to the next candidate 703 # 704 if config.skipNewOrbit() then 705 _IMAGES_IncCount(skippedorbit); 706 continue; 707 fi; 708 _IMAGES_StartTimer(orbit); 709 num := make_orbit(x); 710 _IMAGES_StopTimer(orbit); 711 rep := config.getQuality(num); 712 if rep < upb then 713 _IMAGES_StartTimer(improve); 714 ### CAJ - Support bailing out early when a smaller 715 # set is found 716 if early_exit[1] and rep < early_exit[2][depth] then 717 return [MinImage.Smaller, l^(savedArgs.perminv)]; 718 fi; 719 ### END of bailing out early 720 upb := rep; 721 node2 := node.prev; 722 while node2 <> fail do 723 delete_node(node2); 724 node2 := node2.prev; 725 od; 726 _IMAGES_IncCount(ShallowNode); 727 node.validkids := [y]; 728 Info(InfoNSI,3,"Best down to ",upb); 729 _IMAGES_StopTimer(improve); 730 fi; 731 else 732 _IMAGES_IncCount(check2); 733 rep := config.getQuality(num); 734 if rep = upb then 735 _IMAGES_IncCount(ShallowNode); 736 Add(node.validkids,y); 737 fi; 738 fi; 739 od; 740 if node.validkids = [] then 741 _IMAGES_StartTimer(prune); 742 delete_node(node); 743 _IMAGES_StopTimer(prune); 744 fi; 745 node := next_node(node); 746 od; 747 ### CAJ - Support bailing out early when a larger set is found 748 if early_exit[1] and upb > early_exit[2][depth] then 749 return [MinImage.Larger, l^(savedArgs.perminv)]; 750 fi; 751 ### 752 Info(InfoNSI,2,"Layer ",depth," pass 1 complete. Best is ",upb); 753 _IMAGES_StopTimer(pass1); 754 # 755 # Second pass. Actually make all the red nodes and turn them blue 756 # 757 lastupb := upb; 758 Info(InfoNSI, 2, "Branch on ", upb); 759 _IMAGES_StartTimer(changeStabChain); 760 ChangeStabChain(s,[config.getBasePoint(upb)],false); 761 _IMAGES_StopTimer(changeStabChain); 762 if Length(s.orbit) = 1 then 763 # 764 # In this case nothing much can happen. Each surviving node will have exactly one child 765 # and none of the imsets will change 766 # so we mutate the nodes in-place 767 # 768 _IMAGES_StartTimer(shortcut); 769 node := leftmost_node(depth); 770 while node <> fail do 771 if not IsBound(node.selectedbaselength) then 772 node.selectedbaselength := Length(node.selected); 773 fi; 774 Assert(1, Length(node.validkids)=1); 775 Add(node.selected, node.validkids[1]); 776 node := next_node(node); 777 od; 778 Info(InfoNSI,2,"Nothing can happen, short-cutting"); 779 s := s.stabilizer; 780 _IMAGES_StopTimer(shortcut); 781 if Size(skip_func(leftmost_node(depth+1).selected)) = m then 782 Info(InfoNSI,2,"Skip would skip all remaining points"); 783 break; 784 fi; 785 786 continue; # to the next depth 787 fi; 788 _IMAGES_StartTimer(pass2); 789 node := leftmost_node(depth); 790 prevnode := fail; 791 nodect := 0; 792 changed := false; 793 while node <> fail do 794 node.children := []; 795 for x in node.validkids do 796 _IMAGES_IncCount(DeepNode); 797 newnode := rec( selected := Concatenation(node.selected,[x]), 798 substab := Stabilizer(node.substab,x), 799 parent := node, 800 childno := Length(node.children)+1, 801 next := fail, 802 prev := prevnode, 803 deleted := false); 804 nodect := nodect+1; 805 if prevnode <> fail then 806 prevnode.next := newnode; 807 fi; 808 prevnode := newnode; 809 Add(node.children,newnode); 810 811 image := node.image; 812 if image[x] <> config.getBasePoint(upb) then 813 repeat 814 image := OnTuples(image, s.transversal[image[x]]); 815 until image[x] = config.getBasePoint(upb); 816 newnode.image := image; 817 newnode.imset := Set(image); 818 MakeImmutable(newnode.imset); 819 changed := true; 820 else 821 newnode.image := image; 822 newnode.imset := node.imset; 823 fi; 824# Print("Made a node ",newnode.selected, " ",newnode.image,"\n"); 825 od; 826 node := next_node(node); 827 od; 828 _IMAGES_StopTimer(pass2); 829 Info(InfoNSI,2,"Layer ",depth," pass 2 complete. ",nodect," new nodes"); 830 # 831 # Third pass detect stabilizer elements 832 # 833 834 _IMAGES_StartTimer(pass3); 835 if changed and config.tryImproveStabilizer then 836 node := leftmost_node(depth+1); 837 if nodect > _IMAGES_NSI_HASH_LIMIT then 838 dict := SparseHashTable(hash); 839 seen := []; 840 while node <> fail do 841 he := GetHashEntry(dict,node.imset); 842 if fail <> he then 843 handle_new_stabilizer_element(node, he); 844 else 845 AddHashEntry(dict, node.imset, node); 846# if hash(node.imset) in seen then 847# Error(""); 848# fi; 849# AddSet(seen, hash(node.imset)); 850 fi; 851 node := next_node(node); 852 od; 853 Info(InfoNSI,2,"Layer ",depth," pass 3 complete. Used hash table"); 854 s := s.stabilizer; 855 if Length(s.generators) = 0 then 856 Info(InfoNSI,2,"Run out of group, return best image"); 857 node := leftmost_node(depth+1); 858 bestim := node.imset; 859 bestnode := node; 860 node := next_node(node); 861 while node <> fail do 862 if node.imset < bestim then 863 bestim := node.imset; 864 bestnode := node; 865 fi; 866 node := next_node(node); 867 od; 868 _IMAGES_StopTimer(pass3); 869 return [OnTuples(bestnode.image,savedArgs.perminv),l^savedArgs.perminv]; 870 fi; 871 else 872 while node <> fail do 873 imset := node.imset; 874 p := PositionSorted(imsets, imset); 875 if p <= Length(imsets) and imsets[p] = imset then 876 handle_new_stabilizer_element(node, imsetnodes[p]); 877 else 878 Add(imsets,imset,p); 879 Add(imsetnodes,node,p); 880 fi; 881 node := next_node(node); 882 od; 883 Info(InfoNSI,2,"Layer ",depth," pass 3 complete. ",Length(imsets)," images"); 884 s := s.stabilizer; 885 if Length(s.generators) = 0 then 886 Info(InfoNSI,2,"Run out of group, return best image"); 887 _IMAGES_StopTimer(pass3); 888 return [OnTuples(imsetnodes[1].image,savedArgs.perminv),l^savedArgs.perminv]; 889 fi; 890 fi; 891 else 892 s := s.stabilizer; 893 fi; 894 _IMAGES_StopTimer(pass3); 895 if Size(skip_func(leftmost_node(depth+1).selected)) = m then 896 Info(InfoNSI,2,"Skip would skip all remaining points"); 897 break; 898 fi; 899 od; 900 return [OnTuples(leftmost_node(depth+1).image,savedArgs.perminv),l^savedArgs.perminv]; 901end; 902 903 904 905 906 907