/***
aim: to allocate patients to linear accelerators for treatment.
patients may require linacs with certain capabilities.
patients may prefer to avoid certain times of the day.
each patient needs to be scheduled on exactly one machine (with sufficient capabilities)
a linac can only be used by up to one person at a time.
each linac only operates during a certain window during the day.
where not all patients' time preferences can be met, preference should be given to those with more flexibility
***/

set L; /* linacs */
set P; /* patients */
set C; /* capabilities */

param duration{p in P}, >0;  /* in minutes */
param firstTime{L}, >0;      /* in minutes after midnight */
param lastTime{L}, >0;       /* in minutes after midnight */

param nTimes, integer, > 0;

set TIMES, default 1..nTimes; /* set of preferred times */

param patient{i in TIMES};
param startWindow{i in TIMES};
param endWindow{i in TIMES};

param capabilities{L, C}, >= 0, binary;
param requirements{P, C}, >= 0, binary;

var lp{p in P, l in L}, >= 0, binary;    /* selected linac for this patient */

s.t. hasCapability{p in P, l in L, c in C}: 
    lp[p, l] * capabilities[l, c] >= lp[p, l] * requirements[p, c];

var x{p in P, l in L}, >= 0;   /* start time for this patient */

s.t. starty{p in P, l in L}:
    x[p, l] >= firstTime[l] * lp[p, l];

s.t. endy{p in P, l in L}:
    x[p, l] <= (lastTime[l] - duration[p]) * lp[p, l];

var Sequence{p in P, q in P, l in L : p < q}, binary;
/* 1 is p is before q */

param K := sum{l in L} lastTime[l];
/* big enough number */

s.t. foo{p in P, q in P, l in L: p < q}:
    x[p, l] >= (x[q, l] + duration[q]) - K * Sequence[p, q, l] - K * (1 - lp[p, l]) - K * (1 - lp[q, l]);

s.t. bar{p in P, q in P, l in L: p < q}:
    x[q, l] >= (x[p, l] + duration[p]) - K * (1 - Sequence[p, q, l]) - K * (1 - lp[p, l]) - K * (1 - lp[q, l]);

s.t. onemachine{p in P}:
    sum{l in L} lp[p, l] = 1;

param patient_weight{p in P}, >= 0 := 10000 / (1 + sum{i in TIMES} if p = patient[i] then (endWindow[i] - startWindow[i]) else 0);
/* the more times they want to avoid, the less we care about them */

var badness{p in P}, >= 0;

var afterMean{l in L, i in TIMES}, >= 0, binary;

/* true if x[p, l] > mean[i] */
s.t. xx1{l in L, i in TIMES}: x[patient[i], l] >= (startWindow[i] + endWindow[i]) / 2
                                                    - K * (1 - afterMean[l, i]);
s.t. xx2{l in L, i in TIMES}: x[patient[i], l] <= (startWindow[i] + endWindow[i]) / 2
                                                    + K * afterMean[l, i];

s.t. early{l in L, i in TIMES}: badness[patient[i]] >= (startWindow[i] + endWindow[i]) / 2
                                                       + (endWindow[i] - startWindow[i] ) / 2
                                                       - x[patient[i], l]
                                                       - K * (1 - afterMean[l, i])
                                                       - K * (1 - lp[patient[i], l]);
s.t. late{l in L, i in TIMES}: badness[patient[i]]  >= x[patient[i], l] + duration[patient[i]]
                                                       + (endWindow[i] - startWindow[i] ) / 2
                                                       - (startWindow[i] + endWindow[i]) / 2
                                                       - K * afterMean[l, i]
                                                       - K * (1 - lp[patient[i], l]);

var z;
s.t. weighted_badness{p in P}: z >= badness[p] * patient_weight[p];

minimize obj: z;

solve;

printf("\nResults\n");
for {l in L} {
    printf "%s (%dm utilisation):",
            l,
            sum{p in P} duration[p] * lp[p, l]
    ;
    for {c in C} {
        for {{0}: capabilities[l, c] == 1} {
            printf " %s", c;
        }
    }
    printf "\n";
    for {p in P} {
        for {{0}: lp[p, l] == 1} {
            printf "%s: %02d:%02d-%02d:%02d",
                    p,
                    round(x[p, l]) div 60,
                    round(x[p, l]) mod 60,
                    round(x[p, l] + duration[p]) div 60,
                    round(x[p, l] + duration[p]) mod 60
            ;
            printf ", %dm, weight %d, %f\n",
                    badness[p],
                    patient_weight[p],
                    weighted_badness[p]
                    ;
        }
    }
    printf "\n";
}
printf "z: %d\n", z;

data;

/* number of minutes required to administer treatment */
param : P : duration :=
    1 10
    2 15
    3 10
    4 15
    5 15
    6 10
    7 15
    8 10
    9 15
   10 15
   11 10
   12 15
   13 10
   14 15
   15 15;
   
/* linacs' available from now. 600=10am, 480=8am */
param : L : firstTime lastTime :=
    LINAC1 600  1200
    LINAC2 480  1200;

param : TIMES : patient startWindow endWindow :=
        1       1       800         900  /* ie: patient 1 should AVOID between midnight and 3pm */
        2       4       0           604
        3       4       610         1500
;

set C:= IMRT VMAT MRI;

param requirements : IMRT VMAT MRI :=
                 1   0    1    0
                 2   0    1    0
                 3   0    1    0
                 4   0    1    0
                 5   1    0    0
                 6   1    0    0
                 7   1    0    0
                 8   0    1    0
                 9   0    1    0
                 10  0    1    0
                 11  1    0    0
                 12  1    0    0
                 13  1    0    0
                 14  1    0    0
                 15  1    0    0
;

param capabilities : IMRT VMAT  MRI :=
                 LINAC1   1  0  0   /* linac 1 is IMRT only */
                 LINAC2   0  1  0   /* linac 2 is IMRT and VMAT capable */
;

end;