%-----------------------------------------------------------------------------
%	Huffman lattices -- suite of exploration tools
%	D. Stott Parker, June 1996.
%
%	Revised June 1997 for SICStus3;
%	also includes an improved version of the top-down Huffman algorithm.
%-----------------------------------------------------------------------------
% $Id: huffman.pl,v 1.14 1997/07/06 23:16:02 stott Exp stott $

%-----------------------------------------------------------------------------
%  Examples in the paper
%-----------------------------------------------------------------------------

examples_in_the_paper :-
	nl,
	writenl('Huffman path lengths for a given weight sequence, as well as'),
	writenl('     the space of all possible path lengths and their costs.'),
	nl,
	forAll( weight_sequence(W),
		(
			huffman_tree_path_length_sequence(W,L,V),
			write('Huffman sequence for '),writenl(W),
			write(L), write(' with cost '),writenl(V),
			novel_huffman_tree_path_length_sequence(W,Lnovel,Vnovel),
			((Vnovel = V, Lnovel = L) ->
				writenl('novel Huffman algorithm gives the same solution.')
			;
				writenl('novel Huffman algorithm gives different solution:'),
				write(Lnovel), write(' with cost '), writenl(Vnovel)
			),
			nl,
			writenl('The space of all solutions and their costs:'),
			nl,
			write_seqs_and_values(W),
			nl
		)
	),
	nl,
	writenl('GLBs and LUBs of interesting pairs of trees (of size n=9):'),
	nl,
	forAll( interesting_pair(S,T),
		(
			huffman_bounds(S,T,GLB,LUB),
			writenl( S /\ T = GLB ),
			writenl( S \/ T = LUB ),
			nl
		)
	).

weight_sequence( [10,9,8,7,3,2,1] ).
weight_sequence( [189,95,73,71,23,21,18,9,1] ).

interesting_pair( [1,2,4,5,5,5,5,5,5], [2,2,2,3,4,5,6,7,7] ).
interesting_pair( [1,2,4,5,5,5,5,5,5], [2,2,3,3,3,4,5,6,6] ).
interesting_pair( [1,3,4,4,4,4,4,5,5], [2,2,3,3,3,4,5,6,6] ).
interesting_pair( [1,4,4,4,4,4,4,4,4], [2,2,3,3,3,4,5,6,6] ).
interesting_pair( [1,4,4,4,4,4,4,4,4], [2,2,2,3,5,5,5,6,6] ).



%-----------------------------------------------------------------------------
%	Huffman tree construction
%-----------------------------------------------------------------------------

huffman_tree_path_length_sequence(W,L,V) :-
	huffman_tree(W,(_-(_,T))),
	path_length_sequence(T,L),
	code_value(L,W,V).

huffman_tree(W,T) :-
	keylist(W,L),
	keysort(L,LS),
	huffman_merge(LS,T).

keylist([],[]).
keylist([X|Xs],[Y|Ys]) :- Y=(X-(1,[X])), keylist(Xs,Ys).

huffman_merge([T],T) :- !.
huffman_merge([A,B|L],T) :-
	huffman_sum(A,B,C),
	keysort([C|L],L1),
	huffman_merge(L1,T).

huffman_sum( W1-(Level1,T1), W2-(Level2,T2), W-(Level,T) ) :-
	W is W1+W2,
	Level1plus1 is Level1+1,
	Level2plus1 is Level2+1,
	maximum(Level1plus1,Level2plus1,Level),
	T = [T1|T2].

path_length_sequence([_],[0]) :- !.
path_length_sequence([T1|T2],L) :-
	path_length_sequence(T1,L1),
	path_length_sequence(T2,L2),
	combine(L1,L2,L).

combine(L1,L2,L) :- plus1(L1,P1), plus1(L2,P2), append(P1,P2,P12), qsort(P12,L).

plus1([],[]).
plus1([X|Xs],[Y|Ys]) :- Y is X+1, plus1(Xs,Ys).


test_huffman :-
	forAll( random_weight_sequence(W), test_huffman(W) ).

test_huffman(W) :-
	huffman_tree_path_length_sequence(W,L,V),
	write('Huffman sequence for '),writenl(W),
	write(L), write(' with cost '),writenl(V),
	novel_huffman_tree_path_length_sequence(W,Lnovel,Vnovel),
	((Vnovel = V, Lnovel = L) ->
		writenl('novel Huffman algorithm gives the same solution.')
	;
		writenl('novel Huffman algorithm gives different solution:'),
		write(Lnovel), write(' with cost '), writenl(Vnovel)
	),
	nl.

random_weight_sequence( [3,2,2,1] ).
random_weight_sequence( [3,2,2,2,1,1] ).
random_weight_sequence( [4,2,2,2,1,1] ).
random_weight_sequence( [5,2,2,2,1,1] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,49] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,48] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,47] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,44] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,43] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,42] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,36] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,34] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,31] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,30] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,20] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,11] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,10] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,9] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,8] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,7] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,6] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,5] ).
random_weight_sequence( [97,93,92,85,80,79,63,59,57,56,55,50,49,1] ).
% next 20 are random 100*exponential(3).
% next 10 are random 100*chi_square(5).
% next 10 are random uniform(1..1000).
random_weight_sequence( [89,659,303,256,50,117,473,363,252,175,5,1465,99,26,51,899] ).
random_weight_sequence( [369,294,272,831,675,129,174,18,183,404,7,154,180,130,603,416] ).
random_weight_sequence( [633,53,294,160,219,286,497,471,177,542,402,374,410,29,193,498] ).
random_weight_sequence( [112,754,185,16,131,550,55,218,608,267,30,796,200,369,566,1028] ).
random_weight_sequence( [376,167,87,570,319,362,81,415,29,157,14,263,819,322,380,11] ).
random_weight_sequence( [552,7,661,289,55,750,252,463,143,97,162,684,191,90,949,75] ).
random_weight_sequence( [260,271,215,753,71,306,40,177,204,369,473,213,149,588,369,262] ).
random_weight_sequence( [155,13,265,104,521,308,308,947,15,67,242,134,923,474,438,61] ).
random_weight_sequence( [46,65,180,866,636,305,8,163,40,315,213,97,185,60,416,194] ).
random_weight_sequence( [80,623,540,196,230,47,31,614,114,74,136,100,587,163,370,243] ).
random_weight_sequence( [390,348,97,97,199,429,621,217,388,390,256,115,62,458,41,792] ).
random_weight_sequence( [378,645,382,254,28,222,22,266,407,183,313,409,103,139,96,34] ).
random_weight_sequence( [332,399,118,853,163,109,1381,283,436,349,44,206,274,526,458,218] ).
random_weight_sequence( [319,174,296,495,55,4,128,79,72,7,107,111,314,86,130,66] ).
random_weight_sequence( [26,315,173,431,37,259,136,53,444,609,66,458,1384,74,143,266] ).
random_weight_sequence( [172,204,229,111,163,183,104,1,90,103,17,153,72,79,22,176] ).
random_weight_sequence( [69,1,23,265,167,529,20,6,115,342,246,402,496,760,64,201] ).
random_weight_sequence( [256,228,168,454,83,378,90,394,44,38,93,318,316,154,89,427] ).
random_weight_sequence( [668,530,897,436,12,32,91,1223,203,387,677,438,1693,442,213,225] ).
random_weight_sequence( [194,589,131,250,402,391,1063,646,662,558,0,1178,32,30,518,107] ).
random_weight_sequence( [209,456,485,332,277,626,1266,328,128,483,425,135,355,413,97,685] ).
random_weight_sequence( [399,378,606,150,252,703,204,172,194,382,645,461,215,308,439,490] ).
random_weight_sequence( [366,203,436,153,363,267,475,431,386,730,427,383,262,206,234,643] ).
random_weight_sequence( [459,423,209,578,393,474,244,1036,116,551,522,537,241,238,244,748] ).
random_weight_sequence( [100,573,224,463,128,144,191,426,348,826,303,480,597,130,277,498] ).
random_weight_sequence( [167,541,107,886,498,327,128,311,358,871,314,658,603,1260,1054,409] ).
random_weight_sequence( [1195,360,412,464,278,326,50,161,441,552,762,623,180,109,128,531] ).
random_weight_sequence( [537,512,681,463,404,386,322,620,701,425,566,238,310,324,416,674] ).
random_weight_sequence( [195,452,109,165,598,468,296,441,563,605,508,1093,483,367,443,450] ).
random_weight_sequence( [400,107,437,638,223,428,227,178,247,528,100,145,366,117,339,147] ).
random_weight_sequence( [264,268,175,258,259,693,873,290,680,958,812,357,184,945,681,951] ).
random_weight_sequence( [374,953,479,229,540,507,510,135,177,669,221,337,871,231,469,334] ).
random_weight_sequence( [115,377,187,168,334,920,707,312,402,146,727,420,143,366,724,225] ).
random_weight_sequence( [397,265,350,837,586,150,195,528,389,160,295,115,268,500,206,119] ).
random_weight_sequence( [576,967,151,703,367,595,633,219,183,133,209,678,389,616,597,680] ).
random_weight_sequence( [807,554,935,695,867,161,783,473,215,114,207,921,844,149,780,300] ).
random_weight_sequence( [855,308,867,415,546,738,401,338,408,469,879,946,242,161,590,840] ).
random_weight_sequence( [557,868,517,906,607,473,313,940,152,568,698,416,149,514,666,797] ).
random_weight_sequence( [163,174,345,311,904,331,295,397,877,832,948,529,325,929,614,730] ).
random_weight_sequence( [693,236,570,279,443,711,263,184,557,102,352,680,972,137,440,590] ).


%-----------------------------------------------------------------------------
%	Top-down (backwards) Huffman tree construction
%-----------------------------------------------------------------------------

novel_huffman_tree_path_length_sequence(W,L,V) :-
	qsort(W,Wascending),
	optimal_tree_path_length_rev_sequence(Wascending,Lr),
	reverse(Lr,L),
	code_value(L,W,V).

optimal_tree_path_length_rev_sequence([_],[0]) :- !.
optimal_tree_path_length_rev_sequence([Wn,Wn1|Ws],OptL) :-
	W0 is Wn+Wn1,
	insertion(W0,Ws,Wascending),
	optimal_tree_path_length_rev_sequence(Wascending, ContractionOptL),
	(   (
	    rev_suffix_increment(ContractionOptL,J),
	    rev_suffix_length(ContractionOptL,TwoK),
	    RemnantLength is TwoK-2,
	    rev_suffix_split(RemnantLength,Ws,Wn_2k,Wn_2k_1),
	    Wn_2k_1 - (J-1)*Wn_2k  =<  Wn1 + Wn
	    ) ->
		rev_expansion(lower,ContractionOptL,OptL)
	;
		rev_expansion(upper,ContractionOptL,OptL)
	).


insertion(X,[],[X]) :- !.
insertion(X,[Y|Ys],[X,Y|Ys]) :- X<Y, !.
insertion(X,[Y|Ys],[Y|XYs])  :- X>=Y, insertion(X,Ys,XYs).


%-----------------------------------------------------------------------------
%	Majorization lattice glb and lub
%-----------------------------------------------------------------------------

majorization_bounds(S,T,GLB,LUB) :-
	minuslog2seq(X,S),
	minuslog2seq(Y,T),
	running_sum(X,X_Z),
	running_sum(Y,Y_Z),
	element_wise_minimum(X_Z,Y_Z,GLB_Z),
	running_sum(GLB_Z_M,GLB_Z),
	minuslog2seq(GLB_Z_M,GLB),
	element_wise_maximum(X_Z,Y_Z,LUB_Z),
	running_sum(LUB_Z_M,LUB_Z),
	minuslog2seq(LUB_Z_M,LUB).


%-----------------------------------------------------------------------------
%	Huffman lattice glb and lub (via the recursive algorithm in the paper)
%-----------------------------------------------------------------------------

huffman_bounds(S,T,GLB,LUB) :-
	reverse(S,SR),
	reverse(T,TR),
	huffman_glb_lub(SR,TR,GLBR,LUBR),
	reverse(GLBR,GLB),
	reverse(LUBR,LUB).

huffman_glb_lub(X,Y,X,Y) :- pow2minus_rev_majorizes(Y,X), !.
huffman_glb_lub(X,Y,Y,X) :- pow2minus_rev_majorizes(X,Y), !.
huffman_glb_lub([X,X|Xs],[Y,Y|Ys],GLB,LUB) :-
	X1 is X-1,
	Y1 is Y-1,
	insert(X1,Xs,Xc),
	insert(Y1,Ys,Yc),
	huffman_glb_lub(Xc,Yc,GLB1,LUB1),
	extend_glb(GLB1,[X,X|Xs],[Y,Y|Ys],GLB),
	extend_lub(LUB1,[X,X|Xs],[Y,Y|Ys],LUB).

extend_glb(GLB1,X,Y,GLB) :-
	setof(A, member(A,GLB1), P), % all possibilities in ascending order
	reverse(P,Possibilities),
	member(A,Possibilities), % consider extensions in descending order
		remove_first(A,GLB1,GLB2),
		A1 is A+1,
		insert(A1,GLB2,GLB3),
		insert(A1,GLB3,GLB),
	pow2minus_rev_majorizes(X,GLB),
	pow2minus_rev_majorizes(Y,GLB),
	!.	% both X and Y dominate GLB, as desired

extend_lub(LUB1,X,Y,LUB) :-
	setof(A, member(A,LUB1), P), % all possibilities in ascending order
	member(A,P), % consider extensions in ascending order
		remove_first(A,LUB1,LUB2),
		A1 is A+1,
		insert(A1,LUB2,LUB3),
		insert(A1,LUB3,LUB),
	pow2minus_rev_majorizes(LUB,X),
	pow2minus_rev_majorizes(LUB,Y),
	!.	% LUB dominates both X and Y, as desired

insert(X,[],[X]) :- !.
insert(X,[Y|Ys],[X,Y|Ys]) :- X>=Y, !.
insert(X,[Y|Ys],[Y|XYs])  :- X<Y, insert(X,Ys,XYs).

remove_first(X,[X|Xs],Xs) :- !.
remove_first(X,[Y|Ys],[Y|XYs]) :- remove_first(X,Ys,XYs).

pow2minus_rev_majorizes(X,Y) :-
	minuslog2seq(PX,X),
	minuslog2seq(PY,Y),
	rev_majorizes(PX,PY,_,_).

rev_majorizes([],[], 0.0, 0.0).
rev_majorizes([X],[Y],X,Y) :- !, X>=Y.
rev_majorizes([X|Xs],[Y|Ys],XXsum,YYsum) :-
	rev_majorizes(Xs,Ys,Xsum,Ysum),
	XXsum is Xsum+X,
	YYsum is Ysum+Y,
	XXsum >= YYsum.
	

minuslog2seq([],[]).
minuslog2seq([X|Xs],[Y|Ys]) :- nonvar(X), !, log2(1.0/X,Y), minuslog2seq(Xs,Ys).
minuslog2seq([X|Xs],[Y|Ys]) :- nonvar(Y), !, pow2(-Y,X), minuslog2seq(Xs,Ys).


running_sum(X,T) :- running_sum(X,T,0).

running_sum([],[],_).
running_sum([X|Xs],[T|Ts],T0) :- nonvar(X), !, T is X+T0, running_sum(Xs,Ts,T).
running_sum([X|Xs],[T|Ts],T0) :- nonvar(T), !, X is T-T0, running_sum(Xs,Ts,T).


element_wise_minimum([],[],[]). 
element_wise_minimum([X|Xs],[Y|Ys],[XY|XYs]) :-
	minimum(X,Y,XY),
	element_wise_minimum(Xs,Ys,XYs). 

element_wise_maximum([],[],[]). 
element_wise_maximum([X|Xs],[Y|Ys],[XY|XYs]) :-
	maximum(X,Y,XY),
	element_wise_maximum(Xs,Ys,XYs). 

element_wise_leq([],[]).
element_wise_leq([X|Xs],[Y|Ys]) :- X =< Y, element_wise_leq(Xs,Ys).


%-----------------------------------------------------------------------------
%	Huffman lattice glb and lub by lower and upper expansion analysis
%-----------------------------------------------------------------------------

alternative_huffman_bounds(S,T,GLB,LUB) :-
	reverse(S,SR),
	reverse(T,TR),
	alternative_huffman_glb_lub(SR,TR,GLBR,LUBR),
	reverse(GLBR,GLB),
	reverse(LUBR,LUB),
	!.
alternative_huffman_bounds(S,T,_,_) :-
	write('alternative_huffman_bounds failure:  S,T = '),write((S,T)),nl.

alternative_huffman_glb_lub(X,Y,X,Y) :- pow2minus_rev_majorizes(Y,X), !.
alternative_huffman_glb_lub(X,Y,Y,X) :- pow2minus_rev_majorizes(X,Y), !.
alternative_huffman_glb_lub(X,Y,GLB,LUB) :-
	rev_contraction(X,Xc,Xtype),
	rev_contraction(Y,Yc,Ytype),
	alternative_huffman_glb_lub(Xc,Yc,GLBc,LUBc),
	(Xtype = Ytype ->
		rev_expansion(Xtype, GLBc, GLB),
		rev_expansion(Xtype, LUBc, LUB)
	;
		alternative_extend_glb(GLBc,X,Y,GLB),
		alternative_extend_lub(LUBc,X,Y,LUB)
	).

rev_contraction([Q,Q,Q_j|Qs],Contraction,Type) :-
	Q_1 is Q-1,
	(Q_1 >= Q_j ->
		Type = upper,
		Contraction=[Q_1,Q_j|Qs]
	;
		Type = lower,
		insert(Q_1,[Q_j|Qs],Contraction)
	).

rev_expansion( upper, [Q|Qs], [Q1,Q1|Qs] ) :- !, Q1 is Q+1.
rev_expansion( lower, [Q|Qs], [Q|Ys] ) :- expand_after(Q,Qs,Ys), !.
rev_expansion( lower, [Q|Qs], [Q1,Q1|Qs] ) :- !, Q1 is Q+1.
% use upper expansion in the rare situation that lower is impossible.

expand_after(Q,[Q|Xs],[Q|Ys]) :- !, expand_after(Q,Xs,Ys).
expand_after(_,[X|Xs],[X1,X1|Xs]) :- X1 is X+1.

rev_suffix_length([Q|Qs],N1) :- rev_suffix_length(Q,Qs,N), N1 is N+1.
rev_suffix_length(Q,[Q|Xs],N1) :- !, rev_suffix_length(Q,Xs,N), N1 is N+1.
rev_suffix_length(_,_,0).

rev_suffix_increment([Q|Qs],J) :- rev_suffix_increment(Q,Qs,J).
rev_suffix_increment(Q,[Q|Xs],J) :- !, rev_suffix_increment(Q,Xs,J).
rev_suffix_increment(Q,[Q_J|_],J) :- !, J is Q-Q_J.
rev_suffix_increment(_,[],0).

rev_suffix_split(0,[Wn_2k,Wn_2k_1|_],Wn_2k,Wn_2k_1) :- !.
rev_suffix_split(I,[_|Ws],Wn_2k,Wn_2k_1) :-
	I1 is I-1,
	rev_suffix_split(I1,Ws,Wn_2k,Wn_2k_1).

alternative_extend_glb(GLB1,X,Y,GLB) :-
	setof(A, member(A,GLB1), P), % all possibilities in ascending order
	reverse(P,[V,V1|_]),
	member(A,[V,V1]), % consider extensions in descending order
		remove_first(A,GLB1,GLB2),
		A1 is A+1,
		insert(A1,GLB2,GLB3),
		insert(A1,GLB3,GLB),
	pow2minus_rev_majorizes(X,GLB),
	pow2minus_rev_majorizes(Y,GLB),
	!.	% both X and Y dominate GLB, as desired

alternative_extend_lub(LUB1,X,Y,LUB) :-
	setof(A, member(A,LUB1), P), % all possibilities in ascending order
	reverse(P,[V,V1|_]),
	member(A,[V1,V]), % consider extensions in ascending order
		remove_first(A,LUB1,LUB2),
		A1 is A+1,
		insert(A1,LUB2,LUB3),
		insert(A1,LUB3,LUB),
	pow2minus_rev_majorizes(LUB,X),
	pow2minus_rev_majorizes(LUB,Y),
	!.	% LUB dominates both X and Y, as desired

verify_alternative(N) :-
	forAll( tree_path_length_sequence(N,S),
		forAll( tree_path_length_sequence(N,T),
			verify_alternative(S,T) ) ).

verify_alternative(S,T) :-
	huffman_bounds(S,T,GLB,LUB),
	alternative_huffman_bounds(S,T,GLBr,LUBr),
	(GLB==GLBr ->
		true
	;
		write('for S,T = '),write((S,T)),write(' GLB = '),write(GLB),
		write('but alternative GLB = '),write(GLBr),nl
	),
	(LUB==LUBr ->
		true
	;
		write('for S,T = '),write((S,T)),write(' LUB = '),write(LUB),
		write('but alternative LUB = '),write(LUBr),nl
	).


%-----------------------------------------------------------------------------
%	Huffman bounds by truncation (alas, this approach does not always work)
%-----------------------------------------------------------------------------

heuristic_huffman_bounds(S,T,GLB,LUB) :-
	minuslog2seq(X,S),
	minuslog2seq(Y,T),
	running_sum(X,X_Z),
	running_sum(Y,Y_Z),
	running_minorizer(X_Z,Y_Z,0.0,GLB),
	running_majorizer(X_Z,Y_Z,0.0,LUB).

running_minorizer([],[],_,[]).
running_minorizer([X|Xs],[Y|Ys],Sum,[L|Ls]) :-
	minimum(X-Sum,Y-Sum,Bound),
	log2(Bound,LBound),
	ceiling_value(-LBound,L),
	pow2(-L,PL),
	SumPL is Sum+PL,
	running_minorizer(Xs,Ys,SumPL,Ls).

running_majorizer([],[],_,[]).
running_majorizer([X|Xs],[Y|Ys],Sum,[L|Ls]) :-
	maximum(X-Sum,Y-Sum,Bound),
	log2(Bound,LBound),
	floor_value(-LBound,L),
	pow2(-L,PL),
	SumPL is Sum+PL,
	running_majorizer(Xs,Ys,SumPL,Ls).

huffman_truncation(W,L) :-
	normalization(W,Wn),
	minuslog2seq(Wn,Ln),
	higher_truncation(Ln,L).

lower_truncation([],[]).
lower_truncation([V|Vs],[T|Ts]) :-
	ceiling_value(V,T),
	pow2(-V,PV),
	pow2(-T,PT),
	lower_truncation(Vs,PV,PT,Ts).

lower_truncation([],_,_,[]).
lower_truncation([V|Vs],PVsum,PTsum,[T|Ts]) :-
	pow2(-V,PV),
	PVsumPV is PVsum+PV,
	log2(PVsumPV-PTsum,LT),
	ceiling_value(-LT,T),
	pow2(-T,PT),
	PTsumPT is PTsum+PT,
	lower_truncation(Vs,PVsumPV,PTsumPT, Ts).

higher_truncation(X,Y) :-
	reverse(X,XR),
	lower_truncation(XR,YR),
	reverse(YR,Y).


%-----------------------------------------------------------------------------
%	Exhaustive testing of submodularity (using a fixed weight sequence W).
%-----------------------------------------------------------------------------

verify_submodularity(T) :-
	given_weight_sequence(T,W),
	verify_submodularity(W).

given_weight_sequence( 1, [10,9,8,7,3,2,1] ).
given_weight_sequence( 2, [16,16,16,8,8,4,2,2,1] ).
given_weight_sequence( 3, [128,64,32,32,16,16,16,4,1] ).
given_weight_sequence( 4, [189,95,73,71,23,21,18,9,1] ).
given_weight_sequence( 5, [512,512,128,128,128,64,64,32,16,16,16,8,8,4,2,2,1] ).


verify_submodularity(W) :-
	length(W,N),
	forAll( tree_path_length_sequence(N,S),
		forAll( tree_path_length_sequence(N,T),
			verify_submodularity(W,S,T) ) ).

verify_submodularity(W,S,T) :-
	huffman_bounds(S,T,GLB,LUB),
	code_value(S,W,S_value),
	code_value(T,W,T_value),
	code_value(GLB,W,GLB_value),
	code_value(LUB,W,LUB_value),
	Submodularity is (S_value+T_value) - (GLB_value+LUB_value),
	(Submodularity < 0 ->
	    (write('non-submodular '),write((S,T,GLB,LUB)),
		write(' value = '),write(Submodularity),nl) ; true).

%-----------------------------------------------------------------------------
%	Verify isomorphism between the imbalance and majorization orderings
%	(when restricted to path-length sequences):
%	2^{-s} is majorized by 2^{-t}   ==   s is at least as balanced as t.
%-----------------------------------------------------------------------------

verify_imbalance_majorization_equivalence(N) :-
	forAll( tree_path_length_sequence(N,S),
		forAll( tree_path_length_sequence(N,T),
			verify_imbalance_majorization_equivalence(S,T) ) ).


verify_imbalance_majorization_equivalence(S,S) :- !.
verify_imbalance_majorization_equivalence(S,T) :-
	(pow2minus_majorized_by(S,T) ->
		(at_least_as_balanced_as(S,T) ->  
			true
		;
			write((S < T)),
			write(' wrt majorization, but not wrt imbalance'),
			nl
		)
	;
		(at_least_as_balanced_as(S,T) -> 
			write((S < T)),
			write(' wrt imbalance, but not wrt majorization'),
			nl
		;
			true
		)
	).

pow2minus_majorized_by(S,T) :-
	minuslog2seq(X,S),
	minuslog2seq(Y,T),
	running_sum(X,X_Z),
	running_sum(Y,Y_Z),
	element_wise_leq(X_Z,Y_Z).

at_least_as_balanced_as(S,T) :-
	total(S,Stotal),
	total(T,Ttotal),
	ImbalanceLevelDifference is Ttotal-Stotal,
	ImbalanceLevelDifference > 0,
	at_least_as_balanced_as(S,T,ImbalanceLevelDifference),
	!.

at_least_as_balanced_as(S,S,_) :- !.
at_least_as_balanced_as(S,T,Level) :-
	Level > 0,
	succ(Level1,Level),
	balancing_exchange(T1,T),
	at_least_as_balanced_as(S,T1,Level1).

balancing_exchange([X|Xs],[X|Ys])  :- balancing_exchange(Xs,Ys).
balancing_exchange([P1,P1|Xs],[P|Ys]) :- succ(P,P1), remaining_exchange(Xs,Ys).

remaining_exchange([X|Xs],[X|Ys]) :- remaining_exchange(Xs,Ys).
remaining_exchange([Q|Xs],[Q1,Q1|Xs]) :- succ(Q,Q1).


%-----------------------------------------------------------------------------
%	Find all majorization GLBs and LUBs that are not purely integral
%-----------------------------------------------------------------------------

verify_bound_truncation(N) :-
	forAll( tree_path_length_sequence(N,S),
		forAll( tree_path_length_sequence(N,T),
			verify_bound_truncation(S,T) ) ).


verify_bound_truncation(S,S) :- !.
verify_bound_truncation(S,T) :-
	majorization_bounds(S,T,MGLB,MLUB),
	(non_integral_seq(MGLB) ->
		((lower_truncation(MGLB,GLB),
		  at_least_as_balanced_as(GLB,S),
		  at_least_as_balanced_as(GLB,T)) ->
			true
		;
			write((S/\T)=GLB),
			write(' computed, but not verified'),
			nl
		)
	;
		true
	),
	(non_integral_seq(MLUB) ->
		((higher_truncation(MLUB,LUB),
		  at_least_as_balanced_as(S,LUB),
		  at_least_as_balanced_as(T,LUB)) ->
			true
		;
			write((S\/T)=LUB),
			write(' computed, but not verified'),
			nl
		)
	;
		true
	).

find_non_integral_bounds(N) :-
	forAll( tree_path_length_sequence(N,S),
		forAll( tree_path_length_sequence(N,T),
			find_non_integral_bounds(S,T) ) ).


find_non_integral_bounds(S,S) :- !.
find_non_integral_bounds(S,T) :-
	majorization_bounds(S,T,MGLB,MLUB),
	(non_integral_seq(MGLB) ->
	    (write(S),write(' /\\ '),write(T),write(' = '),write(MGLB),
	     lower_truncation(MGLB,GLB),
	     write(' --- tree glb: '),write(GLB),nl); true),
	(non_integral_seq(MLUB) ->
	    (write(S),write(' \\/ '),write(T),write(' = '),write(MLUB),
	     higher_truncation(MLUB,LUB),
	     write(' --- tree lub: '),write(LUB),nl); true).

%-----------------------------------------------------------------------------
%	Printing of all tree path-length sequences, along with some property.
%-----------------------------------------------------------------------------

write_seqs_and_seq_entropies(N) :-
	writenl('sequence -> entropy value'),
	forAll( tree_path_length_sequence(N,L),
		(entropy_value(L,H),writenl((L->H))) ).

write_seqs_and_values(W) :-
	length(W,N),
	writenl('sequence -> code value'),
	forAll( tree_path_length_sequence(N,L),
		(code_value(L,W,V),writenl((L->V))) ).

write_seqs_and_wpls(W) :-
	length(W,N),
	normalization(W, Wn),
	minuslog2seq(Wn,MinusLogW),
	entropy_value(MinusLogW,Entropy),
	write('entropy lower bound for W ='), writenl(Entropy),
	% write('lower bound for L ='), writenl(MinusLogW),
	writenl('sequence -> normalized weighted path length'),
	forAll( tree_path_length_sequence(N,L),
		(code_value(L,Wn,Vn),
		 % entropy_check(Entropy,L),
		 % weak_majorization_check(MinusLogW,L),
		 writenl((L->Vn))) ).

normalization(Ws, Xs) :- total(Ws,S), normalize_by(Ws,S,Xs).

total([],0).
total([X|Xs],T) :- total(Xs,S), T is X+S.

normalize_by([],_,[]).
normalize_by([W|Ws],D,[X|Xs]) :- X is W/D, normalize_by(Ws,D,Xs).


code_value([],[],0).
code_value([L1|Ls],[W1|Ws],V) :-
	code_value(Ls,Ws,Vs),
	V is L1*W1 + Vs.

entropy_value([],0).
entropy_value([L|Ls],H) :- entropy_value(Ls,T), pow2(-L,PL), H is T+L*PL.


entropy_check(Entropy,L) :-
	entropy_value(L,H),
	H >= Entropy -> write('above entropy: ') ; write('below entropy: ').
	

weak_majorization_check(X,Y) :-
	weakly_majorized_by(X,Y) -> write('yes: ') ; write(' no: ').

weakly_majorized_by(X,Y) :- weakly_majorized_by(X,Y,0,0).

weakly_majorized_by([],[],_,_).
weakly_majorized_by([X|Xs],[Y|Ys],XT,YT) :-
	NXT is X+XT, NYT is Y+YT, NXT =< NYT,
	weakly_majorized_by(Xs,Ys,NXT,NYT).


%-----------------------------------------------------------------------------
%  Generation of tree path-length sequences
%
%  L is a tree sequence for a binary tree t with N leaves
%  if t can be partitioned into two immediate subtrees t1 and t2
%  such that:
%	t1 has N1 leaves,
%	t2 has N2 leaves,
%	the tree sequence for t1 is L1,
%	the tree sequence for t2 is L2,
%	every value in L1 is =< every value in L2,
%	L is essentially L1 appended with L2, but every element in this
%		list is incremented by one to reflect the increase in
%		depth caused by combining the two trees into a larger tree.
%-----------------------------------------------------------------------------

write_seqs(N) :-
	forAll( tree_path_length_sequence(N,L), writenl(L) ).

tree_path_length_sequence(1,[0]).
tree_path_length_sequence(2,[1,1]).
tree_path_length_sequence(3,[1,2,2]).
tree_path_length_sequence(N,L) :- 
	N >= 4,
	partition(N,N1,N2),
	tree_path_length_sequence(N1,L1),
	last_member(L1,X),
	tree_path_length_sequence(N2,L2),
	first_member(L2,Y),
	X =< Y,
	deepen_one_level(L1,L1d),
	deepen_one_level(L2,L2d),
	append(L1d,L2d,L).

deepen_one_level([],[]).
deepen_one_level([N|L],[N1|Ld]) :- N1 is N+1, deepen_one_level(L,Ld).

partition(N,N1,N2) :-
	peano(S,N),
	plus(S1,S2,S),
	nonzero(S1),
	nonzero(S2),
	peano(S1,N1),
	peano(S2,N2).

% we use Peano arithmetic so we can backtrack easily through all solutions

plus(0,X,X).
plus(s(X),Y,s(Z)) :- plus(X,Y,Z).

nonzero(s(_)).

peano(0,0) :- !.
peano(s(X),N1) :- nonvar(N1), N1>=0, !, N is N1-1, peano(X,N).
peano(s(X),N1) :- peano(X,N), N1 is N+1.

%-----------------------------------------------------------------------------
%	Procedural utilities
%-----------------------------------------------------------------------------

forAll(G,T) :- \+ (G, \+ T).

writenl(T) :- write(T), nl.

%-----------------------------------------------------------------------------
%	List utilities
%-----------------------------------------------------------------------------

qsort(L,R) :- qsort(L,R,[]).

qsort([],R,R).
qsort([X|L],R,R0) :-
	partition(L,X,L1,L2),
	qsort(L2,R1,R0),
	qsort(L1,R,[X|R1]).

partition([],_,[],[]).
partition([X|L],Y,L1,[X|L2]) :- X > Y, !, partition(L,Y,L1,L2).
partition([X|L],Y,[X|L1],L2) :- X =< Y, partition(L,Y,L1,L2).

last_member([X],X) :- !.
last_member([_|L],X) :- last_member(L,X).

first_member([X|_],X).

member(X,[X|_]).
member(X,[_|L]) :- member(X,L).

append([],L,L).
append([X|L1],L2,[X|L3]) :- append(L1,L2,L3).

reverse(L,R) :- reverse(L,[],R).

reverse([],R,R).
reverse([X|Xs],S,R) :- reverse(Xs,[X|S],R).

%-----------------------------------------------------------------------------
%	Numeric computation utilities
%-----------------------------------------------------------------------------

minimum(X,Y,X) :- X<Y, !.
minimum(_,Y,Y).

maximum(X,Y,X) :- X>Y, !.
maximum(_,Y,Y).

succ(X,X1) :- number(X), !, X1 is X+1.
succ(X,X1) :- number(X1), X is X1-1.

log2(X,Y) :-
	FloatX is 1.0*X,
	log(FloatX,LogX),
	log(2.0,Log2),
	Log2X is LogX/Log2,
	(integral(Log2X) -> integer_value(Log2X,Y) ; Y = Log2X).

pow2(X,Y) :-
	FloatX is 1.0*X,
	pow(2.0, FloatX, TwoX),
	(integral(TwoX) -> integer_value(TwoX,Y) ; Y = TwoX).


ceiling_value(X,N) :- floor_value(-X,MFloorX), N is -MFloorX.

floor_value(X,N) :- integral(X), !, integer_value(X,N).
floor_value(X,N) :- F is 1.0*X, floor(F,Y), integer_value(Y,N).

integer_value(X,N) :- X<0, !, integer_value(-X,MN), N is -MN.
integer_value(X,0) :- essentially_zero(X), !.
integer_value(X,N) :- X>0, X1 is X-1, integer_value(X1,N1), N is N1+1.

% hacks are used here to circumvent Prolog's conversion of ints to floats.

integral(X) :- integer(X), !.
integral(X) :- F is 1.0*X, floor(F,FloorX), D is X-FloorX, essentially_zero(D).

essentially_zero(X) :- X =:= 0, !.
essentially_zero(X) :- -0.0000000001 < X, X <  0.0000000001.

non_integral_seq([X|_]) :- \+ integral(X), !.
non_integral_seq([_|Xs]) :- non_integral_seq(Xs).

%-----------------------------------------------------------------------------
%	Foreign function interface to standard C math library
%
%	Note all inputs to these functions must be floats;
%	failure results if they are ints.
%-----------------------------------------------------------------------------

foreign( sqrt,	sqrt(+float,[-float])	).
foreign( log1p,	log1p(+float,[-float])	).  % log(1+X)  [for X near zero]
foreign( log10,	log10(+float,[-float])	).
foreign( log,	log(+float,[-float])	).
% foreign( exp,	exp(+float,[-float])	).
foreign( expm1,	expm1(+float,[-float])	).  % exp(X)-1  [for X near zero]
foreign( sin,	sin(+float,[-float])	).
foreign( cos,	cos(+float,[-float])	).
foreign( tan,	tan(+float,[-float])	).
foreign( asin,	asin(+float,[-float])	).
foreign( acos,	acos(+float,[-float])	).
foreign( atan,	atan(+float,[-float])	).
foreign( sinh,	sinh(+float,[-float])	).
foreign( cosh,	cosh(+float,[-float])	).
foreign( tanh,	tanh(+float,[-float])	).
foreign( floor,	floor(+float,[-float])	).
foreign( pow,	pow(+float,+float,[-float])	).

foreign_file('math.o',
	[ sqrt, log1p, log10, log, expm1,
	sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, floor, pow ]).

:- load_foreign_files(['math.o'],['-lm']).