1:- module(count_loops, []).    2
    3% (start at copy of minato-r4)
    4		/*****************************************************
    5		*     count cycles in the rectangular grid graphs    *
    6		*****************************************************/
    7:- use_module(util(math)).    8:- use_module(util(meta2)).    9:- use_module(pac(basic)).   10:- use_module(pac(meta)).   11:- use_module(util(misc)).   12:- use_module(pac('expand-pac')).   13:- use_module(zdd('zdd-array')).   14:- use_module(zdd(zdd)).   15:- use_module(zdd('zdd-misc')).   16:- use_module(zdd('zdd-plain')).   17:- use_module(pac(op)).   18:- op(1060, xfy, #).		% exor
   19:- op(1000, xfy, &).   20
   21% term_expansion --> pac:expand_pac.
   22term_expansion --> expand_pac.
   23
   24		/******************
   25		*     Switches    *
   26		******************/
   27% ORTHOGONAl
   28% active_node(X, Y):- active_node_orth(X, Y).
   29% rect_nodes(X, Y):- rect_nodes_orth(X, Y).
   30
   31% SKEW.
   32active_node(X, Y):- active_node_skew(X, Y).
   33rect_nodes(X, Y):- rect_nodes_skew(X, Y).
   34
   35% some tiny.
   36w(Opt):- memberchk(out(X), Opt), memberchk(state(S), Opt), psa(X, S).
   37
   38% ?- count_rect_cycles(rect(0,0), C, _).
   39% ?- count_rect_cycles(rect(1,0), C, _).
   40% ?- count_rect_cycles(rect(0,1), C, _).
   41% ?- count_rect_cycles(rect(1,1), C, U).
   42% ?- count_rect_cycles(rect(1,2), C, _).
   43% ?- count_rect_cycles(rect(2,1), C, _).
   44% ?- count_rect_cycles(rect(2,2), C, _).
   45% ?- count_rect_cycles(rect(1,10), C, _).
   46% ?- count_rect_cycles(rect(10,1), C, _).
   47% ?- count_rect_cycles(rect(100,1), C, _).
   48% ?- count_rect_cycles(rect(1,100), C, _).
   49% ?- repeat(between(1, 6, N), (check(N, C), writeln( count_for(N) =  C))).
   50check(N, C):- findall((I,J), (between(0, N, I), between(I, N, J), I < J),  S),
   51	 length(S, C).
   52
   53% ?- repeat(between(1, 6, N), (count_rect_cycles(rect(N, 1), C), writeln(count_rect_for(N) =  C))).
   54
   55% ?- N = 3, findall((I,J), (between(0, N, I), between(I, N, J), I < J),  S),
   56%	length(S, C).
   57
   58% ?- time(count_rect_cycles(rect(6,6), C, _)).
   59%@ % 25,642,240 inferences, 3.558 CPU in 4.465 seconds (80% CPU, 7207373 Lips)
   60%@ C = 575780564.
   61% ?- time(count_rect_cycles(rect(7,7), C, _)).
   62%@ % 106,041,498 inferences, 12.348 CPU in 12.454 seconds (99% CPU, 8587850 Lips)
   63%@ C = 789360053252.
   64% ?- time(count_rect_cycles(rect(8,8), C, _)).
   65%% [2022/08/10] change of rehash.  faster, but correct ?
   66%@ % 332,129,654 inferences, 33.514 CPU in 33.728 seconds (99% CPU, 9910274 Lips)
   67%@ C = 2318527339461265.
   68%@ % 436,367,556 inferences, 52.415 CPU in 52.981 seconds (99% CPU, 8325235 Lips)
   69%@ C = 3266598486981642.
   70% ?- time(count_rect_cycles(rect(9,9), C, _)).
   71%@ % 1,467,957,762 inferences, 167.674 CPU in 168.988 seconds (99% CPU, 8754851 Lips)
   72%@ C = 27359264067916806101.
   73%@ % 1,793,021,796 inferences, 224.894 CPU in 227.462 seconds (99% CPU, 7972753 Lips)
   74%@ C = 41044208702632496804.
   75% ?- call_with_time_limit(1800, time(count_rect_cycles(rect(10,10), C, S))).
   76%@ % 12,782,653,306 inferences, 1164.030 CPU in 1174.499 seconds (99% CPU, 10981375 Lips)
   77%@ C = 1568758030464750013214100,
   78% ?- call_with_time_limit(36000, count_rect_cycles(rect(11,11), C, _)).
   79%@ C = 182413291514248049241470885236.
   80% ?- convert_skew(p(1,1), Q), convert_skew(P, Q), convert_skew(P, R).
   81%@ Q = R, R = p(2, 1),
   82%@ P = p(1, 1).
   83orth_skew(p(I, J), Y):- var(Y), !,
   84	Y = p(K, I),
   85	K is I + J.
   86orth_skew(p(I, J), p(K, I)):- J is K - I.
   87
   88% ?- arrow_symbol(_->_, F).
   89% ?- arrow_symbol(a->b, F, X, Y).
   90change_symbol(A-B, A->B).
   91arrow_symbol( _ -> _).
   92arrow_symbol(A, A0):- functor(A, A0, 2).
   93arrow_symbol(A, A0, A1, A2):- functor(A, A0, 2),
   94		arg(1, A, A1),
   95		arg(2, A, A2).
   96%
   97is_link(->(_, _)).
   98%
   99convert_normal_pair(P, Q, S):- orth_skew_pair(P, Q0),
  100	normal_pair(Q0, Q, S).
  101%
  102normal_arrow(X, Y):- @(normal_arrow(X, Y)).
  103%
  104normal_arrow(A->B, Arr, S):-!,
  105	(	zdd_compare(>, A, B, S) -> Arr = (B->A)
  106	;	Arr = (A->B)
  107	).
  108normal_arrow(A-B, Arr, S):-normal_arrow(A->B, Arr, S).
  109%
  110
  111lt_pair(A - _, B -_, S):- zdd_compare(<, A, B, S).
  112
  113
  114		/*************************************
  115		*     basics of incrementing DG      *
  116		*************************************/
  117
  118composable_pairs(A-B, A-C, B, C).
  119composable_pairs(A-B, C-A, B, C).
  120composable_pairs(B-A, A-C, B, C).
  121composable_pairs(B-A, C-A, B, C).
  122%
  123normal_pair(P, U):- @(normal_pair(P, U)).
  124%
  125normal_pair(A-B, B-A, S):- zdd_compare(>, A, B, S), !.
  126normal_pair(A->B, B->A, S):- zdd_compare(>, A, B, S), !.
  127normal_pair(P, P, _).
  128
  129% ?- orth_skew_pair(p(1,1)-p(2,2), X), orth_skew_pair(Y, X).
  130%@ X = p(2, 1)-p(4, 2),
  131%@ Y = p(1, 1)-p(2, 2).
  132orth_skew_pair(P, Q):-
  133	(	P = (A-B) -> Q = (A0-B0)
  134	;	P = (A->B),
  135		Q = (A0->B0)
  136	),
  137	orth_skew(A, A0),
  138	orth_skew(B, B0).
  139
  140% ?- link_sort_opp([loop, a, loop], X).
  141link_sort_opp(X, Y):- sort(X, X0), reverse(X0, Y).
  142
  143% ?- node_ends(p(0,0), En, rect(1,1)).
  144% ?- node_ends(p(0,0), En, rect(0,0)).
  145% ?- node_ends(p(1,1), En, rect(2,2)).
  146
  147node_ends(P, En, rect(W,H)):- P = p(I,J),
  148	(	I < W ->
  149		K is I + 1,
  150		En0 = [p(K,J)]
  151	;	En0 = []
  152	),
  153	(	J < H ->
  154		L is J + 1,
  155		En = [p(I, L)|En0]
  156	;   En = En0
  157	).
  158
  159
  160% ?- node_ends_skew(p(0,0), En, rect(1,1)).
  161% ?- node_ends_skew(p(0,0), En, rect(0,0)).
  162% ?- node_ends_skew(p(1,1), En, rect(2,2)).
  163% ?- node_ends_skew(p(1,1), En, rect(1,2)).
  164node_ends_skew(P, En, rect(W,H)):- P = p(I,J),
  165	(	J < W ->
  166		K is J + 1,
  167		M is I + 1,
  168		En0 = [p(M, K)]
  169	;	En0 = []
  170	),
  171	(	I < H + J ->
  172		L is I + 1,
  173		En = [p(L, J)|En0]
  174	;   En = En0
  175	).
  176
  177% ?- R=rect(2,2), repeat(in_rect(P, R),
  178%	(node_ends_shadows(P, En, Sh, R), writeln(P-En-Sh))).
  179in_rect(p(I, J), rect(W, H)):- between(0, W, I), between(0, H, J).
  180
  181% ?- rect_nodes_orth(rect(1,1), PES), maplist(writeln, PES).
  182% ?- rect_nodes_orth(rect(3,3), PES), maplist(writeln, PES).
  183% ?- rect_nodes_orth(rect(10,10), PES), maplist(writeln, PES).
  184% ?- rect_nodes_orth(rect(30,30), PES), length(PES, L).
  185% ?- rect_nodes_orth(rect(100,100), PES), length(PES, L).
  186% ?- time((rect_nodes_orth(rect(1000,1000), PES), length(PES, L))).
  187
  188rect_nodes_orth(R, PES):-
  189	findall(P-E, (in_rect(P, R), node_ends(P, E, R)), PES).
  190
  191% ?- rect_nodes_skew(rect(100,100), PES), length(PES, L).
  192% ?- time((rect_nodes_skew(rect(1000,1000), PES), length(PES, L))).
  193rect_nodes_skew(R, PES_skew):- rect_nodes_orth(R, PES_orth),
  194	convert_orth_skew(PES_orth, PES_skew).
  195%
  196convert_orth_skew([], []).
  197convert_orth_skew([A-B|C], [A0-B0|C0]):-
  198	orth_skew(A, A0),
  199	maplist(orth_skew, B, B0),
  200	convert_orth_skew(C, C0).
  201
  202		/****************************************************
  203		*     replace_end/bypass:  Most basic operations.   *
  204		****************************************************/
  205
  206select_cycles(P, [A, B], X, Cyc, S):-!,
  207	select_cycles_([P->A, P->B], A-B, X, Cyc, S).
  208select_cycles(_, _, _, 0, _).
  209
  210% ?- zdd((X<<{[a-b, c-c, a->d, b->d], [a-b, c-c, c->f, b->f], [u-v]}, psa(X), @(select_cycles_([1->a, 1->b], a-b, X, Y)), psa(Y))).
  211select_cycles_(_, _, X, 0, _):- X<2, !.
  212select_cycles_(B, A, X, Y, S):- cofact(X, t(U, L, R), S),
  213	select_cycles_(B, A, L, L0, S),
  214 	select_cycles_(B, A, U, R, R0, S),
  215	zdd_join(L0, R0, Y, S).
  216
  217%
  218select_cycles_(B, A, A, X, Y, S):-!,
  219	collect_link_free(X, Y0, S),
  220	zdd_ord_insert(B, Y0, Y, S).
  221select_cycles_(B, A, U-U, X, Y, S):-!, select_cycles_(B, A, X, Y, S).
  222select_cycles_(_, _, _, _, 0, _).
  223
  224
  225% ?- zdd X<< {[a-b, c->d], [e->f]}, collect_link_free(X, Y), psa(Y).
  226% ?- zdd X<<{[a-b, c->d], [a-a, e->f]}, collect_link_free(X, Y), psa(Y).
  227collect_link_free(X, X, _S):- X<2, !.
  228collect_link_free(X, Y, S):- cofact(X, t(A, L, R), S),
  229	collect_link_free(L, L0, S),
  230	(	A = (U-U) -> collect_link_free(R, R0, S)
  231	;	A = (_-_) -> R0 = 0
  232	;	zdd_insert(A, R, R0, S)
  233	),
  234	zdd_join(L0, R0, Y, S).
  235
  236		/*************************
  237		*     count_rect_path    *
  238		*************************/
  239
  240% ?- test(rect(1,1), Y).
  241test(X, Y):- count_rect_cycles(X, Y, _, _).
  242
  243
  244% ?- N = 3,  count_rect_cycles(rect(0, N), C).
  245% ?- N = 3, count_rect_cycles(rect(1, N), C).
  246% ?- N = 10, count_rect_cycles(rect(1, N), C).
  247%% Qcompile: /Users/cantor/devel/zdd/prolog/zdd/count-loops.pl
  248%@ true.
  249%@ Context module: count_loops.
  250%@ true.
  251% ?- spy(node_one_by_one).
  252% ?- spy(psa).
  253% ?- count_rect_cycles(rect(0,0), C, S), w(S).
  254% ?- count_rect_cycles(rect(1,0), C, S), w(S).
  255% ?- count_rect_cycles(rect(0,1), C, S), w(S).
  256% ?- count_rect_cycles(rect(1,1), C, S).
  257% ?- count_rect_cycles(rect(2,1), C, S).
  258% ?- count_rect_cycles(rect(1,2), C, S).
  259% ?- count_rect_cycles(rect(3,1), C, S).
  260% ?- count_rect_cycles(rect(2,2), C, S).
  261% ?- count_rect_cycles(rect(2,1), C, S).
  262% ?- count_rect_cycles(rect(3,1), C, S).
  263% ?- count_rect_cycles(rect(3,3), C, S).
  264% ?- count_rect_cycles(rect(4,4), C, S).
  265
  266% [2022/08/08]
  267% ?- time(count_rect_cycles(rect(10,10), C, S)).
  268%@ % 6,626,002,510 inferences, 613.462 CPU in 618.651 seconds (99% CPU, 10801003 Lips)
  269%@ C = 988808811046283595068099,
  270%@ S = [state(..), out(2473748)|_].
  271% [2022/08/07]
  272% ?- time(count_rect_cycles(rect(10,10), C, S)).
  273%@ % 6,626,002,510 inferences, 639.629 CPU in 644.225 seconds (99% CPU, 10359131 Lips)
  274%@ C = 988808811046283595068099,
  275%@ S = [state(..), out(2473748)|_].
  276
  277
  278% ?- zdd((X<<pow([1,2]), card(X, C))).
  279% ?- time(count_rect_cycles(rect(8,8), C, S)).
  280%@ % 436,367,556 inferences, 35.323 CPU in 36.801 seconds (96% CPU, 12353625 Lip
  281%@ C = 3266598486981642,
  282%@ S = [state(..), out(4921865)|_].
  283% ?- call_with_time_limit(200, time(count_rect_cycles(rect(8,8), C, _))).
  284%@ % 908,479,618 inferences, 89.804 CPU in 91.150 seconds (99% CPU, 10116254 Lips)
  285%@ C = 3266598486981642.
  286% ?- count_rect_cycles(rect(10,10), C, _).
  287%@ C = 988808811046283595068099.
  288
  289count_rect_cycles(Rect, C):- count_rect_cycles(Rect, C, _).
  290%
  291count_rect_cycles(Rect, C, Opts):-
  292	memberchk(state(S), Opts),
  293	memberchk(out(X), Opts),
  294	open_state(S),
  295	set_key(rect, Rect, S),
  296	rect_nodes(Rect, A),
  297	sort(A, A1),
  298	reverse(A1, B),
  299	B =[Max-_|_],
  300	arg(1, Max, K),
  301	set_key(prune_index, K, S),
  302	node_one_by_one(B, 1-_, 0-X, S),
  303%	psa(X, S),
  304	card(X, C, S).
  305
  306% ?- zdd(@(node_one_by_one([], 1-X, 0-Y))).
  307node_one_by_one([], X-X, Y-Y, _).
  308node_one_by_one([P-En|R], X-X0, Y-Y0, S):-
  309	writeln(P),
  310	select_cycles(P, En, X, Cyc, S),
  311	zdd_join(Cyc, Y, Y1, S),
  312	add_node_one(P, En, X, X1, S),
  313	prune_inactive_node(P, X1, X2, S),
  314	get_key(prune_index, I, S),
  315	arg(1, P, J),
  316	(	J = I ->  X3 = X2,  Y2 = Y1
  317	;	set_key(prune_index, J, S),
  318		writeln("Slimming..."),
  319		zdd_slim([X2, Y1], [X3, Y2], S),
  320		garbage_collect
  321	),
  322	node_one_by_one(R, X3-X0, Y2-Y0, S).
  323
  324
  325% ?- zdd((X<<{[b-c, b->c]}, psa(X),  @(add_node_one(a, [b], X, Y)), psa(Y))).
  326add_node_one(P, [], X, Y, S):-!, zdd_insert(P-P, X, Y, S).
  327add_node_one(P, [A], X, Y, S):-!,
  328	subst_node([P->A], A, P, X, Y0, S),
  329	zdd_insert(P-P, X, Y1, S),
  330	zdd_join(Y0, Y1, Y, S).
  331add_node_one(P, [A, B], X, Y, S):-
  332	subst_node([P->A], A, P, X, Y0, S),
  333	subst_node([P->B], B, P, X, Y1, S),
  334	bypass_via_node([P->A, P->B], B-A, X, Y2, S),
  335	zdd_insert(P-P, X, Y3, S),
  336	zdd_join_list([Y0, Y1, Y2, Y3], Y, S).
  337
  338
  339% ?- zdd((X<<{[c-d]}, @(subst_node([a->c], c, a, X, Y)), psa(X), psa(Y))).
  340subst_node(_, _, _, X, 0, _):- X < 2, !.
  341subst_node(Es, A, P, X, Y, S):- % memo is useless here
  342	/* replace A with P */
  343	cofact(X, t(U, L, R), S),
  344	arrow_symbol(U, Arrow, Lu, Ru),
  345	(	Arrow = (->) ->	Y = 0
  346	;	zdd_compare(<, A, Lu, S) ->	Y = 0
  347	;	subst_node(Es, A, P, L, L0, S),
  348		(	Lu = A	->
  349			normal_pair(Ru-P, V, S),
  350			zdd_ord_insert([V|Es], R, R0, S)
  351		;	Ru = A	->
  352			normal_pair(Lu-P, V, S),
  353			zdd_ord_insert([V|Es], R, R0, S)
  354		;	subst_node(Es, A, P, R, R1, S),
  355			zdd_insert(U, R1, R0, S)
  356		),
  357		zdd_join(L0, R0, Y, S)
  358	).
  359
  360% ?- zdd((X<<{[a-b, c-d]}, @(bypass_via_node(b-c, X, Y)), psa(X), psa(Y))).
  361bypass_via_node(E, X, Y, S):- bypass_via_node([], E, X, Y, S).
  362
  363% ?- zdd((X<<{[a-b, c-d]}, @(bypass_via_node([b->c], b-c, X, Y)), psa(X), psa(Y))).
  364bypass_via_node(_, _, X, 0, _):- X < 2,!.
  365bypass_via_node(Bypass, E, X, Y, S):- % memo is useless here
  366	cofact(X, t(P, L, R),  S),
  367	(	P = (_->_)  ->  Y = 0
  368	;	lt_pair(E, P, S) -> Y = 0
  369	;  	bypass_via_node(Bypass, E, L, L0, S),
  370		(	E = P ->	R0 = 0
  371		;  	composable_pairs(E, P, U, V) ->
  372			subst_node(Bypass, U, V, R, R0, S)
  373		;	bypass_via_node(Bypass, E, R, R1, S),
  374			zdd_insert(P, R1, R0, S)
  375		),
  376		zdd_join(L0, R0, Y, S)
  377	).
  378
  379% ?- zdd((X<<{[p(1)-p(2)]}, @(prune_inactive_node(p(1), X, Y, p(2))))).
  380% ?- zdd((X<<{[p(1)-p(2)]}, @(prune_inactive_node(p(0), X, Y, p(3))))).
  381% ?- zdd((X<<{[p(1)-p(2), p(3)-p(4)]}, @(prune_inactive_node(p(2), X, Y, p(5))))).
  382% ?- zdd((X<<{[p(1)-p(5), p(3)-p(4)]}, @(prune_inactive_node(p(2), X, Y, p(5))))).
  383% ?- zdd((X<<{[p(1)-p(5)], [p(3)-p(4)]}, @(prune_inactive_node(p(2), X, Y, p(5))), psa(Y))).
  384% ?- zdd((X<<{[p(1)-p(5), p(2)-p(2)]}, @(prune_inactive_node(p(2), X, Y, p(5))), psa(Y))).
  385% ?- zdd((X<<{[p(1)-p(4), p(5)-p(5)]}, @(prune_inactive_node(p(2), X, Y, p(5))), psa(Y))).
  386% ?- zdd((X<<{[p(1)-p(4)], [p(5)-p(5)]}, @(prune_inactive_node(p(2), X, Y, p(5))), psa(Y))).
  387% ?- zdd((X<<{[p(1)-p(4)], [p(5)-p(5)]}, @(prune_inactive_node(p(5), X, Y, p(5))), psa(Y))).
  388
  389% ?- count_rect_cycles(rect(10, 1), C, S).
  390% ?- count_rect_cycles(rect(1,1), C, S).
  391% ?- R=rect(1,1), findall(on(P,Q), (in_rect(P, R), in_rect(Q, R), active_node(P, Q)), S),
  392%  maplist(writeln, S).
  393
  394% ?- R=rect(10,10), findall(on(P,Q), (in_rect(P, R), in_rect(Q, R),
  395%	active_node(P, Q)), S),
  396%   maplist(writeln, S), fail; true.
  397
  398% ?- active_node_orth(p(0,0),p(0,0)).  % true (conentionally)
  399% ?- active_node_orth(p(0,0),p(0,1)).
  400% ?- active_node_orth(p(0,0),p(1,0)).
  401% ?- active_node_orth(p(0,0),p(1,1)).
  402% ?- active_node_orth(p(0,1),p(1,0)).  % true.
  403% ?- active_node_orth(p(0,1),p(1,1)).
  404% ?- active_node_orth(p(1,0),p(1,1)).  % true
  405
  406% active_node_orth(P, Q) is true if Q is directly linked from some P' < P.
  407active_node_orth(p(0, J), p(0, L)):-!, L = J.
  408active_node_orth(p(I, J), p(I, L)):-!, J =< L.
  409active_node_orth(p(I, J), p(K, L)):- K is I+1, !, L < J.
  410
  411% ?- rect_nodes(rect(10,10), Ns), findall(Q, ( member(P-_, Ns),
  412%	orth_skew(P, Q)), Sk), sort(Sk, SSk),
  413%	findall(active(U, V), (member(U, SSk), member(V, SSk), active_node_skew(U, V)), W), maplist(writeln, W), fail; true.
  414
  415% ?- active_node_skew(p(0,0), p(0, 0)).  % true
  416% ?- active_node_skew(p(0,0), p(0, 1)).  % false (vacuously)
  417% ?- active_node_skew(p(1,1), p(2, 1)).  % true
  418% ?- active_node_skew(p(1,0), p(1, 1)).  % true
  419% ?- active_node_skew(p(1,1), p(2, 2)).  % false
  420% ?- active_node_skew(p(1,1), p(2, 1)).  % true.
  421% ?- active_node_skew(p(1,0), p(2, 1)).  % false
  422% ?- active_node_skew(p(1,1), p(2, 0)).  % true
  423active_node_skew(p(I, J), p(I, L)):- !, J =< L. % Actually not ncessary clause ?
  424active_node_skew(p(I, J), p(K, L)):- K is I + 1,
  425			(	 L = J ->  J > 0
  426			;	 L < J
  427			).
  428%
  429prune_inactive_node(_, X, X, _):- X<2, !.
  430prune_inactive_node(P, X, Y, S):-
  431	cofact(X, t(A, L, R), S),
  432	(	A = (_->_) -> Y = X
  433	; 	prune_inactive_node(P, L, L0, S),
  434		A = (U - V),
  435		(	active_node(P, U) ->
  436			(	active_node(P, V) ->
  437				prune_inactive_node(P, R, R1, S),
  438				zdd_insert(A, R1, R0, S)
  439			;	R0 = 0
  440			)
  441		;	U = V ->
  442			prune_inactive_node(P, R, R0, S)
  443		;	R0 = 0
  444		),
  445		zdd_join(L0, R0, Y, S)
  446	)