1:- module(orth, []).    2
    3:- use_module(zdd('zdd-array')).    4:- use_module(zdd(zdd)).    5:- use_module(pac(op)).    6
    7		/*******************
    8		*	 helpers       *
    9		*******************/
   10
   11% ?- setup_ctrl_opt([rect(2,2)], X), write(X).
   12% ?- setup_ctrl_opt([rect(2,2), gc(all), a-b], X), write(X).
   13setup_ctrl_opt(X, ctrl(GC, S, T, Rect)):-
   14	memberchk(rect(W, H), X),
   15	Rect = rect(W, H),
   16	(	memberchk(gc(GC), X) -> true
   17	;	GC = 0
   18	),
   19	(	memberchk(S-T, X) -> true
   20	;	S = p(0,0),
   21		T = p(W, H)
   22	).
   23
   24% get ctrl option.
   25get_ctrl(X, end, Y)	:-!, arg(3, X, Y).
   26get_ctrl(X, gc, Y)	:-!, arg(1, X, Y).
   27get_ctrl(X, start, Y):-!, arg(2, X, Y).
   28get_ctrl(X, rect, Y):- arg(4, X, Y).
   29
   30% ?- arrow_symbol(_->_, F).
   31% ?- arrow_symbol(a->b, F, X, Y).
   32arrow_symbol( _ -> _).
   33%
   34arrow_symbol(A, A0):- functor(A, A0, 2).
   35arrow_symbol(A, A0, A1, A2):- functor(A, A0, 2),
   36		arg(1, A, A1),
   37		arg(2, A, A2).
   38
   39% One of the most basic helpers.
   40composable_pairs(A-B, A-C, B, C).
   41composable_pairs(A-B, C-A, B, C).
   42composable_pairs(B-A, A-C, B, C).
   43composable_pairs(B-A, C-A, B, C).
   44%
   45normal_pair(A-B, U-V):-!, ( B @< A -> U=B, V=A; U=A, V=B ).
   46normal_pair(A->B, U->V):- ( B @< A -> U=B, V=A; U=A, V=B ).
   47
   48% ?- frontier_node(p(0, 0)-[p(0,1), p(1,0)], X).
   49% ?- frontier_node(p(2, 2)-[], X).
   50% ?- frontier_node(p(1, 1)-[p(2,1)], X).
   51% ?- frontier_node(p(1, 1)-[p(1,2)], X).
   52
   53% For rectangular grid graph
   54frontier_node(_-[_, P], P):-!.
   55frontier_node(p(_, J)-[p(K, J)], p(K, J)):-!.
   56frontier_node(p(I, J)-[p(I, _)], p(K, J)):-!, K is I+1.
   57frontier_node(p(I, J)-[], p(K, J)):- K is I + 1.
   58
   59% ?- rect_nodes(rect(0,2), Ns).
   60% ?- rect_nodes(rect(10,10), Ns), length(Ns, L).
   61rect_nodes(rect(W, H), Ns):-
   62	findall(p(I,J),
   63			 (	between(0, W, I),
   64				between(0, H, J)
   65			 ),
   66			 Ns).
   67
   68% ?- rect_udg(rect(1,1), UDG).
   69% ?- rect_udg(rect(10,10), UDG),length(UDG, L).
   70rect_udg(rect(W, H), UDG):-
   71	findall( p(I,J)-Suc,
   72				 (	between(0, W, I),
   73					between(0, H, J),
   74					findall( p(K,L),
   75						(  L=J, K is I + 1, K =< W
   76						;  K=I, L is J + 1, L =< H
   77						),
   78						Suc0),
   79					sort(Suc0, Suc)
   80				 ),
   81				 UDG).
   82
   83
   84		/************************
   85		*     sample queries    *
   86		************************/
   87% ?- zdd.
   88% ?- rect_mate(rect(0,0), X), card(X, C).
   89% ?- rect_mate(rect(1,0), X), card(X, C).
   90% ?- rect_mate(rect(1,1), X), card(X, C).
   91% ?- rect_mate(rect(1,2), X), card(X, C).
   92% ?- rect_mate(rect(1,3), X), card(X, C).
   93% ?- rect_mate(rect(1,3), X), card(X, C).
   94% ?- rect_mate(rect(1,4), X), card(X, C).
   95% ?- rect_mate(rect(2,0), X), card(X, C).
   96% ?- rect_mate(rect(2,1), X), card(X, C).
   97% ?- rect_mate(rect(2,1), X), show_array.
   98% ?- rect_mate(rect(3,3), X), card(X, C).
   99% ?- rect_mate(rect(4,4), X), card(X, C).
  100% ?- time((rect_mate(rect(5,5), X), card(X, C))).
  101%@ % 18,080,809 inferences, 0.850 CPU
  102%@ X = 10102,
  103%@ C = 1262816 .
  104% ?- time((rect_mate(rect(6,6), [gc(0)], X), card(X, C))).
  105%@ % 91,423,943 inferences, 5.054 CPU
  106%@ X = 40678,
  107%@ C = 575780564 .
  108
  109% ?- time((rect_mate(rect(7,7), [gc(0)], X), card(X, C))).
  110%@ % 459,450,802 inferences, 29.595 CPU
  111%@ X = 115391,
  112%@ C = 789360053252 .

@ % 2,177,476,416 inferences, 148.247 CPU @ X = 548826, @ C = 3266598486981642 .

  119% ?- call_with_time_limit(2000, time((rect_mate(rect(9,9), X), card(X, C)))).
  120%@ % 9,066,262,242 inferences, 811.674 CPU
  121%@ X = 2328421,
  122%@ C = 41044208702632496804.
  123
  124% ?- call_with_time_limit(36000, time((rect_mate(rect(10,10), X), card(X, C)))).
  125%@ % 29,396,165,381 inferences, 3825.651 CPU in 3841.534 seconds (100% CPU, 7683964 Lips)
  126%@ X = 1268341,
  127%@ S = ..,
  128%@ C = 1568758030464750013214100.
  129
  130% [2022/10/01]
  131% ?- call_with_time_limit(36000, time((rect_mate(rect(11,11), X), card(X, C)))).
  132%@ % 129,787,391,715 inferences, 14036.125 CPU in 14092.162 seconds (100% CPU, 9246668 Lips)
  133%@ X = 4207535,
  134%@ S = ..,
  135%@ C = 182413291514248049241470885236.
  136
  137% ?- call_with_time_limit(300, time((rect_mate(rect(12,12), X, [gc(all)]), card(X, C)))).
  138
  139% ?- X is 3600*24*2.
  140%@ X = 172800.
  141
  142%	24 hours. (86400/3600).
  143% ?- call_with_time_limit(172800, time((rect_mate(rect(12,12), X, [gc(all)]), card(X, C)))).
  144
  145% ?- call_with_time_limit(360000, time((rect_mate(rect(13,13), X, [gc(0)], S), card(X, C, S)))).
  146
  147		/**************
  148		*     main    *
  149		**************/
  150%
  151rect_mate(Rect, X):-	udg_mate([gc(0), Rect], X).
  152%
  153rect_mate(Rect, Opts, X):- udg_mate([Rect|Opts], X).
  154%
  155udg_mate(Opts, X):-
  156	setup_ctrl_opt(Opts, Ctrl),
  157	get_ctrl(Ctrl, rect, Rect),
  158	rect_udg(Rect, Dg0),
  159	reverse(Dg0, Dg),
  160	udg_mate(Ctrl, Dg, 1, X0),
  161	get_ctrl(Ctrl, start, Start),
  162	get_ctrl(Ctrl, end, End),
  163	prune_final(Start, End, X0, X).
  164%
  165udg_mate(_, [], X, X).
  166udg_mate(Ctrl, [N-Ps|UDG], X, Y):-
  167	get_ctrl(Ctrl, end, End),
  168	frontier_node(N-Ps, F),
  169	add_node(F-End, N, Ps, X, X0),
  170	prune_by_classify_link(F, End, X0, X1),
  171	ctrl_gc(Ctrl, N, X1, X2),
  172	udg_mate(Ctrl, UDG, X2, Y).
  173
  174		/****************************************************
  175		*     replace_end/bypass:  Most basic operations.   *
  176		****************************************************/
  177
  178add_node(_, N, [], X, Y):-!, zdd_insert(N-N, X, Y).
  179add_node(FE, p(0,0), [P], X, Y):-!, N=p(0,0),
  180	subst_node(FE, [N->P], P, N, X, Y).
  181add_node(FE, N, [P], X, Y):-!,
  182	subst_node(FE, [N->P], P, N, X, X0),
  183	zdd_insert(N-N, X, X1),
  184	zdd_join(X0, X1, Y).
  185add_node(FE, p(0,0), [P, Q], X, Y):-!, N=p(0,0),
  186	subst_node(FE, [N->P], P, N, X, X0),
  187	subst_node(FE, [N->Q], Q, N, X, X1),
  188	zdd_join(X0, X1, Y).
  189add_node(FE, N, [P, Q], X, Y):-
  190	subst_node(FE, [N->P], P, N, X, X0),
  191	subst_node(FE, [N->Q], Q, N, X, X1),
  192 	bypass(FE, [N->P, N->Q], P-Q, X, X2),
  193	zdd_insert(N-N, X, X3),
  194	zdd_join_list([X0, X1, X2, X3], Y).
  195
  196% ?- spy(subst_node).
  197% ?- zdd((X<< +[[c-d]], @(subst_node(c-a, [a->c], c, a, X, Y)), psa(X), psa(Y))).
  198% ?- #(S), X<< +[[c-d]], subst_node(a-x, [u->v], c, b, X, Y, S), psa(X, S), psa(Y, S).
  199subst_node(_, _, _, _, X, 0):- X < 2, !.
  200subst_node(FE, Es, A, P, X, Y):-	 FE = Fr-End,		% replace A with P
  201	cofact(X, t(U, L, R)),
  202 	subst_node(FE, Es, A, P, L, L0),
  203	classify_link(on_frontier_for_subst, Fr, End, U, Case),
  204	arrow_symbol(U, _, Lu, Ru),
  205	(	( Case = 0 ; Case = arrow ) -> R0 = 0
  206	;	Case = ignore ->
  207		subst_node(FE, Es, A, P, R, R0)
  208	;	(	Ru = A	->
  209			normal_pair(Lu-P, V),
  210			zdd_ord_insert([V|Es], R, R0)
  211		;	Lu = A	->
  212			normal_pair(P-Ru, V),
  213			zdd_ord_insert([V|Es], R, R0)
  214		;	subst_node(FE, Es, A, P, R, R1),
  215			zdd_insert(U, R1, R0)
  216		)
  217	),
  218	zdd_join(L0, R0, Y).
  219%
  220bypass(_, _, _, X, 0):- X<2, !.
  221bypass(FE, Bypass, U, X, Y):- FE = Fr-E,
  222	cofact(X, t(A, L, R)),
  223	bypass(FE, Bypass, U, L, L0),
  224	classify_link(on_frontier, Fr, E, A, Case),
  225	(	( Case = 0; Case = arrow ) ->	R0 = 0			% many hits.
  226	;	Case = ignore -> bypass(FE, Bypass, U, R, R0)		% no hits.
  227	;	(	A = U -> R0 = 0		%,  write(.)			% so so hits
  228		; 	composable_pairs(U, A, U0, V0) ->
  229			subst_node(FE, Bypass, U0, V0, R, R0)
  230		;	bypass(FE, Bypass, U, R, R1),
  231			zdd_insert(A, R1, R0)
  232		)
  233	),
  234	zdd_join(L0, R0, Y).
  235
  236% NOT USED HERE. But as a prototype.
  237add_link(_, _, X, 0):- X<2, !.
  238add_link(FE, U, X, Y):- FE = F-E,
  239	cofact(X, t(A, L, R)),
  240	add_link(FE, U, L, L0),
  241	classify_link(on_frontier, F, E, A, Case),
  242	(	( Case = 0; Case = arrow ) ->	R0 = 0			% many hits.
  243	;	Case = ignore -> add_link(FE, U, R, R0)		% no hits.
  244	;	U = Ul-Ur,
  245		(	A = U -> R0 = 0		%,  write(.)			% so so hits
  246		; 	composable_pairs(U, A, U0, V0) ->
  247			subst_node(FE, [Ul->Ur], U0, V0, R, R0)
  248		;	add_link(FE, U, R, R1),
  249			zdd_insert(A, R1, R0)
  250		)
  251	),
  252	zdd_join(L0, R0, Y).
  253
  254		/***************
  255		*     prune    *
  256		***************/
  257
  258% ?- zdd X<< {[a-b, a->b]}, prune_final(a, b, X, Y).
  259prune_final(_, _, X, 0):- X<2, !.
  260prune_final(P, Q, X, Y):- cofact(X, t(A, L, R)),
  261	prune_final(P, Q, L, L0),
  262	(	A = (_->_) -> R0 = 0
  263	;  	A = P-Q -> prune_final0(R, R0)
  264	;	A = U-U -> prune_final(R, R0)
  265	;	R0 = 0
  266	),
  267	zdd_join(L0, R0, Y).
  268%
  269prune_final0(X, X):- X<2, !.
  270prune_final0(X, Y):- cofact(X, t(A, L, R)),
  271	prune_final0(L, L0),
  272	(	A = (_->_) -> zdd_insert(A, R, R0)
  273  	;	A = (B-B) -> prune_final0(R, R0)
  274	;	R0 = 0
  275	),
  276	zdd_join(L0, R0, Y).
  277%
  278ctrl_gc(Ctrl, P, X, Y):- get_ctrl(Ctrl, gc, GC),
  279	ctrl_gc_case(GC, P, X, Y).
  280%
  281ctrl_gc_case(all, P, X, Y):-!,
  282	writeln(all),
  283	format("GC at ~w \n", [P]),
  284	slim_gc(X, Y).		% better
  285%	zdd_gc(X, Y).
  286
  287ctrl_gc_case(K, P, X, Y):- integer(K), !,
  288	P = p(_, J),
  289	(	J = K ->
  290		format("GC at ~w \n", [P]),
  291%		zdd_gc(X, Y)
  292		slim_gc(X, Y)	% better
  293	;	Y=X
  294	).
  295ctrl_gc_case(_, _, X, X).
  296
  297		/***********************************
  298		*     classify_link by frontier    *
  299		***********************************/
  300
  301%  "A node P is on frontier" means that P may be touched by a new link in the future.
  302
  303% ?- on_frontier(p(0,3), p(1,1)). % false.
  304% ?- on_frontier(p(1,0), p(1,1)).
  305
  306%  Left most columnn aware frontier.
  307on_frontier(P, p(1, K)):-!,   % for pruning after adding nodes.
  308	(	P = p(1,J) -> J < K
  309	;	P = p(0,J),
  310		J =< K+1
  311	).
  312% Also works.
  313% on_frontier(P, p(1, K)):-!,   % for pruning after adding nodes.
  314% 	(	P = p(1,J) -> J =< K
  315% 	;	P = p(0,J),
  316% 		J =< K + 1
  317% 	).
  318on_frontier(P, F):- P @< F.
  319%
  320on_frontier_for_subst(P, p(1, K)):-!,   % for during substituting nodes.
  321	(	P=p(1,J) -> J =< K
  322	;	P=p(0,J), J =< K+1
  323	).
  324on_frontier_for_subst(P, F):- P @=< F.
  325
  326:- meta_predicate classify_link(2,  ?, ?, ?, ?).  327classify_link(_, _, _, _->_, arrow):-!.
  328classify_link(Pred, F, End, A-B,  Case):- call(Pred, A, F), !,
  329   	(	call(Pred, B, F) -> Case = keep
  330	;	B = End -> Case = keep
  331	;	Case = 0
  332	).
  333classify_link(_, _, E, E-E, 0):-!.		% added. Many hits, though only at beginning
  334classify_link(_, _, _, A-A, ignore):-!.
  335classify_link(_, _, _, _, 0).
  336
  337% ?- classify_link(p(1,1), p(3,3), p(0,2)-p(3,3), X).  % keep
  338% ?- classify_link(p(0,0), p(3,3), p(1,1)-p(1,1), X).  % ignore
  339% ?- classify_link(p(0,0), p(3,3), p(1,0)-p(1,1), X).  % 0
  340%
  341classify_link_for_bypass(F, E, A, Case):- classify_link(on_frontier, F, E, A, Case).
  342%
  343classify_link_for_subst(X, Y, Z, U):- classify_link(on_frontier_for_subst, X, Y, Z, U).
  344
  345%
  346prune_by_classify_link(_, _, X, X):- X<2, !.
  347prune_by_classify_link(F, End, X, Y):- cofact(X, t(A, L, R)),
  348	prune_by_classify_link(F, End, L, L0),
  349	classify_link(on_frontier, F, End, A, Case),    %  better than classify_link0 as for End-End
  350	(	Case = arrow -> zdd_insert(A, R, R0)
  351	;	Case = keep ->					% many hits.
  352		prune_by_classify_link(F, End, R, R1),
  353		zdd_insert(A, R1, R0)
  354	;	Case = ignore ->				% many hits.
  355		prune_by_classify_link(F, End, R, R0)
  356	;	R0 = 0							% many bits.
  357	),
  358	zdd_join(L0, R0, Y)