Three water jugs problem in MiniZinc

three_water_jugs.mzn

% Description of problem:
%
% https://mathworld.wolfram.com/ThreeJugProblem.html
%
% "Given three jugs with x pints in the first, y in the second, and z in the third,
%  obtain a desired amount in one of the vessels by completely filling up and/or
%  emptying vessels into others."
%
% Standard statement: 
%
% - Three jugs with capacities           8,5,3.
% - At the start they contain quantities 8,0,0.
% - The goal state is quantities         4,4,0.
% - You can only transfer water between the jugs, not refill or dump water on the floor.
%
% We find (using gecode solver, in ~1 sec):
%
% ---
% History step 1: 8 0 0 , transfer 5 gallon from jug 1 to jug 2
% History step 2: 3 5 0 , transfer 3 gallon from jug 2 to jug 3
% History step 3: 3 2 3 , transfer 3 gallon from jug 3 to jug 1
% History step 4: 6 2 0 , transfer 2 gallon from jug 2 to jug 3
% History step 5: 6 0 2 , transfer 5 gallon from jug 1 to jug 2
% History step 6: 1 5 2 , transfer 1 gallon from jug 2 to jug 3
% History step 7: 1 4 3 , transfer 3 gallon from jug 3 to jug 1
% History step 8: 4 4 0 , do nothing
% ---
%
% In this code, one can switch on "filling" and "dumping on floor"
% by setting the parameters can_fill, can_empty.
%
% Problems:
%
% - This code is probably not the best as I am not used to MiniZinc.
% - The ouput is generated with a constraint; this soaks up a lot of time.
%   Much better if one could do output in a loop, but that seems fastidious.
%
% Compare with code by Hakan Kjellerstrand at:
% http://hakank.org/minizinc/water_buckets1.mzn
%
% Author & License
%
% David Tonhofer ([email protected]) says: 
% This code is released under MIT license (https://opensource.org/licenses/MIT)

int  : max_step   = 10;    % try to do it in max steps, possibly less
bool : can_fill   = false;  % whether one can refill a jug fully from an infinite water source
bool : can_empty  = false;  % whether one can dump all the water water from a jug into an infinite water sink

% We have 3 jugs + 1 pseudo-jug: a very large water source / sink
% that exists so that we can have a constant sum over all the jug contents.

set of int : JUGS = {1,2,3,4}; 
set of int : JUGS_WITHOUT_PSEUDO = {1,2,3}; 
int : PSEUDO_JUG = 4;

array[JUGS] of int: start    = [8,0,0,1000];  % starting with these jug contents 
array[JUGS] of int: goal     = [4,4,0,1000];  % try to reach these jug contents 
array[JUGS] of int: capacity = [8,5,3,10000]; % with the jugs having these capacities (pseudo-jug is enormous)

int : total = sum(j in JUGS)(start[j]);

% Detect data entry problems

constraint assert(forall(j in JUGS) (start[j] <= capacity[j]), "some jug's start value is exceeding capacity");
constraint assert(forall(j in JUGS) (goal[j]  <= capacity[j]), "some jug's goal value is exceeding capacity");
constraint assert(sum(start) == sum(goal), "the sum of the start must be the sum of the goal");

% What we are building: a history of jug contents (a 2-D array)

array[1..max_step,JUGS] of var int: jug_history;

% --- Constraints ---

% We already know the start values.

constraint forall(j in JUGS) (jug_history[1,j] = start[j]);

% Jug contents are always within the respective jug's capacity

constraint forall(h in 1..max_step,j in JUGS)
   (0 <= jug_history[h,j] /\ jug_history[h,j] <= capacity[j]);

% Sum-of-jug contents is a constant (this helps in finding solutions)

constraint forall(h in 1..max_step) 
   (sum(j in JUGS)(jug_history[h,j]) == total);

% l is the number of steps until the goal state has been reached; we want to minimize it

var 1..max_step: l;

% Once the goal has been reached at step l, nothing changes anymore for the
% remainder of the history. 
% The constraint that the goal is actually not attained earlier than l is in
% the "all_different" constraint over the flattened[] array.

constraint forall(h in l..max_step, j in JUGS) 
   (jug_history[h,j] = goal[j]);

% We do not want to visit the same state twice, so build a unique int from
% each state (flatten each state into an int) then assert that all those
% integers are different. This should be more efficient than build a loop
% constraint that compares states pairwise.

array[1..max_step] of var int: flattened;

% This is hardcoded for 3 jugs and a max capacity of 10; can it be made flexible?

constraint forall(h in 1..max_step)
   (if (h<=l) then
      flattened[h] = (jug_history[h,1]*10 + jug_history[h,2])*10 + jug_history[h,3]
    else
      flattened[h] = -h % arbitrary but unique
    endif);

include "globals.mzn";
constraint all_different(flattened);

% Valid operations. These constraint us to valid histories only (you get the
% feel of quantum mechanics: disallow allow impossibilites). Note that there is no
% "do nothing" operation. All operations reduce to transfer between jug jx and jug jy,
% with the "emptying" operation a transfer into the pseudo-jug and the "filling"
% operation a transfer from the pseudo-jug.

constraint
   forall(h in 1..l-1)
    (
       exists(jx,jy in JUGS where jx != jy) (
          let { 
             var int: tt = min(capacity[jy] - jug_history[h,jy], jug_history[h,jx])
          } 
          in 
          (
             tt > 0
             /\
             (jy == PSEUDO_JUG -> can_empty)
             /\
             (jx == PSEUDO_JUG -> can_fill)
             /\
             jug_history[h+1,jx] == jug_history[h,jx] - tt
             /\
             jug_history[h+1,jy] == jug_history[h,jy] + tt
          )
          /\
          forall(j in JUGS where (j!=jx /\ j!=jy))
             (jug_history[h,j] == jug_history[h+1,j]) % not-involved jugs don't change
       ));

% For printout, also build a description of actions.
% This should not be done using a constraint but should use a "fixed"
% array marked output_only, how to do this properly???

array[1..max_step] of var ACTION : action;
array[1..max_step] of var int : amount;
array[1..max_step] of var int : from;
array[1..max_step] of var int : to;

enum ACTION = {DO_NOTHING,TRANSFER,FILL,EMPTY};

constraint forall(h in 1..max_step)
   ((h==max_step 
     \/ 
    (h<max_step -> (forall(j in JUGS) (jug_history[h,j] == jug_history[h+1,j]))))
    -> (action[h] = DO_NOTHING /\ amount[h] = 0 /\ from[h] = 0 /\ to[h] = 0));

constraint forall(h in 1..max_step-1)
   (forall(jx in JUGS, jy in JUGS) (
      let { var int : d = jug_history[h,jx] - jug_history[h+1,jx] } 
      in (
         (d > 0 /\ jug_history[h,jy] - jug_history[h+1,jy] = -d)
         ->
         (
            (amount[h] = d /\ from[h] = jx /\ to[h] = jy)
            /\
            ((jx != PSEUDO_JUG /\ jy != PSEUDO_JUG) -> (action[h] = TRANSFER))
            /\
            ((jx == PSEUDO_JUG /\ jy != PSEUDO_JUG) -> (action[h] = FILL)) 
            /\
            ((jx != PSEUDO_JUG /\ jy == PSEUDO_JUG) -> (action[h] = EMPTY)) 
         ))));

% That's it.

solve minimize l;

output [
   "History step \(h): "
   ++
   concat([ "\(jug_history[h,j])" ++ if (j<max(JUGS)) then " " else "" endif | j in JUGS_WITHOUT_PSEUDO ])
   ++ 
   ", "
   ++
   (
   if     (fix(action[h]) == DO_NOTHING) 
   then   "do nothing"
   elseif (fix(action[h]) == TRANSFER)    
   then   "transfer \(amount[h]) gallon from jug \(from[h]) to jug \(to[h])"
   elseif (fix(action[h]) == FILL)    
   then   "fill jug \(to[h]) with \(amount[h]) gallon"
   elseif (fix(action[h]) == EMPTY)    
   then   "empty jug \(from[h]) of \(amount[h]) gallon"   
   else   "???"
   endif
   )
   ++
   "\n"
   | h in 1..fix(l) ];