1:- module(hamiton, []).    2
    3:- use_module(zdd('zdd-array')).    4:- use_module(zdd(zdd)).    5:- use_module(pac(op)).    6		/******************************
    7		*     counting path in UDG    *
    8		******************************/
    9
   10% For test.
   11% ?- trace.
   12% ?- hamilton(a-b,[a-b], C).
   13% ?- hamilton(a-c,[a-b, b-c], C).
   14% ?- hamilton(a-c,[a-b, b-c, a-c], C).
   15
   16hamilton(X, Y, Z):- udg_path_count(X, Y, Z).
   17
   18%
   19udg_path_count(Ends, Links, C):- udg_path_count(Ends, Links, C, _).
   20
   21%
   22udg_path_count(Ends, Links, C, X):-
   23		prepare_udg(Ends, Links),
   24		!,
   25		get_key(coa, Coa),
   26		udg_mate_prune(Coa, 1, X),
   27		card(X, C).
   28%
   29udg_mate_prune(Ls, X, Y):-
   30	add_links(Ls, X, Y0), % psa(Y0),
   31	get_key(ends, A-B),
   32	prune_final(A, B, Y0, Y).
   33
   34%  hamilton paths
   35% ?- zdd.
   36
   37% ?- rect_hamilton(rect(1,1), C).
   38%@ C = 0.
   39% ?- rect_hamilton(rect(1,0), C).
   40%@ C = 1.
   41% ?- rect_hamilton(rect(2,2), C).
   42%@ C = 2.
   43% ?- rect_hamilton(rect(2,4), C).
   44%@ C = 8.
   45% ?- rect_hamilton(rect(2,5), C).
   46%@ C = 16.
   47% ?- rect_hamilton(rect(3,3), C).
   48%@ C = 0.
   49% ?- rect_hamilton(rect(4,4), C).
   50%@ C = 104.
   51% ?- rect_hamilton(rect(5,5), C).
   52%@ C = 0.
   53% ?- rect_hamilton(rect(6,6), C).
   54%@ C = 111712.
   55% ?- time(rect_hamilton(rect(7,7), C)).
   56%@ % 78,168,751 inferences, 8.614 CPU in 9.611 seconds (90% CPU, 9074568 Lips)
   57%@ C = 0.
   58% ?- time(rect_hamilton(rect(8,8), C)).
   59%@ % 420,988,309 inferences, 48.069 CPU in 48.939 seconds (98% CPU, 8757922 Lips)
   60%@ C = 2688307514.
   61%@ % 484,915,416 inferences, 33.085 CPU in 33.630 seconds (98% CPU, 14656544 Lips)
   62%@ C = 2688307514.
   63%@ % 484,915,416 inferences, 44.089 CPU in 44.641 seconds (99% CPU, 10998626 Lips)
   64%@ C = 2688307514.
   65
   66%@ % 499,824,225 inferences, 31.372 CPU in 32.454 seconds (97% CPU, 15932031 Lips)
   67%@ C = 2688307514.
   68
   69rect_hamilton(R, C) :- rect_path_count(R, C).
   70%
   71rect_path_count(R, C):- R = rect(W, H),
   72	rect_path_count( p(0,0) - p(W,H), R, C).
   73%
   74rect_path_count(Ends, R, C):- rect_links(R, Links),
   75	udg_path_count(Ends, Links, C).
   76
   77% ?- rect_links(rect(1,1), Links).
   78% ?- rect_links(rect(10,10), Links),length(Links, L).
   79rect_links(rect(W, H), Links):-
   80	findall( p(I,J) - p(K,L),
   81				 (	between(0, W, I),
   82					between(0, H, J),
   83					(  L = J, K is I + 1, K =< W
   84					;  K = I, L is J + 1, L =< H
   85					)
   86				 ),
   87				 Links).
   88
   89		/********************************
   90		*     Prepare UDG in coalgebra  *
   91		********************************/
   92
   93
   94% ?- zdd udg_path_count(a-c,[a-b, a-d, b-c, d-c], C).
   95% ?- zdd udg_path_count(a-f,[a-b, a-d, b-c, d-c, d-e, e-f, c-f], C).
   96% ?- zdd udg_path_count(a-k,
   97%	[a-b, a-d, b-c,
   98%	d-c, d-e, e-f,
   99%	c-f, b-g, g-h,
  100%	h-k, c-h, f-k], C).
  101
  102% ?- zdd prepare_udg(a-d, [a-b, b-c, c-d]),
  103%	get_key(ends, ST),
  104%	get_key(frontier, F),
  105%	get_key(dom, Dom),
  106%	get_key(coa, U).
  107
  108prepare_udg(ST, Links):-
  109	prepare_udg(Links),
  110	prepare_ends(ST, A-B),
  111	set_key(ends, A-B).
  112%
  113prepare_udg(Links):-
  114	prepare_udg(Links, Succs, D, Vec),
  115	length(D, N),
  116	completing_succs(Succs, Succs0, N),
  117	set_key(coa, Succs0),
  118	set_key(dom, D),
  119	set_key(frontier, Vec).
  120%
  121prepare_udg(Links, Succs, D, Vec):-
  122	intern_links(Links, Links0),
  123	normal_mate_list(Links0, Links1),
  124	sort(Links1, Links2),
  125	domain_of_links(Links2, D),
  126	rel_to_fun(Links2, Succs),
  127	Vec=..[#|D],
  128	setup_frontier(Links1, Vec).
  129%
  130prepare_ends(A-B, R):-!, R = A0-B0,
  131	memo(node_id(A)-I),
  132	memo(node_id(B)-J),
  133	(	nonvar(I), nonvar(J) -> sort([I, J], [A0, B0])
  134	;	format("No link at ~w or ~w\n", [A,B]),
  135		fail
  136	).
  137prepare_ends(E, _):-
  138	format("Unexpected form of end nodes ~w \n", [E]),
  139	fail.
  140
  141% ?-completing_succs([], Y, 2).
  142% ?-completing_succs([2-[a]], Y, 3).
  143completing_succs(X, X, 0):-!.
  144completing_succs([I-A|Ls], [I-A|Ms], I):-!, J is I - 1,
  145	completing_succs(Ls, Ms, J).
  146completing_succs(Ls, [I-[]|Ms], I):- J is I - 1,
  147	completing_succs(Ls, Ms, J).
  148
  149% ?- normal_mate_list([1-2], X).
  150% ?- normal_mate_list([2-1, 1-2], X).
  151normal_mate_list([], []).
  152normal_mate_list([P|R], [P0|R0]):- P = I-J,
  153	(	J @< I -> P0 = J - I
  154	;	P0 = P
  155	),
  156	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]]
  164% ?- rel_to_fun([], R).
  165% ?- rel_to_fun([a-b, a-c, b-d, b-e], R).
  166rel_to_fun(L, R):- sort(L, L0), rel_to_fun(L0, [], R).
  167%
  168rel_to_fun([], X, X).
  169rel_to_fun([A-B|L], [A-U|V], R):-!,
  170	rel_to_fun(L, [A-[B|U]|V], R).
  171rel_to_fun([A-B|L], U, R):-!,
  172	rel_to_fun(L, [A-[B]|U], R).
  173
  174% ?- domain_of_links([a-b, b-c, a-c], Y).
  175domain_of_links(X, Y):-
  176	findall(A , (	member(L, X),
  177					( L = (A - _)
  178					; L = (_ - A)
  179					)
  180				),
  181		   Y0),
  182   sort(Y0, Y).
  183
  184% ?- zdd node_id(a, 0, C).
  185node_id(N, C, C0):- node_id(N, _, C, C0).
  186
  187% ?- zdd time((numlist(1, 100000, Ns), foldl(pred(S, ([I, C, C0]:-
  188%	node_id(st(I), K, C, C0, S))), Ns, 0, R))).
  189node_id(N, I, C, C0):- memo(node_id(N)-I),
  190	(	nonvar(I) -> C0 = C
  191	;	C0 is C+1,
  192		I = C0
  193	).
  194
  195% ?- zdd intern_links([a-b, c-d], R).
  196% ?- zdd intern_links([a-b, b-a], R).
  197intern_links(L, L0):- intern_links(L, L0, 0, _).
  198%
  199intern_links([], [], C, C).
  200intern_links([A-B|L], [I-J|M], C, D):-
  201	node_id(A, I, C, C0),
  202	node_id(B, J, C0, C1),
  203	intern_links(L, M, C1, D).
  204
  205		/******************
  206		*     frontier    *
  207		******************/
 on_frontier(+I, +J, +K) is det
True if node I is accessible directly by a link from a node less than J.
  213% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X), on_frontier(1, 3, X).
  214on_frontier(I, J, F):- arg(I, F, K), K < J.
 on_frontier(+I, +J, +K) is det
True if node I is not accessible directly by a link from a node less than J.
  220% ?- setup_links_frontier(3,[1-2,2-3], _F), write(_F),
  221%	off_frontier(3, 1, _F).
  222% ?- setup_links_frontier(3,[1-2,2-3], _F), write(_F),
  223%	off_frontier(3, 2, _F).
  224
  225% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X).
  226%@ X = f(1, 1, 2).
  227% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X), off_frontier(1, 3, X). %false
  228% ?- zdd,  X=f(1,2,3), setup_frontier([1-2,2-3], X), off_frontier(3, 1, X).
  229% ?- zdd,  X=f(1,2,3), setup_frontier([1-2,2-3], X), off_frontier(1, 1, X).
  230% ?- zdd,  X=f(1,2,3), setup_frontier([1-2,2-3], X), on_frontier(1, 1, X).
  231% ?- zdd,  X=f(1,2,3), setup_frontier([1-2,2-3], X), on_frontier(1, 3, X).
  232
  233
  234
  235off_frontier(I, J, F):- arg(I, F, K), K >= J.
  236
  237% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X).
  238setup_frontier([], _).
  239setup_frontier([I-J|L], F):-
  240	update_frontier(I, J, F),
  241	!,
  242	setup_frontier(L, F).
  243
  244% ?- X=f(1,2), update_frontier(1,2, X).
  245% ?- X=f(1,2,3), update_frontier(2, 3, X), update_frontier(1, 2, X).
  246update_frontier(I, J, V):-
  247	arg(J, V, A),
  248	(	I < A -> setarg(J, V, I)
  249	;	true
  250	).
  251
  252		/*******************
  253		*	 Helpers       *
  254		*******************/
  255
  256% ?- arrow_symbol(_->_, F).
  257% ?- arrow_symbol(a->b, F, X, Y).
  258arrow_symbol( _ -> _).
  259%
  260arrow_symbol(A, A0):- functor(A, A0, 2).
  261arrow_symbol(A, A0, A1, A2):- functor(A, A0, 2),
  262		arg(1, A, A1),
  263		arg(2, A, A2).
  264
  265% The most basic helpers.
  266composable_pairs(A-B, A-C, B, C).
  267composable_pairs(A-B, C-A, B, C).
  268composable_pairs(B-A, A-C, B, C).
  269composable_pairs(B-A, C-A, B, C).
  270%
  271normal_pair(A-B, B-A):- B < A, !.
  272normal_pair(A->B, B->A):- B < A, !.
  273normal_pair(X, X).
  274%
  275strong_less_than(_-A, B-_):- A<B.
  276
  277		/************************
  278		*     core predicates   *
  279		************************/
  280%
  281add_links([], X, X):-!.
  282add_links([A-Ns|Ls], X, Y):-!,
  283	cofact(X0, t(A-A, 0, X)),
  284	add_links(A, Ns, X0, X1),
  285	prune_by_frontier(A, X1, X2),
  286	add_links(Ls, X2, Y).
  287%
  288add_links(_, [], X, X):-!.
  289add_links(A, [B|Ns], X, Y):-
  290	add_link(A-B, X, X0),
  291	zdd_join(X, X0, X1),
  292	add_links(A, Ns, X1, Y).
  293%
  294add_link(_, X, 0):- X<2, !.
  295add_link(U, X, Y):-
  296	cofact(X, t(A, L, R)),
  297	arrow_symbol(A, F),
  298	(	F = (->) -> Y = 0
  299	; 	add_link(U, L, L0),
  300		(	U = A  -> R0 = 0
  301		;   strong_less_than(U, A) -> R0 = 0
  302		;  (	composable_pairs(U, A, V, W) ->
  303				U = (Ul-Ur),
  304			    subst_node([Ul->Ur], V, W, R, R0)
  305		   ;	add_link(U, R, R1),
  306			    zdd_insert(A, R1, R0)
  307		   )
  308		),
  309		zdd_join(L0, R0, Y)
  310	).
  311%
  312subst_node(_, _, _, X, 0):- X<2, !.
  313subst_node(Es, A, P, X, Y):-	cofact(X, t(U, L, R)), % replace A with P
  314	arrow_symbol(U, F, Lu, Ru),
  315	(	F = (->) ->  Y = 0
  316	;	subst_node(Es, A, P, L, L0),
  317		(	Ru = A	->
  318			normal_pair(Lu-P, V),
  319			zdd_ord_insert([V|Es], R, R0)
  320		;	Lu = A	->
  321			normal_pair(P-Ru, V),
  322			zdd_ord_insert([V|Es], R, R0)
  323		;	subst_node(Es, A, P, R, R1),
  324			zdd_insert(U, R1, R0)
  325		),
  326		zdd_join(L0, R0, Y)
  327	).
  328
  329		/***************************
  330		*     Prune by frontier    *
  331		***************************/
  332
  333prune_by_frontier(I, X, Y):-
  334	get_key(frontier, V),
  335	get_key(ends, M),
  336	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.
  343prune_by_frontier(X, X, _, _, _):- X<2, !.
  344prune_by_frontier(X, Y, I, M, V):- cofact(X, t(A, L, R)),
  345	classify_pair(A, I, M, V, C),
  346	(	C = arrow -> Y = X
  347	; 	prune_by_frontier(L, L0, I, M, V),
  348		(	C = keep ->
  349			prune_by_frontier(R, R1, I, M, V),
  350			zdd_insert(A, R1, R0)
  351		;	C = ignore ->
  352			prune_by_frontier(R, R0, I, M, V)
  353		;	R0 = 0
  354		),
  355		zdd_join(L0, R0, Y)
  356	).
  357
  358
  359% helper.
  360on_pair(J, J-_).
  361on_pair(J, _-J).
  362%
  363classify_pair((_->_), _, _, _, arrow):-!.
  364classify_pair(J-J, I, _, V, C):-!,
  365	(	on_frontier(J, I, V) -> C = keep
  366	;	C = 0
  367	).
  368classify_pair(J-K, I, Pair, V, C):-
  369	(	on_frontier(J, I, V) ->
  370		(	on_frontier(K, I, V) -> C = keep
  371		;	(	on_pair(K, Pair) -> C = keep
  372			;	C = 0
  373			)
  374		)
  375	;	on_pair(J, Pair)->
  376		(	on_frontier(K, I, V) -> C = keep
  377		;	on_pair(K, Pair) -> C = keep
  378		;	C = 0
  379		)
  380	;	C = 0
  381	).
  382
  383% ?- zdd, X<< +[*[a-b, a->b]], prune_final(a, b, X, Y), psa(X), psa(Y).
  384prune_final(_, _, X, 0):- X<2, !.
  385prune_final(P, Q, X, Y):- cofact(X, t(A, _L, R)),
  386	(	A = (_->_) -> Y = 0
  387	; 	A = P-Q    -> prune_final0(R, Y)
  388	;	Y = 0
  389	).
  390%
  391prune_final0(X, X):- X<2, !.
  392prune_final0(X, Y):- cofact(X, t(A, _L, _R)),
  393	(	A = (_->_) -> Y = X
  394	; 	Y = 0
  395	)