count-region-by-degree.pl

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

  200% ?- rel_to_fun([], R).
  201% ?- rel_to_fun([a-b, a-c, b-d, b-e], R).
  202rel_to_fun(L, R):- sort(L, L0), rel_to_fun(L0, [], R).
  203%
  204rel_to_fun([], X, X).
  205rel_to_fun([A-B|L], [A-U|V], R):-!,
  206	rel_to_fun(L, [A-[B|U]|V], R).
  207rel_to_fun([A-B|L], U, R):-!,
  208	rel_to_fun(L, [A-[B]|U], R).
  209
  210% ?- domain_of_links([a-b, b-c, a-c], Y).
  211domain_of_links(X, Y):-
  212	findall(A , (	member(L, X),
  213					( L = (A - _)
  214					; L = (_ - A)
  215					)
  216				),
  217		   Y0),
  218   sort(Y0, Y).
  219
  220% ?- node_id(a, 0, C).
  221node_id(N, C, C0):- node_id(N, _, C, C0).
  222
  223% c.f.
  224% ?- foldl(pred([A, B, [A|B]]), [a,b,c], [], R).
  225% ?- foldl(pred(S, [A, B, [S, A|B]]), [a,b,c], [], U), S=1.
  226% ?- foldr(pred([A, B, [A|B]]), [a,b,c], [], R).
  227% ?- show(foldr(pred([A, B, [A|B]]), [a,b,c], [], R)).
  228% ?- foldr(pred(S, [A, B, [S, A|B]]), [a,b,c], [], U), S=1.
  229
  230node_id(N, I, C, C0):- memo(node_id(N)-I),
  231	(	nonvar(I) -> C0 = C
  232	;	C0 is C+1,
  233		I = C0
  234	).
  235
  236% ?- intern_links([a-b, c-d], R).
  237% ?- intern_links([a-b, b-a], R).
  238intern_links(L, L0):- intern_links(L, L0, 0, _).
  239%
  240intern_links([], [], C, C).
  241intern_links([A-B|L], [I-J|M], C, D):-
  242	node_id(A, I, C, C0),
  243	node_id(B, J, C0, C1),
  244	intern_links(L, M, C1, D).
  245
  246		/******************
  247		*     frontier    *
  248		******************/

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

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

  262% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X), off_frontier(1, 3, X). %false
  263off_frontier(I, J, F):- arg(I, F, K), J =< K.
  264
  265% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X).
  266setup_frontier([], _).
  267setup_frontier([I-J|L], F):-
  268	update_frontier(I, J, F),
  269	!,
  270	setup_frontier(L, F).
  271
  272% ?- X=f(1,2), update_frontier(1,2, X).
  273% ?- X=f(1,2,3), update_frontier(2, 3, X), update_frontier(1, 2, X).
  274update_frontier(I, J, V):-
  275	arg(J, V, A),
  276	(	I < A -> setarg(J, V, I)
  277	;	true
  278	).
  279
  280
  281		/************************
  282		*     core predicates   *
  283		************************/
  284add_links([], X, X):-!.
  285add_links([A-Ns|Ls], X, Y):-!,
  286	cofact(X0, t(n(A,0,A), 0, X)),
  287	add_links(A, Ns, X0, X1),
  288	prune_by_frontier(A, X1, X2),
  289	add_links(Ls, X2, Y).
  290%
  291add_links(_, [], X, X):-!.
  292add_links(A, [B|Ns], X, Y):-
  293	add_link(A-B, X, X0),
  294	zdd_join(X, X0, X1),
  295	add_links(A, Ns, X1, Y).
  296%
  297add_link(_, X, 0):- X<2, !.
  298add_link(I-J, X, Y):-
  299	cofact(X, t(A, L, R)),
  300	add_link(I-J, L, L0),
  301	A = n(K,G,C),
  302	(	I @< K -> R0 = 0
  303	;	G = 2 -> R0  = 0
  304	;	K = I ->
  305		update_class_degree(J, C, R, R1),
  306		cleanup_dot_star(R1,  R2),
  307		G1 is G + 1,
  308		A1 = n(I, G1, C),
  309		zdd_insert(A1, R2, R0)
  310	;	add_link(I-J, R, R1),
  311		zdd_insert(A, R1, R0)
  312	),
  313	zdd_join(L0, R0, Y).
  314%
  315update_class_degree(_, _, X, 0):- X < 2, !.
  316update_class_degree(J, C, X, Y):- cofact(X, t(V, L, R)),
  317	update_class_degree(J, C, L, L0),
  318	V = n(K, G, C0),
  319	(	J = K ->
  320		(   G = 2 -> R0 = 0		% must not be a path.
  321		;	C = C0 -> R0 = 0	% link in a already merged cluster
  322		;	G1 is G + 1,
  323			subst_class_id(C0, C, R, R1),  % union cluster by a link
  324			zdd_insert(n(K, G1, C), R1, R2),
  325			cofact(R0, t(*, change(C0, C), R2))
  326		)
  327	;	update_class_degree(J, C, R, R1),
  328		insert_through_dot(V, R1, R0)
  329	),
  330	cofact(Y, t((.), L0, R0)).
  331%
  332subst_class_id(_, _, X, X):-X<2,!.
  333subst_class_id(C, D, X, Y):- cofact(X, t(U,L,R)),
  334	subst_class_id(C, D, L, L0),
  335	subst_class_id(C, D, R, R0),
  336	U = n(I, G, C0),
  337	(   C = C0 -> C1 = D
  338	;   C1 = C0
  339	),
  340	cofact(Y, t(n(I, G, C1), L0, R0)).
  341%
  342cleanup_dot_star(X, X):- X<2, !.
  343cleanup_dot_star(X, Y):- cofact(X, U),
  344	cleanup_dot_star_case(U, Y).
  345%
  346cleanup_dot_star_case(t(., L, R), V):-!,
  347	cleanup_dot_star(R, R0),
  348	cleanup_dot_star(L, L0),
  349	zdd_join(L0, R0, V).
  350cleanup_dot_star_case(t(*, _, R), R):-!.
  351cleanup_dot_star_case(X, X).
  352
  353%
  354insert_through_dot(_, X, 0):- X<2, !.
  355insert_through_dot(A, X, Y):- cofact(X, T),
  356	T = t(U, L, R),
  357	(	U = (.) ->
  358		insert_through_dot(A, L, L0),
  359		insert_through_dot(A, R, R0),
  360		cofact(Y, t(U, L0, R0))
  361	;	U = (*) ->
  362		insert_aside_star(A, T, Y)
  363	).
  364%
  365insert_aside_star(n(I, Deg, C), T,  Y):-
  366	T = t(*, change(C0, C1), R),
  367	( C = C0 -> N = n(I, Deg, C1)
  368	; N = n(I, Deg, C)
  369	),
  370	zdd_insert(N, R, R0),
  371	cofact(Y, t(*, change(C0, C1), R0)).
  372
  373		/***************************
  374		*     Prune by frontier    *
  375		***************************/
  376
  377prune_by_frontier(I, X, Y):-
  378	get_key(frontier, V),
  379	get_key(ends, M),
  380	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.

  387prune_by_frontier(X, X, _, _, _):- X<2, !.
  388prune_by_frontier(X, Y, I, M, V):- cofact(X, t(A, L, R)),
  389	classify_triple(A, I, M, V, C),
  390	prune_by_frontier(L, L0, I, M, V),
  391	(	C = keep ->
  392		prune_by_frontier(R, R1, I, M, V),
  393		zdd_insert(A, R1, R0)
  394	;	R0 = 0
  395	),
  396	zdd_join(L0, R0, Y).
  397
  398
  399% helper.
  400on_pair(J, J-_).
  401on_pair(J, _-J).
  402
  403%
  404classify_triple(n(J, Deg, _), I, Pair, V, C):-!,
  405	(	on_frontier(J, I, V) -> C = keep
  406	;	(	on_pair(J, Pair) ->
  407			(	Deg = 1 -> C = keep
  408			;	C = 0
  409			)
  410		;	(	Deg = 1 -> C = 0
  411			;	C = keep
  412			)
  413		)
  414	).
  415%
  416prune_final(X, Y):-
  417	get_key(ends, Pair),
  418	prune_final(Pair, X, Y).
  419
  420%
  421prune_final(_, X, X):- X<2, !.
  422prune_final(Pair, X, Y):- cofact(X, t(A, L, R)),
  423	prune_final(Pair, L, L0),
  424	A = n(_, _, C),
  425	prune_final(C, Pair, R, R1),
  426	zdd_insert(A, R1, R0),
  427	zdd_join(L0, R0, Y).
  428%
  429prune_final(_, _, X, X):- X<2, !.
  430prune_final(C0, Pair, X, Y):- cofact(X, t(A, L, R)),
  431	prune_final(C0, Pair, L, L0),
  432	A = n(J, Deg, C),
  433	( 	on_pair(J, Pair) ->
  434		(	Deg = 1 -> prune_final(C0, Pair, R, R0)
  435		;   R0  = 0
  436		)
  437	;	(	Deg = 1 -> R0 = 0
  438		;	Deg = 2 ->
  439			(	C0 = C ->
  440				prune_final(C, Pair, R, R1),
  441				zdd_insert(A, R1, R0)
  442			;	R0 = 0
  443			)
  444		;	prune_final(C0, Pair, R, R1),
  445			zdd_insert(A, R1, R0)
  446		)
  447	),
  448	zdd_join(L0, R0, Y)