/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   Solving a bin packing task by delayed column generation.
   Written 2005, 2015, 2022 by Markus Triska (triska@metalevel.at).
   Public domain code. Tested with Scryer Prolog.

   A data file needs to provide weights/1, values/1 and demands/1.

   For example, with data.pl consisting of:

      weights([30,20,50]).
      values([90,70,20]).
      demands([12,7,15]).

   We can run:

      $ scryer-prolog camels.pl data.pl

   and obtain a solution with:

      ?- camels.
      % initial patterns:
      [2,0,0].
      [0,2,0].
      [0,0,2].
      % generated packing patterns:
      [0,2,1].
      [1,1,1].
      % solution vector (coefficient of each pattern):
      [0,0,1,0,13].
         true.

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

:- use_module(library(simplex)).
:- use_module(library(clpz)).
:- use_module(library(format)).
:- use_module(library(dcgs)).
:- use_module(library(lists)).
:- use_module(library(between)).

max_weight(100).
max_value(200).

set_nth0(N, Element, Es0, Es) :-
        length(Prefix, N),
        append(Prefix, [_|Rest], Es0),
        append(Prefix, [Element|Rest], Es).

initial_patterns([], _, _, _) --> [].
initial_patterns([W|Ws], [V|Vs], Item, ZeroPattern) -->
        [Pattern],
        { max_weight(MaxWeight),
          max_value(MaxValue),
          Max #= min(MaxWeight // W, MaxValue // V),
          set_nth0(Item, Max, ZeroPattern, Pattern),
          portray_clause(Pattern),
          Item1 #= Item + 1 },
        initial_patterns(Ws, Vs, Item1, ZeroPattern).

init_constraints(D, P, N0-S0, N-S) :-
        nth0(N0, P, P_j_j),
        Max is ceiling(D / P_j_j),
        constraint([use(N0)] =< Max, S0, S),
        N #= N0 + 1.

integer_use(N, use(N)).

init_demands(D, P, N0-S0, N-S) :-
        phrase(make_sum(P, 0), Sum),
        constraint(demand(N0), Sum >= D, S0, S),
        N #= N0 + 1.

make_sum([], _) --> [].
make_sum([Coeff|Cs], N0) -->
        (   { Coeff =:= 0 } -> []
        ;   [Coeff*use(N0)]
        ),
        { N1 #= N0 + 1 },
        make_sum(Cs, N1).

camels :-
        weights(Weights),
        values(Values),
        demands(Demands),
        length(Demands, NumItems),
        length(ZeroPattern, NumItems),
        maplist(=(0), ZeroPattern),
        format("% initial patterns:~n", []),
        phrase(initial_patterns(Weights, Values, 0, ZeroPattern), Patterns0),
        gen_state(S0),
        foldl(init_constraints, Demands, Patterns0, 0-S0, _-S1),
        MaxUse #= NumItems - 1,
        numlist(0, MaxUse, UseNs),
        maplist(integer_use, UseNs, Obj0),
        foldl(init_demands, Demands, Patterns0, 0-S1, _-S2),
        append(Obj0, NewPatterns, UseVars),
        format("% generated packing patterns:~n", []),
        phrase(camels(Demands, S2, Obj0, NumItems, S3, Obj1), NewPatterns),
        format("% solution vector (coefficient of each pattern):~n", []),
        foldl(make_integral, UseVars, S3, S4),
        minimize(Obj1, S4, Solved),
        maplist(variable_value(Solved), UseVars, Vals),
        portray_clause(Vals).

make_integral(UseVar, S0, S) :- constraint(integral(UseVar), S0, S).

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

adjust_constraints([], _, [], _, Upper, Upper, State, State).
adjust_constraints([Demand|Ds], UV, [P|Ps], Curr, Upper0, Upper, State0, State) :-
        Curr1 #= 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, ceiling(Demand/P)),
            adjust_constraints(Ds, UV, Ps, Curr1, Upper1, Upper, State1, State)
        ).



knapctrl(W, W*x(N0), N0, N) :- N #= N0 + 1.

knap_objective(Solved, _, Coeff*x(N0), N0, N) :-
        shadow_price(Solved, demand(N0), Coeff),
        N #= N0 + 1.

knap_integer(_, N0-S0, N-S) :-
        constraint(integral(x(N0)), S0, S),
        N #= N0 + 1.


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

knap_value(Solved, _, V, N0, N) :-
        variable_value(Solved, x(N0), V),
        N #= N0 + 1.