% Copyright (C) 2005 Markus Triska triska@gmx.at

:- use_module(library(simplex)).

%:- include('v150.txt').
%:- include('v25.txt').
:- include('v15.txt').
%:- include('v3.txt').
%:- include('v21.txt').

max_weight(100).
max_value(200).

all_zero([]).
all_zero([0|As]) :-
	all_zero(As).

set_nth0(0, [_|Rest], Element, [Element|Rest]) :- !.
set_nth0(N, [E|Es], Element, [E|Rest]) :-
	N1 is N - 1,
	set_nth0(N1, Es, Element, Rest).

init_patterns([], [], _, _, []).
init_patterns([W|Ws], [V|Vs], CurrItem, ZeroPattern, [Pattern|Ps]) :-
	max_weight(MaxWeight),
	max_value(MaxValue),
	Max is min(floor(MaxWeight / W), floor(MaxValue  / V)),	
	set_nth0(CurrItem, ZeroPattern, Max, Pattern),
	format("~w\n", [Pattern]),
	CurrItem1 is CurrItem + 1,
	init_patterns(Ws, Vs, CurrItem1, ZeroPattern, Ps).

init_constraints([], [], _, S, S).
init_constraints([D|Ds], [P|Ps], N0, S0, S) :-
	nth0(N0, P, P_j_j),
	Max is ceil(D / P_j_j),
	constraint([use(N0)] =< Max, S0, S1),
	N1 is N0 + 1,
	init_constraints(Ds, Ps, N1, S1, S).

init_objective([], _, []).
init_objective([_|Rest], N0, [use(N0)|Us]) :-
	N1 is N0 + 1,
	init_objective(Rest, N1, Us).

init_demands([], [], _, S, S).
init_demands([D|Ds], [P|Ps], N0, S0, S) :-
	make_sum(P, 0, Sum),
	constraint(demand(N0), Sum >= D, S0, S1),
	N1 is N0 + 1,
	init_demands(Ds, Ps, N1, S1, S).

make_sum([], _, []).
make_sum([Coeff|Cs], N0, Rest0) :-
	( Coeff =:= 0 ->
		Rest0 = Rest
	;
		Rest0 = [Coeff*use(N0)|Rest]
	),
	N1 is N0 + 1,
	make_sum(Cs, N1, Rest).

camels :-
	weights(Weights),
	values(Values),
	demands(Demands),
	length(Demands, NumItems),
	length(ZeroPattern, NumItems),
	all_zero(ZeroPattern),
	init_patterns(Weights, Values, 0, ZeroPattern, Patterns0),
	gen_state(S0),
	init_constraints(Demands, Patterns0, 0, S0, S1),
	init_objective(ZeroPattern, 0, Obj0),
	init_demands(Demands, Patterns0, 0, S1, S2),
	camels(Demands, S2, Obj0, NumItems, S3, Obj1, NumPatts),
	format("Vector:\n"),
	make_integrals(0, NumPatts, S3, S4),
	minimize(Obj1, S4, Solved),
	print_use(0, NumPatts, Solved).

print_use(N, N, _) :- !,
	nl.
print_use(K, N, Solved) :-
	variable_value(Solved, use(K), Val),
	format("~w ", [Val]),
	K1 is K + 1,
	print_use(K1, N, Solved).
	

make_integrals(N, N, S, S) :- !.
make_integrals(N, M, S0, S) :-
	constraint(integral(use(N)), S0, S1),
	N1 is N + 1,
	make_integrals(N1, M, S1, S).

camels(Demands, State0, Obj0, NPatt0, State, Obj, NPatt) :-
	knapsack(State0, Obj0, ZBest, NewPattern),
	( ZBest =:= 1 ->
		% format("no profitable column found\n"),
		NPatt = NPatt0,
		State = State0,
		Obj = Obj0
	;
		format("~w\n", [NewPattern]),
		UseVar = use(NPatt0),
		adjust_constraints(Demands, UseVar, NewPattern, 0, 0, Upper, State0, State1),
		constraint([UseVar] =< Upper, State1, State2),
		NPatt1 is NPatt0 + 1,
		camels(Demands, State2, [UseVar|Obj0], NPatt1, State, Obj, NPatt)
	).

adjust_constraints([], _, [], _, Upper, Upper, State, State).
adjust_constraints([Demand|Ds], UV, [P|Ps], Curr, Upper0, Upper, State0, State) :-
	Curr1 is Curr + 1,
	( P =:= 0 ->
		adjust_constraints(Ds, UV, Ps, Curr1, Upper0, Upper, State0, State)
	;
		constraint_add(demand(Curr), [P*UV], State0, State1),
		Upper1 is max(Upper0, ceil(Demand/P)),
		adjust_constraints(Ds, UV, Ps, Curr1, Upper1, Upper, State1, State)
	).



knapctrl([], _, []).
knapctrl([W|Ws], N0, [W*x(N0)|Rest]) :-
	N1 is N0 + 1,
	knapctrl(Ws, N1, Rest).

knap_objective([], _, _, []).
knap_objective([_|Rest], Solved, N0, [Coeff*x(N0)|Ss]) :-
	shadow_price(Solved, demand(N0), Coeff),
	N1 is N0 + 1,
	knap_objective(Rest, Solved, N1, Ss).

knap_integers([], _, S, S).
knap_integers([_|Rest], N0, S0, S) :-
	constraint(integral(x(N0)), S0, S1),
	N1 is N0 + 1,
	knap_integers(Rest, N1, S1, S).
	

knapsack(Relaxed, RelaxedObj, ObjValue, Pattern) :-
	minimize(RelaxedObj, Relaxed, Solved),
	gen_state(KnapS0),
	weights(Ws),
	knapctrl(Ws, 0, KnapSum1),
	max_weight(MaxWeight),
	constraint(KnapSum1 =< MaxWeight, KnapS0, KnapS1),
	values(Vs),
	knapctrl(Vs, 0, KnapSum2),
	max_value(MaxValue),
	constraint(KnapSum2 =< MaxValue, KnapS1, KnapS2),
	knap_integers(Ws, 0, KnapS2, KnapS3),
	knap_objective(Ws, Solved, 0, KnapObj),
	maximize(KnapObj, KnapS3, KnapSolved),
	objective(KnapSolved, ObjValue),
	knap_values(Ws, KnapSolved, 0, Pattern).

knap_values([], _, _, []).
knap_values([_|As], Solved, N0, [V|Vs]) :-
	variable_value(Solved, x(N0), V),
	N1 is N0 + 1,
	knap_values(As, Solved, N1, Vs).