Saturday, September 8, 2018

Minizinc CP/SAT vs MIP

In a previous post I discussed a simple scheduling model [1]: assign tasks to time slots subject to capacity constraints, such that we minimize the number of time slots with activity. A simple MIP model was developed for this along the lines of: \[\bbox[lightcyan,10px,border:3px solid darkblue]{\begin{align}\min& \sum_t y_t\\ & \sum_{t| \mathit{ok}_{i,t}}  x_{i,t} = 1 && \forall i \\ & \sum_{i| \mathit{ok}_{i,t}} x_{i,t} \le N && \forall t \\ & y_t \ge x_{i,t} &&\forall i,t|\mathit{ok}_{i,t}\\ & x_{i,t}, y_t \in \{0,1\}\end{align}}\] As we can see, not all task-slot combinations \((i,t)\) are allowed. The sample data set to allowed assignments is:


task_availability_map = {
    "T1" : [0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    "T2" : [0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    "T3" : [0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    "T4" : [0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    "T5" : [0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    "T6" : [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
    "T7" : [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0],
    "T8" : [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0],
    "T9" : [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
    "T10": [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0]
}

In this case we see the allowed assignments could be encoded as min, max tuples, as we have contiguous ranges. I will assume the more general case, where we only have unstructured, sparse data \(ok_{i,t}\) indicating if an assignment is allowed.


The question came up if it would be easier to formulate this as a Constraint Programming problem. There are two basic approaches we can use:


  1. Use binary variables \(x_{i,t}\). When we use a max function, we can write:\[\begin{align} \min & \sum_t \max_i \left\{ a_{i,t} x_{i,t}\right\}  \\ & \sum_t a_{i,t}  x_{i,t} = 1 && \forall i \\ & \sum_i a_{i,t} x_{i,t} \le N && \forall t \\ & a_{i,t}=0 \Rightarrow x_{i,t}=0 \\ & x_{i,t} \in \{0,1\}\end{align}\] Here \(a_{i,t}\) is our availability parameter. This generates more variables than needed. 
  2. Use binary variables but squeeze out the decision variables that we know ahead of time to be zero. This leads to \[\begin{align} \min & \sum_t \max_{i| \mathit{ok}_{i,t}} \left\{ x_{i,t} \right\} \\ & \sum_{t| \mathit{ok}_{i,t}}  x_{i,t} = 1 && \forall i \\ & \sum_{i| \mathit{ok}_{i,t}} x_{i,t} \le N && \forall t \\ & x_{i,t} \in \{0,1\}\end{align}\] This is close to the MIP model, but we do without the \(y\) variables.
  3. Use integer variables \(x_i\) to denote when task \(i\) is scheduled. \[\begin{align} \min & \sum_t \max_i \left\{ I(x_i=t) \right\} \\ & \sum_i I(x_i=t) \le N && \forall t \\&x_i \in \{t | ok_{i,t}\} \end{align}\] where \(I(.)\) is the indicator function. I.e. \[I(\delta) = \begin{cases} 1 & \text{if $\delta$ is true} \\ 0 & \text{otherwise}\end{cases}\]

First let's review the reference implementation in GAMS:

sets
  i
'task' /task1*task10/
  t
'time slot' /t1*t17/
;


*
* data
*

table allowed(i,t)
           
t1  t2  t3  t4  t5  t6  t7  t8  t9 t10 t11 t12 t13 t14 t15 t16 t17
   
task1    0   0   0   1   1   1   1   1   0   0   0   0   0   0   0   0   0
   
task2    0   0   0   0   0   1   1   0   0   0   0   0   0   0   0   0   0
   
task3    0   0   0   1   1   1   1   1   0   0   0   0   0   0   0   0   0
   
task4    0   0   0   0   0   1   1   1   0   0   0   0   0   0   0   0   0
   
task5    0   0   0   0   1   1   1   1   0   0   0   0   0   0   0   0   0
   
task6    1   1   1   1   1   1   1   1   1   1   1   1   1   1   0   0   0
   
task7    0   0   0   0   0   0   0   0   0   0   1   1   1   1   1   0   0
   
task8    0   0   0   0   0   0   0   0   0   0   1   1   1   1   0   0   0
   
task9    0   0   0   0   0   0   1   1   1   1   0   0   0   0   0   0   0
   
task10   0   0   0   0   0   0   1   1   1   1   0   0   0   0   0   0   0
;

scalar N 'capacity of a slot' /3/;

*
* derived data
*
set ok(i,t) 'set version of allowed parameter';
ok(i,t) = allowed(i,t);

scalar nok 'number of allowed assignments';
nok =
card(ok);
display nok;


*
* model
*
binary variables
  x(i,t)
'assignment'
  y(t)
'slot is used'
;
variable z 'objective variable';

equations
   assign1(i) 
'assignment constraint'
   assign2(t) 
'assignment constraint'
   ybound(i,t)
'x(i,t)=1 => y(t)=1'
   obj        
'objective'
;

assign1(i)..
sum(ok(i,t), x(i,t)) =e= 1;
assign2(t)$
sum(ok(i,t),1).. sum(ok(i,t), x(i,t)) =l= N;
ybound(ok(i,t))..  y(t) =g= x(i,t);
obj.. z =e=
sum(t,y(t));

model m /all/;

*
* solve model to optimality
*
option optcr=0;
solve m minimizing z using mip;

*
* display solution
*
option x:0,y:0,z:0;
display x.l, y.l, z.l;


Notes:

  • The parameter allowed and the set ok are stored sparsely: only the nonzero elements are stored. Each has 50 nonzero elements. For a small data set like this, this is not an issue. For larger data sets this can become very important. Note that the parameter allowed or set ok will typically come from a database or other data source. A database table will often be organized in sparse "long" format: e.g.for this case a table with columns task and slot with just 50 records. 
  • In GAMS the model size is determined by the variables that occur in the model, so the number of variables is \(n=68\) (\(x\): 50, \(y\): 17, and \(z\):1). The declaration x(i,t) does not allocate any variables.
  • The number of equations is \(m=76\)  (assign1: 10, assign2: 15, ybound: 50, obj: 1). Note that assign2 automatically skips the equations related to the "empty" columns 16 and 17 in the data table. Those columns lead to constraints of the form \(0 \le 3\). 

We have a number of issues in the model to take care of:

  • sparse data 
  • sparse \(x_{i,t}\) variables
  • empty constraints due to empty columns in the data

The GAMS model takes care of these issues automatically.

Sparsity in Minizinc


It is interesting to see what Minizinc does with a sparse model. Here is a small example:


int: N = 100;
set of int: I = 1..N;
set of int: J = 1..N;
array[I,J] of var int: x;

constraint x[1,1] = 1;
solve satisfy;

*---------------------------

Compiling sparse.mzn
Generated FlatZinc statistics:
Variables: 9999 int
Constraints: none

I somewhat expected to see a model with 1 variable and 1 constraint. We get 9999 variables: the one variable that is fixed is removed.

Other modeling tools, especially based on a programming language like Python, can create variables using loops (or list comprehensions). This allows us to skip certain variables.

Binary variable model


The easiest is just to forget about sparsity and just generate and solve fully allocated \(x_{i,t}\) model in Minizinc. This is model 1 mentioned above.


include "globals.mzn";

% number of tasks, time slots
int: ntasks = 10;
int: nslots = 17;

% sets
set of int: I = 1..ntasks;
set of int: T = 1..nslots;

array[I,T] of int: available =
   [| 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    | 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    | 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    | 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    | 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    | 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0,
    | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0,
    | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0,
    | 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
    | 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0 |];
    
% max tasks per time slot
int: N = 3;    

array[I,T] of var 0..1: x;
var int: count;

constraint forall (i in I) ( sum([ available[i,t]*x[i,t] | t in T ] ) = 1);

constraint forall (t in T) ( sum([ available[i,t]*x[i,t] | i in I ] ) <= N );

constraint count = sum([ max([ available[i,t]*x[i,t] | i in I ]) | t in T ]);

% in this dense model we need to add:
%    available[i,t]=0 => x[i,t] = 0
constraint forall (i in I, t in T where available[i,t]=0) (x[i,t] = 0);  

solve minimize count;

Of course this is a small data set, and the standard Minizinc solvers have no problem in solving this.

It becomes a little bit more interesting with somewhat larger data set. Here is a random problem with 50 time slots and 50 tasks. The results look like:

Problem with 50 tasks and 50 time slots

Minizinc's CP solvers seem to have a big problem with this. They can not confirm the optimal solution of 17 time slots in a reasonable amount of time.

When we look at the solution in the above picture we see that the used capacity in each column is:


----     70 PARAMETER use  used capacity

t6  3,    t7  3,    t11 3,    t13 3,    t14 3,    t17 3,    t18 3,    t21 3,    t25 3,    t26 3,    t32 3,    t33 3
t34 3,    t39 3,    t44 2,    t49 3,    t50 3

All but one selected column have a used capacity of 3 (slot t44 has a count of 2). This means we will never be able to get to 16 time slots. In other words the minimum number of time slots is\[\left\lceil \frac{\mathit{ntasks}}{N}\right\rceil\] (the brackets indicate the ceiling function, i.e. rounding upwards). We can add this lower bound as an explicit constraint to the problem to help the solver. [2]

Minizinc also supports MIP solvers. The CBC MIP solver takes about 5 seconds to solve this model to optimality. Minizinc will automatically translate the max function into something the MIP solver can handle.

I have tried to implement the original squeezed binary variable model (i.e. model 2) in Minizinc. Here we try to skip all variables \(x_{i,t}\) where we know it is zero.  My version (with the small data set) looks like:


include "globals.mzn";

% number of tasks, time slots
int: ntasks = 10;
int: nslots = 17;

% sets
set of int: I = 1..ntasks;
set of int: T = 1..nslots;

% this is stored as a dense matrix
array[I,T] of int: available =
   [| 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    | 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    | 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    | 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    | 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    | 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0,
    | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0,
    | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0,
    | 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
    | 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0 |];
    
% max tasks per time slot
int: N = 3;    

% derived data: column sums of available[i,t].
array[T] of int: sumavail = [sum([available[i,t] | i in I]) | t in T]; 
    
            
% I really want variables x[i,t] only when available[i,t]=1. 
% I don't think we can have sparse arrays. So here we
% just declare the whole thing.
array[I,T] of var 0..1: x;
var int: count;

constraint forall (i in I) ( sum([ x[i,t] | t in T where available[i,t]=1 ] ) = 1);

constraint forall (t in T) ( sum([ x[i,t] | i in I where available[i,t]=1 ] ) <= N );

constraint count = sum([ max([ x[i,t] | i in I where available[i,t]=1 ]) | t in T where sumavail[t]>=1]);

solve minimize count;

Note that the count constraints suddenly needs extra attention. The reason is that we cannot take the max of an empty list. We protected against that. A different way to do this is to add a zero to the list:

constraint count = sum([ max([0] ++ [ x[i,t] | i in I where available[i,t]=1 ]) | t in T ]);

The results for the 50 task, 50 time slot data set are:

ToolModelSizeTimeObjNodes
Minizinc/GecodeModel 1Variables: 330 int
Constraints: 144 int
> 3 hours
Minizinc/CBCModel 1Variables: 886 int, 607 float
Constraints: 420 int, 1165 float
4.5 s1710
Minizinc/CBCModel 2Variables: 885 int, 606 float
Constraints: 420 int, 1163 float
3.2 s170
GAMS/CBCMIP Modelsingle equations: 382
single variables: 334
1.2 s170

I am not sure about the Minizinc size statistics. They depend on the solver type being used: the MIP solvers need more variables. We could look at the generated "flatzinc" file which is passed on to the solver to study what is going on. (I looked at it: well, there is much going on, more than meets the eye).  I expected the size differences between model 1 and 2 using CBC to be very different, but they are very close. The solution path is also somewhat different. Not sure if I understand this.

It is interesting to see the differences in approaches between GAMS and Minizinc.

  • GAMS handles large, sparse models easily,
  • Minizinc allows some higher level modeling abstractions, such as max() and count()

Sparse data and sparse variables and constraints arise often in practical models. The fact that Minizinc does not handle sparse data and variables directly, does not mean we have to give up. E.g. the data available can  be organized in a different way:


include "globals.mzn";

% number of tasks, time slots
int: ntasks = 10;
int: nslots = 17;

% sets
set of int: I = 1..ntasks;
set of int: T = 1..nslots;
   
int: nitems = 50; 
set of int: ITEMS = 1..nitems;
array[ITEMS] of int: task = 
   [1,1,1,1,1,2,2,3,3,3,3,3,4,4,4,5,5,5,5,6,6,6,6,6,
    6,6,6,6,6,6,6,6,6,7,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,10];
array[ITEMS] of int: slot = 
  [4,5,6,7,8,6,7,4,5,6,7,8,6,7,8,5,6,7,8,1,2,3,4,5,6,7,
   8,9,10,11,12,13,14,11,12,13,14,15,11,12,13,14,7,8,9,10,7,8,9,10];
   
            
% max tasks per time slot
int: N = 3;      
            
array[ITEMS] of var 0..1: x;
var int: cnt;

constraint forall (i in I) ( sum([ x[k] | k in ITEMS where task[k]=i ] ) = 1);

constraint forall (t in T) ( sum([ x[k] | k in ITEMS where slot[k]=t ] ) <= N );

% bummer: this does not work
% constraint cnt = sum([ max([ x[k] | k in ITEMS where slot[k]=t ]) | t in T where count(slot,t)>=1]);
% work around:
array[T] of int : sumavail = [ sum([1 | k in ITEMS where slot[k]=t ]) | t in T ];
constraint cnt = sum([ max([ x[k] | k in ITEMS where slot[k]=t ]) | t in T where sumavail[t]>=1]);


solve minimize cnt;


We use here two parallel arrays task and slot to encode the sparse matrix available. Of course, we make multiple passes over the data here with these nested list comprehensions. So this becomes inefficient for very large data sets (but solvers may get into troubles then also).

Integer variable formulation


A Minizinc formulation with integer variables (model 3) can look like:

array[I] of var T: x;
var int: cnt;

constraint forall (i in I) ( x[i] in {t | t in T where available[i,t]=1 } );
constraint forall (t in T) ( count(x,t) <= N);

constraint cnt = sum( [count(x, t) >= 1 | t in T] );

solve minimize cnt;

Let's try the 50 task, 50 time slot data set. Some of the Minizinc solvers find the solution with the optimal objective cnt=17 quite quickly. But proving this is indeed the optimal solution requires to prove that a solution with cnt=16 does not exist. The CP and SAT solvers have a difficult time with this. However the built-in MIP solvers actually do a good job on this. To be able to solve this with a MIP solver Minizinc has to do some non-trivial transformations.

When we add the lower bound \[cnt \ge \left\lceil \frac{\mathit{ntasks}}{N}\right\rceil\] most solvers can prove the optimal solution in a reasonable time. The MIP solvers don't seem to need this bound: they do a good job without it.

The model above invokes the count function twice for each \(t\). We could write:


array[I] of var T: x;
array[T] of var 0..N: countx;
var int: cnt;

constraint forall (t in T) ( countx[t] = count(x,t) );

constraint forall (i in I) ( x[i] in {t | t in T where available[i,t]=1 } );

constraint cnt = sum( [ 1 | t in T where countx[t] >= 1] );

solve minimize cnt;

This does not seem to make too much different. The built-in CBC MIP solver can handle these models with ease. In [1] we see that MIP solvers can handle even sizes like 200 by 200 quite easily (that was the size suggested by the original poster).

Update


Not all solvers are created equal. Laurent Perron shows in [2] how to solve this problem very quickly using or-tools in Python (0.9 seconds). Also note the idea of adding a lower bound on the objective.

References



Appendix: 50 tasks, 50 time slots data set



% number of tasks, time slots
int: ntasks = 50;
int: nslots = 50;

% sets
set of int: I = 1..ntasks;
set of int: T = 1..nslots;

array[I,T] of int: available =
  [|0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,0,
   |0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,
   |0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,
   |0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,
   |0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 |];

% max tasks per time slot
int: N = 3;