/***
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 */

param nPatients, >0, integer;
set P := 1..nPatients ; /* patients */

set C; /* capabilities */

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 duration{p in P}, >0;
/* in minutes */

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}, 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;

param nPatients := 15;
set L := LINAC1 LINAC2;

/* linacs' available from now. 600=10am, 480=8am */
param firstTime :=
    [LINAC1] 600
    [LINAC2] 480;

param lastTime :=
    [LINAC1] 1200
    [LINAC2] 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
;

/* number of minutes required to administer treatment */
param 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;


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;