count_path_by_degree

    1:- module(cpbdc, []).    2:- use_module(zdd('zdd-array')).    3:- use_module(zdd(zdd)).    4:- use_module(pac(op)).    5
    6
    7		/************************************************
    8		%     Count the number of paths					*
    9		%	?- rect_path_count(rect(2,2), C).			*
   10		************************************************/
   11
   12% ?- compare0(C, a-b, n(3,4,5)).
   13% ?- compare0(C, n(3,4,5), a-b).
   14% ?- compare0(C, a-b, a-c).
   15% ?- predsort(compare0, [a, a-b, n(1,2,3)], R).
   16
   17compare0(<, n(_,_,_), _-_):-!.
   18compare0(>, _-_, n(_,_,_)):-!.
   19compare0(C, X, Y):- functor(X,F,N), functor(Y, F, N), !,
   20		compare(C, X, Y).
   21compare0(>, _, n(_,_,_)):-!.
   22compare0(>, _, _-_):-!.
   23compare0(<, n(_,_,_),_):-!.
   24compare0(<, _-_, _).
   25
   26% ?- set_compare(compare0).
   27% ?- get_compare(X).
   28
   29% ?- rect_path_count(rect(1,0), C).
   30% ?- rect_path_count(rect(0,1), C).
   31% ?- rect_path_count(rect(1,1), C).
   32% ?- rect_path_count(rect(2,2), C).
   33% ?- rect_path_count(rect(3,3), C).
   34% ?- rect_path_count(rect(4,4), C).
   35% ?- profile(call_with_time_limit(200, rect_path_count(rect(5,5), C))).
   36%@ ERROR: Unhandled exception: Time limit exceeded
   37
   38% ?- rect_path_count(rect(7,7), C).
   39% ?- (X <<( *[n(a,0,1), n(b, 0, 2)])), add_link(a-b, X, Y), psa(Y).
   40% ?- X << *[n(a,0,1), n(b, 0, 2), n(c, 0, 3)],
   41%	add_link(a-b, X, Y), add_link(b-c, Y, Z), psa(Y), psa(Z).
   42% ?- X << *[n(a,0,1), n(b, 0, 2), n(c, 0, 3), n(d, 0, 4)],
   43%	add_link(a-b, X, Y), add_link(b-c, Y, Z), add_link(c-d, Z, U),
   44% psa(Y), psa(Z), psa(U).
   45
   46% ?- X << *[n(a,0,1), n(b, 0, 2), n(c, 0, 3), n(d, 0, 4)],
   47%	add_link(c-d, X, Y), add_link(b-c, Y, Z), add_link(a-b, Z, U),
   48%	psa(Y), psa(Z), psa(U).
   49
   50% ?- X << +[*[n(a,0,1), n(b, 0, 2), n(c, 0, 3), n(d, 0, 3)],
   51%			*[n(a,0,1), n(b, 1, 2), n(c, 0, 3), n(d, 0, 4)]],
   52%	add_link(c-d, X, Y), add_link(b-c, Y, Z), add_link(a-b, Z, U),
   53%	psa(Y), psa(Z), psa(U).
   54
   55% ?- X <<( *[n(a,0,1), n(b, 0, 2)]), add_link(a-b, X, Y).
   56% ?- X <<( *[n(a,0,1), n(b, 0, 2)]), psa(X).
   57
   58udg_path(End, PathSet):- get_key(coa, Coa),
   59	set_key(ends, End),
   60	udg_mate_prune(Coa, 1, PathSet).
   61
   62		/******************************
   63		*     counting path in UDG    *
   64		******************************/
   65% ?- zdd.
   66% ?- udg_path_count(a*d, [a-b, b-c, c-d], C).  % fail.
   67% ?- udg_path_count(a-d, [a-b, b-c, c-d], C).
   68% ?- udg_path_count(a-x, [a-b, b-c, c-d], C).  % fail.
   69
   70% A big circular graph.
   71% ?- zdd.
   72% ?- N = 1000,
   73%	findall(A,	( between(1, N, I), I<N, J is I + 1, A=(I-J)), Ls),
   74%	time(udg_path_count(1-2,[1-N |Ls], C)).
   75%@ % 87,134,142 inferences, 41.767 CPU in 42.025 seconds (99% CPU, 2086215 Lips)
   76
   77% Counting paths in complete graph.
   78% ?- N = 11, findall(I-J, ( between(1, N, I), between(1, N, J), I < J), Ls),
   79%  time(zdd udg_path_count(1-N, Ls, C)).
   80%@ % 1,835,060,694 inferences, 233.677 CPU in 238.492 seconds (98% CPU, 7852993 Lips)
   81%@ N = 11,
   82%@ Ls = [1-2, 1-3, 1-4, 1-5, 1-6, 1-7, 1-8, 1-9, ... - ...|...],
   83%@ C = 986410.
   84
   85% [2024/01/04] Is global variable efficient ?
   86% ?- N = 11, findall(I-J, ( between(1, N, I), between(1, N, J), I < J), Ls),
   87%  time(zdd udg_path_count(1-N, Ls, C)).
   88%@ % 1,775,162,524 inferences, 211.846 CPU in 216.653 seconds (98% CPU, 8379506 Lips)
   89%@ N = 11,
   90%@ Ls = [1-2, 1-3, 1-4, 1-5, 1-6, 1-7, 1-8, 1-9, ... - ...|...],
   91%@ C = 986410.
   92%
   93udg_path_count(Ends, Links, C):- udg_path_count(Ends, Links, C, _).
   94%
   95udg_path_count(Ends, Links, C, X):-
   96		prepare_udg(Ends, Links),
   97		!,
   98		get_key(coa, Coa),
   99		udg_mate_prune(Coa, 1, X),
  100		card(X, C).
  101%
  102udg_mate_prune(Ls, X, Y):-
  103	add_links(Ls, X, Y0),
  104	prune_final(Y0, Y).
  105
  106% ?- rect_path_count(rect(1,0), C).
  107% ?- rect_path_count(rect(0,1), C).
  108% ?- rect_path_count(rect(1,1), C).
  109% ?- rect_path_count(rect(2,1), C).
  110% ?- rect_path_count(rect(2,2), C).
  111% ?- rect_path_count(rect(3,3), C).
  112% ?- profile(rect_path_count(rect(4,4), C)).
  113% ?- profile(call_with_time_limit(200, rect_path_count(rect(5,5), C))).
  114% ?- udg_path_count(a-c,[a-b, b-c, a-c], C).
  115% ?- call_with_time_limit(120, rect_path_count(rect(5,5), C)). % time out
  116
  117rect_path_count(R, C):- R = rect(W, H),
  118	rect_path_count( p(0,0) - p(W,H), R, C, _).
  119%
  120rect_path_count(Ends, R, C, X):- rect_links(R, Links),
  121	udg_path_count(Ends, Links, C, X).
  122
  123% ?- rect_links(rect(1,1), Links).
  124% ?- rect_links(rect(10,10), Links),length(Links, L).
  125rect_links(rect(W, H), Links):-
  126	findall( p(I,J) - p(K,L),
  127				 (	between(0, W, I),
  128					between(0, H, J),
  129					(  L = J, K is I + 1, K =< W
  130					;  K = I, L is J + 1, L =< H
  131					)
  132				 ),
  133				 Links).
  134
  135		/********************************
  136		*     Prepare UDG in coalgebra  *
  137		********************************/
  138
  139% ?- udg_path_count(a-b,[a-b], C).
  140% ?- udg_path_count(a-c,[a-b, b-c], C).
  141% ?- udg_path_count(a-c,[a-b, b-c, a-c], C).
  142
  143% ?- N=2, findall(I-J, (between(1, N, I), J is I + 1, J=<N ), Ls),
  144%	udg_path_count(1-N, Ls, C).
  145% ?- N=10, findall(I-J, (between(1, N, I), J is I + 1, J=<N ), Ls),
  146%	udg_path_count(1-N, Ls, C).
  147% ?- N=1500, findall(I-J, (between(1, N, I), J is I + 1, J=<N ), Ls),
  148%	udg_path_count(1-N, Ls, C).
  149% ?- N=1500, findall(I-J, (between(1, N, I), J is I + 1, J=<N ), Ls),
  150%	udg_path_count(1-2, Ls, C).
  151
  152% ?- udg_path_count(a-c,[a-b, a-d, b-c, d-c], C).
  153% ?- udg_path_count(a-f,[a-b, a-d, b-c, d-c, d-e, e-f, c-f], C).
  154% ?- udg_path_count(a-k,
  155%	[a-b, a-d, b-c,
  156%	d-c, d-e, e-f,
  157%	c-f, b-g, g-h,
  158%	h-k, c-h, f-k], C).
  159% ?- prepare_udg([a-b, b-c, c-d]), memo(node_id(a)-Id), get_key(dom, X).
  160
  161% ?- prepare_udg(a-d, [a-b, b-c, c-d]),
  162%	get_key(ends, ST),
  163%	get_key(frontier, F),
  164%	get_key(dom, Dom),
  165%	get_key(coa, U).
  166
  167prepare_udg(ST, Links):-
  168	prepare_udg(Links),
  169	prepare_ends(ST, A-B),
  170	set_key(ends, A-B).
  171%
  172prepare_udg(Links):-
  173	prepare_udg(Links, Succs, D, Vec),
  174	length(D, N),
  175	completing_succs(Succs, Succs0, N),
  176	set_key(coa, Succs0),
  177	set_key(dom, D),
  178	set_key(frontier, Vec).
  179%
  180prepare_udg(Links, Succs, D, Vec):-
  181	intern_links(Links, Links0),
  182	normal_mate_list(Links0, Links1),
  183	sort(Links1, Links2),
  184	domain_of_links(Links2, D),
  185	rel_to_fun(Links2, Succs),
  186	Vec=..[#|D],
  187	setup_frontier(Links1, Vec).
  188%
  189prepare_ends(A-B, R):-!, R = A0-B0,
  190	memo(node_id(A)-I),
  191	memo(node_id(B)-J),
  192	(	nonvar(I), nonvar(J) -> sort([I, J], [A0, B0])
  193	;	format("No link at ~w or ~w\n", [A,B]),
  194		fail
  195	).
  196prepare_ends(E, _):-
  197	format("Unexpected form of end nodes ~w \n", [E]),
  198	fail.
  199
  200% ?-completing_succs([], Y, 2).
  201% ?-completing_succs([2-[a]], Y, 3).
  202completing_succs(X, X, 0):-!.
  203completing_succs([I-A|Ls], [I-A|Ms], I):-!, J is I - 1,
  204	completing_succs(Ls, Ms, J).
  205completing_succs(Ls, [I-[]|Ms], I):- J is I - 1,
  206	completing_succs(Ls, Ms, J).
  207
  208% ?- normal_mate_list([1-2], X).
  209% ?- normal_mate_list([2-1, 1-2], X).
  210normal_mate_list([], []).
  211normal_mate_list([P|R], [P0|R0]):- P = I-J,
  212	(	J @< I -> P0 = J - I
  213	;	P0 = P
  214	),
  215	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]]

  223% ?- rel_to_fun([], R).
  224% ?- rel_to_fun([a-b, a-c, b-d, b-e], R).
  225rel_to_fun(L, R):- sort(L, L0), rel_to_fun(L0, [], R).
  226%
  227rel_to_fun([], X, X).
  228rel_to_fun([A-B|L], [A-U|V], R):-!,
  229	rel_to_fun(L, [A-[B|U]|V], R).
  230rel_to_fun([A-B|L], U, R):-!,
  231	rel_to_fun(L, [A-[B]|U], R).
  232
  233% ?- domain_of_links([a-b, b-c, a-c], Y).
  234domain_of_links(X, Y):-
  235	findall(A , (	member(L, X),
  236					( L = (A - _)
  237					; L = (_ - A)
  238					)
  239				),
  240		   Y0),
  241   sort(Y0, Y).
  242
  243% ?- node_id(a, 0, C).
  244node_id(N, C, C0):- node_id(N, _, C, C0).
  245
  246% c.f.
  247% ?- foldl(pred([A, B, [A|B]]), [a,b,c], [], R).
  248% ?- foldl(pred(S, [A, B, [S, A|B]]), [a,b,c], [], U), S=1.
  249% ?- foldr(pred([A, B, [A|B]]), [a,b,c], [], R).
  250% ?- show(foldr(pred([A, B, [A|B]]), [a,b,c], [], R)).
  251% ?- foldr(pred(S, [A, B, [S, A|B]]), [a,b,c], [], U), S=1.
  252
  253node_id(N, I, C, C0):- memo(node_id(N)-I),
  254	(	nonvar(I) -> C0 = C
  255	;	C0 is C+1,
  256		I = C0
  257	).
  258
  259% ?- intern_links([a-b, c-d], R).
  260% ?- intern_links([a-b, b-a], R).
  261intern_links(L, L0):- intern_links(L, L0, 0, _).
  262%
  263intern_links([], [], C, C).
  264intern_links([A-B|L], [I-J|M], C, D):-
  265	node_id(A, I, C, C0),
  266	node_id(B, J, C0, C1),
  267	intern_links(L, M, C1, D).
  268
  269		/******************
  270		*     frontier    *
  271		******************/

on_frontier(+I, +J, +K) is det: True if node I is accessible directly by a link from a node less than J.

  277% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X), on_frontier(1, 3, X).
  278%@ X = f(1, 1, 2).
  279on_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.

  285% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X), off_frontier(1, 3, X). %false
  286off_frontier(I, J, F):- arg(I, F, K), J =< K.
  287
  288% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X).
  289setup_frontier([], _).
  290setup_frontier([I-J|L], F):-
  291	update_frontier(I, J, F),
  292	!,
  293	setup_frontier(L, F).
  294
  295% ?- X=f(1,2), update_frontier(1,2, X).
  296% ?- X=f(1,2,3), update_frontier(2, 3, X), update_frontier(1, 2, X).
  297update_frontier(I, J, V):-
  298	arg(J, V, A),
  299	(	I < A -> setarg(J, V, I)
  300	;	true
  301	).
  302
  303
  304		/************************
  305		*     core predicates   *
  306		************************/
  307add_links([], X, X):-!.
  308add_links([A-Ns|Ls], X, Y):-!,
  309	cofact(X0, t(n(A,0,A), 0, X)),
  310	add_links(A, Ns, X0, X1),
  311	prune_by_frontier(A, X1, X2),
  312	add_links(Ls, X2, Y).
  313%
  314add_links(_, [], X, X):-!.
  315add_links(A, [B|Ns], X, Y):-
  316	b_setval(current_link, A-B),
  317	add_link(X, X0),
  318	zdd_join(X, X0, X1),
  319	add_links(A, Ns, X1, Y).
  320
  321% ?- spy(cleanup_dot_star).
  322% ?- spy(subst_class_id).
  323% ?- rect_path_count(rect(1,0), C).
  324% ?- rect_path_count(rect(0, 1), C).
  325% ?- rect_path_count(rect(1, 1), C).
  326% ?- rect_path_count(rect(4, 4), C).
  327%@ C = 8512.
  328% ?- time(rect_path_count(rect(5,5), C)).
  329%@ % 13,151,026,211 inferences, 2239.205 CPU in 2269.526 seconds (99% CPU, 5873078 Lips)
  330%@ C = 1262816.
  331
  332add_link(X, 0):- X<2, !.
  333add_link(X, Y):-
  334	cofact(X, t(A, L, R)),
  335	add_link(L, L0),
  336	b_getval(current_link, I-J),
  337	A = n(K,G,C),
  338	(	I @< K -> R0 = 0
  339	;	G = 2 -> R0  = 0
  340	;	K = I ->
  341		update_class_degree(J, C, R, R1),
  342		cleanup_dot_star(R1,  R2),
  343		G1 is G + 1,
  344		A1 = n(I, G1, C),
  345		zdd_insert(A1, R2, R0)
  346	;	add_link(R, R1),
  347		zdd_insert(A, R1, R0)
  348	),
  349	zdd_join(L0, R0, Y).
  350
  351%
  352update_class_degree(_, _, X, 0):- X < 2, !.
  353update_class_degree(J, C, X, Y):- cofact(X, t(V, L, R)),
  354	(	V=(_-_) -> Y = X
  355	;	update_class_degree(J, C, L, L0),
  356		V = n(K, G, C0),
  357		(	J = K ->
  358			(   G = 2 -> R0 = 0		% must not be a path.
  359			;	C = C0 -> R0 = 0	% link in a already merged cluster
  360			;	G1 is G + 1,
  361				subst_class_id(C0, C, R, R1),  % union cluster by a link
  362				b_getval(current_link, Link),
  363				zdd_ord_insert([n(K, G1, C),Link], R1, R2),
  364				cofact(R0, t(*, change(C0, C), R2))
  365			)
  366		;	update_class_degree(J, C, R, R1),
  367			insert_through_dot(V, R1, R0)
  368		),
  369		cofact(Y, t((.), L0, R0))
  370	).
  371%
  372subst_class_id(_, _, X, X):- X<2, !.
  373subst_class_id(C, D, X, Y):- memo(subst_cla_id(C,D,X)-Y),
  374	(	nonvar(Y)->true				%, write(.)
  375	;	cofact(X, t(U,L,R)),
  376		(	U = (_-_) -> Y = X
  377		;	subst_class_id(C, D, L, L0),
  378			subst_class_id(C, D, R, R0),
  379			U = n(I, G, C0),
  380			(   C = C0 -> C1 = D
  381			;   C1 = C0
  382			),
  383			cofact(Y, t(n(I, G, C1), L0, R0))
  384		)
  385	).
  386%
  387cleanup_dot_star(X, X):- X<2, !.
  388cleanup_dot_star(X, Y):- cofact(X, U),
  389	cleanup_dot_star_case(U, Y).
  390%
  391cleanup_dot_star_case(t(., L, R), V):-!,
  392	cleanup_dot_star(R, R0),
  393	cleanup_dot_star(L, L0),
  394	zdd_join(L0, R0, V).
  395cleanup_dot_star_case(t(*, _, R), R):-!.
  396cleanup_dot_star_case(X, Y):- cofact(Y, X).
  397
  398%
  399insert_through_dot(_, X, 0):- X<2, !.
  400insert_through_dot(A, X, Y):- cofact(X, T),
  401	T = t(U, L, R),
  402	(	U = (.) ->
  403		insert_through_dot(A, L, L0),
  404		insert_through_dot(A, R, R0),
  405		cofact(Y, t(U, L0, R0))
  406	;	U = (*) ->
  407		insert_aside_star(A, T, Y)
  408	).
  409%
  410insert_aside_star(n(I, Deg, C), T,  Y):-
  411	T = t(*, change(C0, C1), R),
  412	( C = C0 -> N = n(I, Deg, C1)
  413	; N = n(I, Deg, C)
  414	),
  415	zdd_insert(N, R, R0),
  416	cofact(Y, t(*, change(C0, C1), R0)).
  417
  418		/***************************
  419		*     Prune by frontier    *
  420		***************************/
  421
  422prune_by_frontier(I, X, Y):-
  423	get_key(frontier, V),
  424	get_key(ends, M),
  425	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.

  432prune_by_frontier(X, X, _, _, _):- X<2, !.
  433prune_by_frontier(X, Y, I, M, V):- memo(prune_by_frontier(X, I)-Y),
  434	(	nonvar(Y)->true, write(.);
  435	cofact(X, t(A, L, R)),
  436	(	A = (_-_) -> Y = X
  437	;	classify_triple(A, I, M, V, C),
  438		prune_by_frontier(L, L0, I, M, V),
  439		(	C = keep ->
  440			prune_by_frontier(R, R1, I, M, V),
  441			zdd_insert(A, R1, R0)
  442		;	R0 = 0
  443		),
  444		zdd_join(L0, R0, Y)
  445	)).
  446
  447
  448% helper.
  449on_pair(J, J-_).
  450on_pair(J, _-J).
  451
  452%
  453classify_triple(n(J, Deg, _), I, Pair, V, C):-!,
  454	(	on_frontier(J, I, V) -> C = keep
  455	;	(	on_pair(J, Pair) ->
  456			(	Deg = 1 -> C = keep
  457			;	C = 0
  458			)
  459		;	(	Deg = 1 -> C = 0
  460			;	C = keep
  461			)
  462		)
  463	).
  464%
  465prune_final(X, Y):-
  466	get_key(ends, Pair),
  467	prune_final(Pair, X, Y).
  468
  469%
  470prune_final(_, X, X):- X<2, !.
  471prune_final(Pair, X, Y):- cofact(X, t(A, L, R)),
  472	(	A= (_-_) -> Y = X
  473	;	prune_final(Pair, L, L0),
  474		A = n(_, _, C),
  475		prune_final(C, Pair, R, R1),
  476		zdd_insert(A, R1, R0),
  477		zdd_join(L0, R0, Y)
  478	).
  479%
  480prune_final(_, _, X, X):- X<2, !.
  481prune_final(C0, Pair, X, Y):- cofact(X, t(A, L, R)),
  482	(	A = (_-_) -> Y = X
  483 	;   prune_final(C0, Pair, L, L0),
  484		A = n(J, Deg, C),
  485		( 	on_pair(J, Pair) ->
  486			(	Deg = 1 -> prune_final(C0, Pair, R, R0)
  487			;   R0  = 0
  488			)
  489		;	(	Deg = 1 -> R0 = 0
  490			;	Deg = 2 ->
  491				(	C0 = C ->
  492					prune_final(C, Pair, R, R1),
  493					zdd_insert(A, R1, R0)
  494				;	R0 = 0
  495				)
  496			;	prune_final(C0, Pair, R, R1),
  497				zdd_insert(A, R1, R0)
  498			)
  499		),
  500		zdd_join(L0, R0, Y)
  501	)