1:- module(mwp, []).    2
    3		/*******************************************************************
    4		*     Experimental for mate with its possible paths s.             *
    5		*     mp(i-j, s). However this is slow unexpectedly. 10 minutes
    6	    *	  for rect(4,4).
    7		*     subst_node is critical which uses memo with address of s,    *
    8		*     which could be large integers. Not sure whether term_hash    *
    9		*     is effective for such large integer.                         *
   10		*******************************************************************/
   11
   12
   13:- use_module(zdd('zdd-array')).   14:- use_module(zdd(zdd)).   15:- use_module(pac(op)).   16
   17udg_path(End, PathSet):- get_key(coa, Coa),
   18	set_key(ends, End),
   19	udg_mate_prune(Coa, 1, PathSet).
   20
   21
   22% ?- X<< (@a), Y<< (@c), connect_mate_path(X, b, Y, Z), psa(Z).
   23connect_mate_path(X , A, Y, Z):- zdd_insert(A, X, X0),
   24	zdd_merge(X0, Y, Z).
   25
   26		/******************************
   27		*     counting path in UDG    *
   28		******************************/
   29% ?- set_compare(compare).
   30
   31% ?- udg_path_count(a-d, [a-b, b-c, c-d], C).
   32% A big circular graph.
   33% ?- N = 10,
   34%	findall(A,	( between(1, N, I), I<N, J is I + 1, A=(I-J)), Ls),
   35%	time(udg_path_count(1-2,[1-N |Ls], C)).
   36%
   37udg_path_count(Ends, Links, C):-	udg_path_count(Ends, Links, C, _).
   38
   39%
   40udg_path_count(Ends, Links, C, X):-
   41		prepare_udg(Ends, Links),
   42		!,
   43		get_key(coa, Coa),
   44		udg_mate_prune(Coa, 1, X),
   45		card(X, C).
   46%
   47udg_mate_prune(Ls, X, Y):-
   48	add_links(Ls, X, Y0),
   49	get_key(ends, A-B),
   50	prune_final(A, B, Y0, Y).
   51
   52% ?- get_compare(X).
   53% ?- time((rect_path_count(p(3,3)-p(4,4), rect(6,6), C))).
   54% ?- time(rect_path_count(rect(1,1), C)).
   55% ?- time(rect_path_count(rect(3,3), C)).
   56% ?- time(rect_path_count(rect(4,4), C)).
   57% ?- time(rect_path_count(rect(6,6), C)).
   58% ?- time(rect_path_count(rect(7,7), C)).
   59
   60rect_path_count(R, C):- R = rect(W, H),
   61	rect_path_count( p(0,0) - p(W,H), R, C, _).
   62%
   63rect_path_count(Ends, R, C, X):- rect_links(R, Links),
   64	udg_path_count(Ends, Links, C, X).
   65
   66% ?- rect_links(rect(1,1), Links).
   67% ?- rect_links(rect(10,10), Links),length(Links, L).
   68rect_links(rect(W, H), Links):-
   69	findall( p(I,J) - p(K,L),
   70				 (	between(0, W, I),
   71					between(0, H, J),
   72					(  L = J, K is I + 1, K =< W
   73					;  K = I, L is J + 1, L =< H
   74					)
   75				 ),
   76				 Links).
   77
   78		/********************************
   79		*     Prepare UDG in coalgebra  *
   80		********************************/
   81
   82% ?- trace.
   83% ?- spy(add_links).
   84
   85% ?- udg_path_count(a-b,[a-b], C).
   86% ?- udg_path_count(a-c,[a-b, a-d, b-c, d-c], C).
   87% ?- udg_path_count(a-f,[a-b, a-d, b-c, d-c, d-e, e-f, c-f], C).
   88% ?- udg_path_count(a-k,
   89%	[a-b, a-d, b-c,
   90%	d-c, d-e, e-f,
   91%	c-f, b-g, g-h,
   92%	h-k, c-h, f-k], C).
   93
   94% ?- prepare_udg([a-b, b-c, c-d]), memo(node_id(a)-Id),
   95% get_key(dom, X).
   96
   97% ?- prepare_udg(a-d, [a-b, b-c, c-d]),
   98%	get_key(ends, ST),
   99%	get_key(frontier, F),
  100%	get_key(dom, Dom),
  101%	get_key(coa, U).
  102
  103prepare_udg(ST, Links):-
  104	prepare_udg(Links),
  105	prepare_ends(ST, A-B),
  106	set_key(ends, A-B).
  107%
  108prepare_udg(Links):-
  109	prepare_udg(Links, Succs, D, Vec),
  110	length(D, N),
  111	completing_succs(Succs, Succs0, N),
  112	set_key(coa, Succs0),
  113	set_key(dom, D),
  114	set_key(frontier, Vec).
  115%
  116prepare_udg(Links, Succs, D, Vec):-
  117	intern_links(Links, Links0),
  118	normal_mate_list(Links0, Links1),
  119	sort(Links1, Links2),
  120	domain_of_links(Links2, D),
  121	rel_to_fun(Links2, Succs),
  122	Vec=..[#|D],
  123	setup_frontier(Links1, Vec).
  124%
  125prepare_ends(A-B, R):-!, R = A0-B0,
  126	memo(node_id(A)-I),
  127	memo(node_id(B)-J),
  128	(	nonvar(I), nonvar(J) -> sort([I, J], [A0, B0])
  129	;	format("No link at ~w or ~w\n", [A,B]),
  130		fail
  131	).
  132prepare_ends(E, _):-
  133	format("Unexpected form of end nodes ~w \n", [E]),
  134	fail.
  135
  136% ?-completing_succs([], Y, 2).
  137% ?-completing_succs([2-[a]], Y, 3).
  138completing_succs(X, X, 0):-!.
  139completing_succs([I-A|Ls], [I-A|Ms], I):-!, J is I - 1,
  140	completing_succs(Ls, Ms, J).
  141completing_succs(Ls, [I-[]|Ms], I):- J is I - 1,
  142	completing_succs(Ls, Ms, J).
  143
  144% ?- normal_mate_list([1-2], X).
  145% ?- normal_mate_list([2-1, 1-2], X).
  146normal_mate_list([], []).
  147normal_mate_list([P|R], [P0|R0]):- P = I-J,
  148	(	J @< I -> P0 = J - I
  149	;	P0 = P
  150	),
  151	normal_mate_list(R, R0).
 rel_to_fun(+R, -F) is det
convert set R of links to a list F of successor lists. In other words: F is a function derived from the relation R such that F(x) = P (x in dom(R)) if P = { y | R(x,y)} e.g. R=[a-b, a-c, b-d, b-e] => F=[a-[b,c], b-[d,e]]
  159% ?- rel_to_fun([], R).
  160% ?- rel_to_fun([a-b, a-c, b-d, b-e], R).
  161rel_to_fun(L, R):- sort(L, L0), rel_to_fun(L0, [], R).
  162%
  163rel_to_fun([], X, X).
  164rel_to_fun([A-B|L], [A-U|V], R):-!,
  165	rel_to_fun(L, [A-[B|U]|V], R).
  166rel_to_fun([A-B|L], U, R):-!,
  167	rel_to_fun(L, [A-[B]|U], R).
  168
  169% ?- domain_of_links([a-b, b-c, a-c], Y).
  170domain_of_links(X, Y):-
  171	findall(A , (	member(L, X),
  172					( L = (A - _)
  173					; L = (_ - A)
  174					)
  175				),
  176		   Y0),
  177   sort(Y0, Y).
  178
  179% ?- zdd node_id(a, 0, C).
  180node_id(N, C, C0):- node_id(N, _, C, C0).
  181
  182% ?- time((numlist(1, 100000, Ns), foldl(pred(S, ([I, C, C0]:-
  183%	node_id(st(I), K, C, C0, S))), Ns, 0, R))).
  184node_id(N, I, C, C0):- memo(node_id(N)-I),
  185	(	nonvar(I) -> C0 = C
  186	;	C0 is C+1,
  187		I = C0
  188	).
  189
  190% ?- intern_links([a-b, c-d], R).
  191% ?- intern_links([a-b, b-a], R).
  192intern_links(L, L0):- intern_links(L, L0, 0, _).
  193%
  194intern_links([], [], C, C).
  195intern_links([A-B|L], [I-J|M], C, D):-
  196	node_id(A, I, C, C0),
  197	node_id(B, J, C0, C1),
  198	intern_links(L, M, C1, D).
  199
  200		/******************
  201		*     frontier    *
  202		******************/
 on_frontier(+I, +J, +K) is det
True if node I is accessible directly by a link from a node less than J.
  208% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X), on_frontier(1, 3, X).
  209on_frontier(I, J, F):- arg(I, F, K), J > K.
 on_frontier(+I, +J, +K) is det
True if node I is not accessible directly by a link from a node less than J.
  215% ?- setup_links_frontier(3,[1-2,2-3], _F), write(_F),
  216%	off_frontier(3, 1, _F).
  217% ?- setup_links_frontier(3,[1-2,2-3], _F), write(_F),
  218%	off_frontier(3, 2, _F).
  219
  220% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X), off_frontier(1, 3, X). %false
  221off_frontier(I, J, F):- arg(I, F, K), J =< K.
  222
  223% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X).
  224setup_frontier([], _).
  225setup_frontier([I-J|L], F):-
  226	update_frontier(I, J, F),
  227	!,
  228	setup_frontier(L, F).
  229
  230% ?- X=f(1,2), update_frontier(1,2, X).
  231% ?- X=f(1,2,3), update_frontier(2, 3, X), update_frontier(1, 2, X).
  232update_frontier(I, J, V):-
  233	arg(J, V, A),
  234	(	I < A -> setarg(J, V, I)
  235	;	true
  236	).
  237
  238		/*******************
  239		*	 Helpers       *
  240		*******************/
  241
  242% ?- arrow_symbol(_->_, F).
  243% ?- arrow_symbol(a->b, F, X, Y).
  244arrow_symbol( _ -> _).
  245%
  246arrow_symbol(A, A0):- functor(A, A0, 2).
  247arrow_symbol(A, A0, A1, A2):- functor(A, A0, 2),
  248		arg(1, A, A1),
  249		arg(2, A, A2).
  250
  251% The most basic helpers.
  252composable_pairs(A-B, A-C, B, C).
  253composable_pairs(A-B, C-A, B, C).
  254composable_pairs(B-A, A-C, B, C).
  255composable_pairs(B-A, C-A, B, C).
  256%
  257normal_pair(A-B, B-A):- B < A, !.
  258normal_pair(X, X).
  259%
  260strong_less_than(_-A, B-_):- A<B.
  261
  262		/************************
  263		*     core predicates   *
  264		************************/
  265
  266% ?- rect_path_count(rect(3,3), C).
  267
  268add_links([], X, X):-!.
  269add_links([A-Ns|Ls], X, Y):-!,
  270	cofact(X0, t( mp(A-A,1), 0, X)),
  271	add_links(A, Ns, X0, X1),
  272	prune_by_frontier(A, X1, X2),
  273%	zdd_gc(X2, X3),
  274	add_links(Ls, X2, Y).
  275%
  276add_links(_, [], X, X):-!.
  277add_links(A, [B|Ns], X, Y):-
  278	add_link(A-B, X, X0),
  279	zdd_join(X, X0, X1),
  280	add_links(A, Ns, X1, Y).
  281%
  282add_link(_, X, 0):- X<2, !.
  283add_link(U, X, Y):-
  284	cofact(X, t(MP, L, R)),
  285	MP = mp(M, _),
  286	add_link(U, L, L0),
  287	(	U = M  -> R0 = 0	% cycle found
  288	;   strong_less_than(U, M) -> R0 = 0  %
  289	;	(	composable_pairs(U, M, V, W) ->
  290			subst_node(MP, V, W, R, R0)
  291		;	add_link(U, R, R1),
  292			zdd_insert(MP, R1, R0)
  293		)
  294	),
  295	zdd_join(L0, R0, Y).
  296
  297% ?- rect_path_count(rect(2,2), C).
  298% ?- rect_path_count(rect(2,3), C).
  299% ?- profile(rect_path_count(rect(3,3), C)).
  300% ?- rect_path_count(rect(3,4), C).
  301% ?- profile(call_with_time_limit(300, time(rect_path_count(rect(4,4), C)))).
  302%@ % 2,384,256,849 inferences, 297.572 CPU in 301.708 seconds (99% CPU, 8012366 Lips)
  303%@ ERROR: Unhandled exception: Time limit exceeded
  304%@ % 2,384,256,849 inferences, 300.743 CPU in 305.001 seconds (99% CPU, 7927900 Lips)
  305%@ ERROR: Unhandled exception: Time limit exceeded
  306
  307%@ % 3,548,584,576 inferences, 636.603 CPU in 647.391 seconds (98% CPU, 5574248 Lips)
  308%@ C = 8512.
  309
  310subst_node(_, _, _, X, 0):- X<2, !.
  311subst_node(MP, A, P, X, Y):- memo(subst_node(A, P, MP, X)-Y),
  312	(	nonvar(Y) -> true	%		, write(.)
  313	;	cofact(X, t(MP0, L, R)),	% replace A with P
  314		subst_node(MP, A, P, L, L0),
  315		MP = mp(M, Ps),
  316		MP0 = mp(Lu-Ru, Ps0),
  317		(	A < Lu -> R0 = 0		%, write(.)
  318		;	(	Ru = A	->
  319				normal_pair(Lu-P, V),
  320				zdd_merge(Ps, Ps0, Ps1),
  321				zdd_insert(M, Ps1, Ps2),
  322				zdd_insert(mp(V, Ps2), R, R0)
  323			;	Lu = A	->
  324				normal_pair(P-Ru, V),
  325				zdd_merge(Ps, Ps0, Ps1),
  326				zdd_insert(M, Ps1, Ps2),
  327				zdd_insert(mp(V, Ps2), R, R0)
  328			;	subst_node(MP, A, P, R, R1),
  329				zdd_insert(MP0, R1, R0)
  330			)
  331		),
  332		zdd_join(L0, R0, Y)
  333	).
  334
  335		/***************************
  336		*     Prune by frontier    *
  337		***************************/
  338
  339prune_by_frontier(I, X, Y):-
  340	get_key(frontier, V),
  341	get_key(ends, M),
  342	prune_by_frontier(X, Y, I, M, V).
 prune_by_frontier(+X, -Y, +I, +M, +V) is det
Y is unified with pruned X. 1) A-A is removed from path that has A-A when A is off_frontier. 2) Path which has A-B with off_frontier A or B is removed.
  349prune_by_frontier(X, X, _, _, _):- X<2, !.
  350prune_by_frontier(X, Y, I, M, V):- cofact(X, t(MP, L, R)),
  351	MP = mp(A, _),
  352	classify_pair(A, I, M, V, C),
  353	prune_by_frontier(L, L0, I, M, V),
  354	(	C = keep ->
  355		prune_by_frontier(R, R1, I, M, V),
  356		zdd_insert(MP, R1, R0)
  357	;	C = ignore ->
  358		prune_by_frontier(R, R0, I, M, V)
  359	;	R0 = 0
  360	),
  361	zdd_join(L0, R0, Y).
  362
  363% helper.
  364on_pair(J, J-_).
  365on_pair(J, _-J).
  366
  367% works! [2024/01/11]
  368classify_pair(J-J, I, Pair, V, C):-!,
  369	(	on_frontier(J, I, V) -> C = keep
  370	;	(	on_pair(J, Pair) -> C = 0
  371		;	C = ignore
  372		)
  373	).
  374classify_pair(J-K, I, Pair, V, C):-
  375	(	on_frontier(J, I, V) ->
  376		(	on_frontier(K, I, V) -> C = keep
  377		;	(	on_pair(K, Pair) -> C = keep
  378			;	C = 0
  379			)
  380		)
  381	;	(	on_frontier(K, I, V) -> C = keep
  382		;	(	on_pair(J, Pair) -> C = keep
  383			;	C = 0
  384			)
  385		)
  386	).
  387
  388% ?- zdd, X<< +[*[a-b, a->b]], prune_final(a, b, X, Y), psa(X), psa(Y).
  389prune_final(_, _, X, 0):- X<2, !.
  390prune_final(P, Q, X, Y):- cofact(X, t(MP, L, R)),
  391	MP = mp(A, Ps),
  392	prune_final(P, Q, L, L0),
  393	(  	A = P-Q ->
  394		(	prune_final0(R)->
  395			R0 = Ps
  396		;	R0 = 0
  397		)
  398	;	A = V-V -> prune_final(P, Q, R, R0)
  399	;	R0 = 0
  400	),
  401	zdd_join(L0, R0, Y).
  402%
  403prune_final0(1):-!.
  404prune_final0(X):- X>2,
  405	cofact(X, t(MP, L, R)),
  406	(	prune_final0(L) -> true
  407	;	MP = mp(V-V, _) ->
  408		prune_final0(R)
  409	)