044: Steiner triple systems

/*
 * SICSTUS CLPFD DEMONSTRATION PROGRAM
 * Purpose   : Langford's Number Problem
 * Author    : Mats Carlsson
 *
 * | ?- steiner(dual,[],33,bound).
 */

:- module(steiner, [steiner/4]).
:- use_module(library(lists)).
:- use_module(library(clpfd)).

steiner(dual, _Opt, N, Consistency) :-
	problem(N, Triples, Consistency),
	append(Triples, Vars),
	dual_labeling(0, N, Vars),
	format_steiner(Triples, N).
steiner(byrow, Opt, N, Consistency) :-
	problem(N, Triples, Consistency),
	append(Triples, Vars),
	labeling(Opt, Vars),
	format_steiner(Triples, N).
steiner(bycol, Opt, N, Consistency) :-
	problem(N, Triples, Consistency),
	transpose(Triples, Transpose),
	append(Transpose, Vars),
	labeling(Opt, Vars),
	format_steiner(Triples, N).


format_steiner(Triples, N) :-
	format('Steiner instance of order ~d:\n', [N]),
	(   foreach(T,Triples)
	do  format('~t~d~+~t~d~+~t~d~+\n', T)
	).

dual_labeling(N, N, _) :- !.
dual_labeling(N1, N3, Vars) :-
	N2 is N1+1,
	split_by_min(Vars, N2, Cands, Rest, Rest2),
	M is N3>>1,
	dual_choose(0, M, N2, Cands, Rest2),
	dual_labeling(N2, N3, Rest).

split_by_min([], _, []) --> [].
split_by_min([X|L1], N, L2) --> [X],
	{fd_min(X, Xmin)},
	{Xmin=\=N}, !,
	split_by_min(L1, N, L2).
split_by_min([X|L1], N, [X|L2]) -->
	split_by_min(L1, N, L2).
	
dual_choose(M, M, _, Cands, Cands) :- !.
dual_choose(I, M, Val, [Val|Cands], Rest) :-
	J is I+1,
	dual_choose(J, M, Val, Cands, Rest).
dual_choose(I, M, Val, [X|Cands], [X|Rest]) :-
	X #\= Val,
	dual_choose(I, M, Val, Cands, Rest).

problem(N, Triples, Consistency) :-
	M is N mod 6,
	(M=:=1 ; M=:=3), !,
	NTrip is N*(N-1)//6,
	length(Triples, NTrip),
	(   foreach([A,B,C],Triples),
	    foreach([A,B],Tuples),
	    param(N)
	do  domain([A,B,C], 1, N),
	    A #< B, B #< C
	),
	lex_chain(Tuples,[increasing]),
	pair_constraints(Triples, N, Consistency),
	card_constraint(Triples, N).
	
pair_constraints(Triples, N, Consistency) :-
	(   foreach([A,B,C],Triples),
	    fromto(Codes,[AB,AC,BC|S],S,[]),
	    param([N,Consistency])
	do  scalar_product([N,1], [A,B], #=, AB, [consistency(Consistency)]),
	    scalar_product([N,1], [A,C], #=, AC, [consistency(Consistency)]),
	    scalar_product([N,1], [B,C], #=, BC, [consistency(Consistency)])
	),
	all_distinct(Codes, [consistency(Consistency)]).

card_constraint(Triples, N) :-
	M is N>>1,
	(   for(J,1,N),
	    foreach(J-M,Cs),
	    param(M)
	do  true
	),
	append(Triples, Vars),
	global_cardinality(Vars, Cs).