% --- Data quantities ---
enum Patient; % Patients for which samples have been taken
set of int: Horizon = -5..2; % (Relative) sampling times

% --- Data input ---
float: DLMO_threshold = 13.043;
array[Patient, Horizon] of float: sample;

% --- Derived data ---

% For each patient: after which sample could we stop and have found dlmo? (round must then include `critical_test[PT]` and its successor to catch it).
array[Patient] of var Horizon: critical_sample = [ min([i | i in -5..1 where sample[pt, i] <= DLMO_threshold /\ sample[pt, i+1] >= DLMO_threshold]) | pt in Patient];

% --- Decisions ---
% First round:
var Horizon: first_sample; % When to start sampling on the first round 
var 2..length(Horizon): first_num_samples; % How many samples to take on the first round 

% Second round:
var 2..length(Horizon): second_num_samples; % How many samples to take on the second round

% (Assumption: If the levels were to high in all samples during the first round, then we sample `second_num_samples` ending with the `first_sample`.
%              If the levels were to low in all samples during the first round, then we sample `second_num_samples` from `first_sample + first_num_samples`.)

% --- Decision Consequences ---

% For each patient: will their DLMO be detected on the first round?
array[Patient] of var bool: first_detected = [
	critical_sample[pt] >= first_sample /\ critical_sample[pt] + 1 <= first_sample + (first_num_samples - 1)
	| pt in Patient
];

% --- Constraint ---

% The DLMO should be detected on the second round (if not detected on the first round)
constraint forall (pt in Patient) (
	first_detected[pt] 
	% Before first round 
	\/ (critical_sample[pt] >= first_sample - (second_num_samples - 1) /\ critical_sample[pt] + 1 <= first_sample)
	% After first round
	\/ (critical_sample[pt] >= first_sample + (first_num_samples - 1) /\ critical_sample[pt] + 1 <= first_sample + (first_num_samples - 1) + (second_num_samples - 1))
);

% --- Objective ---
solve minimize 
	% Number of tests on first round
	(length(Patient) * first_num_samples)
	% Number of test on second round
	+ (count(pt in Patient) (not first_detected[pt]) * second_num_samples);

% --- Output ---

output ["First Round:\n"]
++ [format(2, i) ++ "   " | i in Horizon] ++ ["\n |"]
++ [if i in fix(first_sample)..fix(first_sample + (first_num_samples-1) - 1) then "====|" else "----|" endif | i in -5..1]
++ ["\n\nFirst try (%): \(count(first_detected) / length(Patient) * 100)\n"]
++ ["\n\nSecond Round:\n"]
++ [format(2, i) ++ "   " | i in Horizon] ++ ["\n |"]
++ [if i in fix(first_sample-(second_num_samples-1))..fix(first_sample-1) \/ i in fix(first_sample + (first_num_samples - 1))..fix(first_sample + (first_num_samples - 1) + (second_num_samples - 1) - 1) then "====|" else "----|" endif | i in -5..1]
++ ["\n\nTotal number of tests = \(_objective);"]


++ ["\n\n\nfirst_sample = \(first_sample);\n"]
++ ["first_num_samples = \(first_num_samples);\n"]
++ ["second_num_samples = \(second_num_samples);"];