plate-optimizer.mzn

% Copyright 2026 COMPD Authors.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     http://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.
%
%
% Description:  A microplate layout designer represented as an optimization problem
%
% Authors: Ramiz GINDULLIN (ramiz.gindullin@it.uu.se)
% Version: 1.3.4
% Last Revision: February 2026
%
% =======================================================================================
% COORDINATE SYSTEM NOTATION - READ THIS FIRST!
% =======================================================================================
% This model uses TWO coordinate systems:
%
% 1. INNER COORDINATES (used in all constraints):
%    - Ranges: 1..num_rows_line × 1..num_cols_line
%    - Excludes outer edge wells (size_empty_edge) and corner empty wells
%    - ALL constraint logic operates in this coordinate space
%
% 2. FULL PLATE COORDINATES (used only in output):
%    - Ranges: 1..num_rows × 1..num_cols  
%    - Includes edge and corner wells
%    - Conversion: calculate_line_index() maps inner to full coordinates
%
% CRITICAL: When modifying constraints, always work in INNER coordinates!
%           Only the output section uses full plate coordinates.
% =======================================================================================

include "globals.mzn";

%-------------------COMPD parameters--------------------------
%%% Tune to change performance / quality of a solution %%%
bool:  flag_consolidate_controls = true;
       % if true  - all controls are placed in a single criteria set (unless they cover more than half of the plate line)
       % if false - each control is placed in an individual criteria set 
       % RATIONALE - when set to false, controls of various types can be grouped together,  e.g., positive and negative controls 
       % can be neighbours, capturing the same plate effect which can be advantegous in some situations
bool:  randomize_order_of_material_concentrations = true;
       % if true  - the order of concentrations within compounds is randomized. Helps with an edge case when the
       %            replicate numbers for compounds are divisible by 4
       %            (see emptywells_controls_compounds_names_concentrations declaration)
       % if false - default order.
int:   num_threshold_use_fake_edge_wells = 0;
       % if the criteria set for a compound contains from 2 to num_threshold_use_fake_edge_wells replicates, 
       % apply fake edge wells when calculating the distances for the criteria
       % Set to 4 or 6 when used with sparse layouts (i.e. layouts with a relatively small number of compounds and a lot of empty space)
       % Set to 0 otherwise
int:   num_wells_set_criteria            = if num_total_wells_line < 300 then 90   else 120  endif;
float: num_threshold_min_dist_criteria   = if num_total_wells_line < 300 then 0.35 else 0.25 endif;
       % RATIONALE for num_wells_set_criteria and num_threshold_min_dist_criteria: 
       %            With many replicates (e.g., >25% of plate line), distance optimization has
       %            diminishing returns because:
       %            1. Quadrant allocation already forces spatial spread (not clustered)
       %            2. Computational cost (O(n²) distances) grows rapidly
       %            3. Downstream statistical correction (LOESS, GAM) can handle residual
       %               plate effects more effectively than micro-optimizing well placement
       %            NOTE: Users can increase threshold via num_wells_set_criteria if needed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%-------------------FUNCTIONS--------------------------
function int: compound_index(int: cmpds, string: cmpd, array[int] of string: cmpds_names) =
                            min([i | i in 1..cmpds where cmpd = cmpds_names[i]]);

function array[int] of int: random_permutation(int: n) = arg_sort([cauchy(0,2) | i in 1..n]);

function var int: distance_manhattan(array[1..2] of var int: P1, array[1..2] of var int: P2) =
                            abs(P1[1] - P2[1]) + abs(P1[2] - P2[2]);

function var int: calculate_line_index(bool: edge_type, int: edge_size, int: line_size, var int: line, var int: coordinate) =
         if edge_type then (line - 1) * (2 * edge_size + line_size) + edge_size + coordinate
                      else (line - 1) *                  line_size  + edge_size + coordinate endif;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%-------------------INPUT DATA--------------------------
%% Information about constraints %%
bool: allow_empty_wells;
bool: replicates_on_different_plates;
bool: replicates_on_same_plate;
bool: concentrations_on_different_rows;
bool: concentrations_on_different_columns;

%% Information about the layout %%
int: horizontal_cell_lines;
int: vertical_cell_lines;
int: size_empty_edge;
opt bool: inner_empty_edge_input;
bool: inner_empty_edge = inner_empty_edge_input default true;

%% Compounds %%
int: compounds; %% number of drugs/compounds
array [1..compounds] of int: compound_replicates;
array [1..compounds] of int: compound_concentrations;
int: max_compound_concentrations = max(compound_concentrations++[0]);

array[1..compounds] of string: compound_names;
array[1..compounds,1..max_compound_concentrations] of string: compound_concentration_names;

%% Combinations (DEPRECATED) %%
% These variables are maintained solely for compatibility with PLAID data file format
% They are NOT used anywhere in the current model logic
% Removing them would break existing data files generated by PLAID
int: combinations;
int: combination_concentrations;
array[1..combinations] of string: combination_names;
array[1..combination_concentrations] of string: combination_concentration_names;

%% Information about controls %%
int: num_controls;
array [1..num_controls] of int: control_replicates; 
array [1..num_controls] of int: control_concentrations;
int: max_control_concentrations = max(control_concentrations++[0]);
array[1..num_controls,1..max_control_concentrations] of string: control_concentration_names;
array[1..num_controls] of string: control_names;

int: total_controls = sum([control_concentrations[i]*control_replicates[i] | i in 1..num_controls]);

%% Plate size / number of wells
int: num_rows;
int: num_cols;
opt int: size_corner_empty_wells;
int: num_corner_empty_wells = size_corner_empty_wells default 0;

%% Deprecated
array[int] of string: compound_concentration_indicators;

%%% Testing %%%
opt bool: testing = false;
opt bool: print_all = false;
bool: debugging = testing /\ print_all default false;
opt int: swap_search;
opt bool: sorted_compounds;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Datafile validation %%%
constraint assert(compounds >= 0,"Invalid datafile: Number of compounds cannot be less than zero.");
constraint assert(combinations >= 0,"Invalid datafile: Number of combinations cannot be less than zero.");
constraint assert(num_controls >= 0,"Invalid datafile: Number of controls should not be less than zero.");
constraint assert(vertical_cell_lines > 0,"Invalid datafile: Number of cell lines should be larger than zero.");
constraint assert(horizontal_cell_lines > 0,"Invalid datafile: Number of cell lines should be larger than zero.");
constraint assert(num_rows_line > 0,"Invalid datafile: Number of rows should be larger than zero.");
constraint assert(num_cols_line > 0,"Invalid datafile: Number of columns should be larger than zero.");
constraint assert((replicates_on_different_plates /\ replicates_on_same_plate) == false,"Invalid datafile: replicates cannot be both on the same plate and on different plates");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Calculating dimensions of a line %%%
% These define the INNER COORDINATE SYSTEM used by all constraints
% num_rows_line, num_cols_line = usable wells per plate (excluding edges/corners)
% This is distinct from num_rows, num_cols (full plate dimensions)
int: num_rows_line = if inner_empty_edge
                   then floor(num_rows / horizontal_cell_lines) - 2 * size_empty_edge
                   else floor((num_rows - 2 * size_empty_edge) / horizontal_cell_lines) endif;
int: num_cols_line = if inner_empty_edge
                   then floor(num_cols / vertical_cell_lines) - 2 * size_empty_edge
                   else floor((num_cols - 2 * size_empty_edge) / vertical_cell_lines) endif;

constraint assert((2*num_corner_empty_wells <= max([num_rows_line, num_cols_line])),
                   "Size of corners exceeds dimension/-s if the plate:\nnumber of rows per line = \(num_rows_line)\nnumber of columns per line = \(num_cols_line)\nthe size of a corner = \(num_corner_empty_wells)");
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%-------------------PREPROCESSING - plate and quadrant distribution--------------------------

% calculate the number of forced empty wells both per quadrant and per line
int: num_forced_empty_wells_quadrant = num_corner_empty_wells * num_corner_empty_wells;
int: num_forced_empty_wells_total = 4 * num_forced_empty_wells_quadrant;

int: num_total_wells_line = num_rows_line * num_cols_line - num_forced_empty_wells_total;
constraint assert(num_total_wells_line > 1, "Model ERROR! Number of wells within a plate line is \(num_total_wells_line) < 1. This should never happen!");

%% Number of wells needed. Note that plates might not be full
int: num_total_wells = sum([compound_concentrations[i]*compound_replicates[i] | i in 1..compounds]) + total_controls; 

%% Number of plates needed
%% max is used to avoid division-by-zero errors
int: num_plates_lines = max(ceil(num_total_wells / num_total_wells_line), 1);

int: num_empty_wells = num_total_wells_line * num_plates_lines - num_total_wells;

%%%%% Data validation %%%%%
constraint assert(num_empty_wells >= 0, "Model ERROR! Inner empty wells is negative. This should never happen!");
constraint assert(((not allow_empty_wells) -> (num_empty_wells = 0)), "The number of inner empty wells is \(num_empty_wells) while none are allowed.\n");
constraint assert(num_total_wells > 0, "Invalid data: the plates cannot be completely empty.");
constraint assert(num_total_wells_line > 0, "Invalid data: There are no wells on the plate.");
constraint assert(min(compound_concentrations++[0]) <= num_total_wells_line, "Invalid data: Number of concentrations does not fit in one plate. If you think this is a mistake, please contact the development team.");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

int: len_emptywells_controls_compounds = 1 + sum(control_concentrations) + sum(compound_concentrations);


% Pre-compute random permutations for each control's concentrations
% This ensures names and quantities use the SAME permutation.
% See the declation of emptywells_controls_compounds_categories for explanation
array[1..num_controls] of array[int] of int: control_permutations_shuffled = [random_permutation(control_concentrations[i]) | i in 1..num_controls];
array[1..num_controls] of array[int] of int: control_permutations_ordered =  [[j | j in 1..control_concentrations[i]] | i in 1..num_controls];
array[1..num_controls, 1..max(control_concentrations)] of int: control_concentration_permutations =
  array2d(1..num_controls, 1..max(control_concentrations),
          [if j <= control_concentrations[i] 
           then (if randomize_order_of_material_concentrations 
                 then control_permutations_shuffled[i][j]
                 else control_permutations_ordered[i][j]
                 endif)
           else 0
           endif
           | i in 1..num_controls, j in 1..max(control_concentrations)]);

% Pre-compute random permutations for each compound's concentrations
array[1..compounds] of array[int] of int: compound_concentration_permutations_shuffled = [random_permutation(compound_concentrations[i]) | i in 1..compounds];
array[1..compounds] of array[int] of int: compound_concentration_permutations_ordered =  [[j | j in 1..compound_concentrations[i]] | i in 1..compounds];
array[1..compounds, 1..max(compound_concentrations)] of int: compound_concentration_permutations =
  array2d(1..compounds, 1..max(compound_concentrations),
          [if j <= compound_concentrations[i] 
           then (if randomize_order_of_material_concentrations 
                 then compound_concentration_permutations_shuffled[i][j]
                 else compound_concentration_permutations_ordered[i][j]
                 endif)
           else 0
           endif
           | i in 1..compounds, j in 1..max(compound_concentrations)]);


% IMPORTANT: We assume that controls and compounds are placed sequentially. i.e. we can't have ["ctrs", "comp", "ctrs"] or ["comp 1", "comp 2", "comp 1"] etc
%   A reasonable assumption, as these arrays are generated authomatically
%   BUT! If this assumption is broken then some of the later logic will be broken. Namely WRT array "emptywells_controls_compounds_criteria_ordered"
array[1..len_emptywells_controls_compounds] of string: emptywells_controls_compounds_categories =
         ["ew"] ++ ["ctrs" | i in 1..sum(control_concentrations)] ++ ["comp" | i in 1..sum(compound_concentrations)];
array[1..2,1..len_emptywells_controls_compounds] of string: emptywells_controls_compounds_names_concentrations =
         array2d(1..2,1..len_emptywells_controls_compounds,
                 ["ew"] ++
                 [control_names[i]  | i in 1..num_controls, j in 1..control_concentrations[i]] ++
                 [compound_names[i] | i in 1..compounds,    j in 1..compound_concentrations[i]] ++
                 [""] ++
                 [control_concentration_names[i,  control_concentration_permutations[i,j]]  | i in 1..num_controls, j in 1..control_concentrations[i]] ++
                 [compound_concentration_names[i, compound_concentration_permutations[i,j]] | i in 1..compounds,    j in 1..compound_concentrations[i]]
                 % random_permutation(*_concentrations[i]) is because we shuffle the list of concentrationsto randomize their assignments within quadrants.
                 % It is useful in the situations when the number of replicates for each material-concentration is divisible by 4:
                 % in that situation, all concentrations of the same dose would be in the same quadrant, which is not desirable.
                 % It is a very specific edge case, but it needs to be accounted for.
                 % The resulting list would like like this: [dr1-ct2, dr1-ct1, dr1-ct3, dr2-ct3, dr2-ct2, dr2-ct1, dr3-ct1, dr3-ct3, ...]
                 % This order must be matched in the declaration of array emptywells_controls_compounds_qtys_total (if a change is introduced,
                 % where the number of replicates between concentrations differs.
                 );

% collect all quantitites before assigning them between plates / lines
array[1..len_emptywells_controls_compounds] of int: emptywells_controls_compounds_qtys_total =
                 [num_empty_wells] ++
                 
                 % If the number of replicates depends on the concentration (unlike the current data files),
                 % then swap "control_replicates[i]" to "control_replicates[i, control_concentration_permutations[i,j]]" in the line below
                 [control_replicates[i]  | i in 1..num_controls, j in 1..control_concentrations[i]] ++
                 
                 % If the number of replicates depends on the concentration (unlike the current data files),
                 % then swap "compound_replicates[i]" to "compound_replicates[i, compound_concentration_permutations[i,j]]" in the line below
                 [compound_replicates[i] | i in 1..compounds, j in 1..compound_concentrations[i]];

% these parameters are only used for replicates_on_same_plate strategy when num_plates_lines > 1
int: compounds_per_plate_min = floor(compounds div num_plates_lines);
int: compounds_per_plate_mod = compounds mod num_plates_lines;
int: emptywells_controls_num = 1 + sum(control_concentrations);
array[1..compounds] of 1..num_plates_lines: compounds_on_plates_lines =
                       [p | p in 1..num_plates_lines, i in 1..(compounds_per_plate_min + (compounds_per_plate_mod >= p))];
array[1..num_plates_lines] of int: compounds_per_plates_lines_total =
                       [sum([compound_replicates[i] * compound_concentrations[i]
                             | i in 1..compounds where compounds_on_plates_lines[i] = p])
                        | p in 1..num_plates_lines];
constraint assert((replicates_on_same_plate /\ (num_plates_lines > 1)) ->
                  (forall(i in 1..num_plates_lines)(compounds_per_plates_lines_total[i] <= num_total_wells_line)),
                  "Model ERROR! Default behavior for compounds_on_plates_lines strategy assigns too many replicates on one plate.\n\(show2d(emptywells_controls_compounds_names_concentrations))\n\(show(emptywells_controls_compounds_qtys_total))\n");
array[1..num_plates_lines] of int: available_wells_per_plates_lines_total =
                       [num_total_wells_line - compounds_per_plates_lines_total[p] | p in 1..num_plates_lines];
int: available_wells_per_plates_lines_total_sum = sum(available_wells_per_plates_lines_total);
int: emptywells_controls_compounds_first_non_zero =
                 min([i | i in 1..len_emptywells_controls_compounds
                          where emptywells_controls_compounds_qtys_total[i] > 0]);

% this list is only used for randomized strategy when num_plates_lines > 1
array[int] of int: emptywells_controls_compound_randomized =
      if (num_plates_lines = 1) \/ replicates_on_same_plate \/ replicates_on_different_plates then []
      else random_permutation(num_plates_lines*num_total_wells_line) endif;

% Divide materials between plates depending on the strategy
array[1..num_plates_lines, 1..len_emptywells_controls_compounds] of int: emptywells_controls_compounds_qtys =
     if num_plates_lines = 1
   then array2d(1..1, 1..len_emptywells_controls_compounds, emptywells_controls_compounds_qtys_total)
  elseif replicates_on_same_plate
      % currently this strategy is done in a simplified way - with an assumption that all compounds have equal total number of replicates
      % thus it's fairly trivial to split them between plates/lines.
      % If this assumption is broken, then bugs are possible and even expected
      % An alternative to somehow ameliorate this problem is to do a presolving a subproblem separately, but currently we dont do it
      % A python script and a separate model file should be easy enough to create
   then array2d(1..num_plates_lines, 1..len_emptywells_controls_compounds,
                [if i <= emptywells_controls_num
               then % this segment is an attempt to write a generalized formula that tries to spread out empty wells and controls
                    % more or less proportional to the available free space (i.e. everything except compounds)
                    % we check every column from emptywells_controls_compounds_first_non_zero (in case if there are no empty wells) to the last control
                    % first column, i.e. emptywells_controls_compounds_first_non_zero, is calculated by substracting quantities from every other columns on the same row (i.e. plate)
                    % every control column j uses floor(nb_replicates[j] * plate_space[i] / sum(plate_space) for plates i = 1..num_plates_lines-1.
                    % and plate i = num_plates_lines is the remaining number of nb_replicates[j].
                    %
                    % NOTE: in certain circumstances it can fail and produce negative qualities. e,g,:
                    % We have 5 plates with free spaces [2, 2, 2, 8, 8] (i.e. we have 22 free wells in total) and we need to place
                    % 2 empty wells and 20 controls. The current version produces the following distribution:
                    % [|  1,  1    % plate 1
                    %  |  1,  1    % plate 2
                    %  |  1,  1    % plate 3
                    %  |  1,  7    % plate 4
                    %  | -2, 10 |] % plate 5
                    % instead of, say:
                    % [|  0,  2    % plate 1
                    %  |  0,  2    % plate 2
                    %  |  0,  2    % plate 3
                    %  |  1,  7    % plate 4
                    %  |  1,  7 |] % plate 5
                    % There is no obvious way how to modify the current algorithm. The best course of action --
                    % set this pre-calculation as a separate constraint problem, solve it and use the result as an imput file.
                    % Likely as specialized Python script is required.
                    
                         % We calculate first column on the last plate:
                      if p = num_plates_lines /\ i = emptywells_controls_compounds_first_non_zero
                    then emptywells_controls_compounds_qtys_total[i] -
                         sum(p1 in 1..num_plates_lines-1)
                            (available_wells_per_plates_lines_total[p1] -
                             sum(i1 in emptywells_controls_compounds_first_non_zero + 1..emptywells_controls_num)
                                (floor((emptywells_controls_compounds_qtys_total[i1] * available_wells_per_plates_lines_total[p1]) /
                                       available_wells_per_plates_lines_total_sum)))
                         % We calculate remaining columns on the last plate:
                  elseif p = num_plates_lines /\ i > emptywells_controls_compounds_first_non_zero
                    then emptywells_controls_compounds_qtys_total[i] -
                         sum(p1 in 1..num_plates_lines-1)
                            (floor((emptywells_controls_compounds_qtys_total[i] * available_wells_per_plates_lines_total[p1]) /
                                   available_wells_per_plates_lines_total_sum))
                         % We calculate first column on every plate except the last one:
                  elseif i = emptywells_controls_compounds_first_non_zero /\ p < num_plates_lines
                    then available_wells_per_plates_lines_total[p] -
                         sum(i1 in emptywells_controls_compounds_first_non_zero + 1..emptywells_controls_num)
                            (floor((emptywells_controls_compounds_qtys_total[i1] * available_wells_per_plates_lines_total[p]) /
                             available_wells_per_plates_lines_total_sum))
                         % We calculate remaining columns on every plate except the last one:
                         % NOTE - we put this branch last in case if there's only one column to distribute,
                         % meaning this case will be covered by previous branches expicitly making
                         % sure that all allocated plate line space is used.
                         % If this branch is called earlier - it's not necessarily to be the case
                    else floor((emptywells_controls_compounds_qtys_total[i] * available_wells_per_plates_lines_total[p]) /
                               available_wells_per_plates_lines_total_sum)
                    endif
                 % and here the compounds are placed in accordance with the list emptywells_controls_compounds_names_concentrations
             elseif compounds_on_plates_lines[compound_index(compounds,
                                                             emptywells_controls_compounds_names_concentrations[1,i],
                                                             compound_names)] = p
               then emptywells_controls_compounds_qtys_total[i]
               else 0
              endif
                 | p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds]
                )
  elseif replicates_on_different_plates
      % strategy where everything is spread out more or less evenly between the plates
      % e.g. if we have emptywells_controls_compounds_qtys_total = [1, 3, 9, 7] corresponding to ['ew', 'ctr', 'cmp1', 'cmp2']
      % which we want to split between 4 plates. To do this we "create" vector:
      % ['ew',   'ctr',  'ctr',   'ctr', 'cmp1', 'cmp1', 'cmp1', 'cmp1', 'cmp1', 'cmp1',
      %  'cmp1', 'cmp1', 'cmp1', 'cmp2', 'cmp2', 'cmp2', 'cmp2', 'cmp2', 'cmp2', 'cmp2']
      % and then take 1st, 5th,  9th, etc elements for plate 1 (i.e. i mod 4 = 1)
      %               2nd, 6th, 10th, etc elements for plate 2 (i.e. i mod 4 = 2)
      %               3nd, 7th, 11th, etc elements for plate 3 (i.e. i mod 4 = 3)
      %               4nd, 8th, 12th, etc elements for plate 4 (i.e. i mod 4 = 0)
      % and aggregate the results to construct:
      % emptywells_controls_compounds_qtys = [| 0, 1, 2, 2
      %                                       | 1, 0, 3, 1
      %                                       | 0, 1, 2, 2
      %                                       | 0, 1, 2, 2
      %                                       |]
   then array2d(1..num_plates_lines, 1..len_emptywells_controls_compounds,
                [sum(j in sum(k in 1..i - 1)(emptywells_controls_compounds_qtys_total[k]) + 1
                            ..sum(k in 1..i)(emptywells_controls_compounds_qtys_total[k]))
                    (((j mod num_plates_lines) = p - 1))
                 | p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds]
                )

      % If both flages are set to false, we, currently, use randomized strategy
   else array2d(1..num_plates_lines, 1..len_emptywells_controls_compounds,
                [sum(j in sum(k in 1..i - 1)(emptywells_controls_compounds_qtys_total[k]) + 1
                            ..sum(k in 1..i)(emptywells_controls_compounds_qtys_total[k]))
                    (((emptywells_controls_compound_randomized[j] mod num_plates_lines) = p - 1))
                 | p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds]
                )
        endif;

% Check that material totals match (detects allocation failures when negatives were clamped to 0)
constraint assert(forall(p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds)
                        (emptywells_controls_compounds_qtys[p,i] >= 0),
				  "Model ERROR! Material distribution failed with 'replicates_on_same_plate' strategy.\n" ++
				      "The available space distribution is too imbalanced to allocate empty wells and controls.\n\n" ++
				      "Diagnostic Information:\n" ++
				      "  Number of plates: " ++ show(num_plates_lines) ++ "\n" ++
				      "  Available space per plate (after compounds): " ++ show(available_wells_per_plates_lines_total) ++ "\n" ++
				      "  Total available space: " ++ show(available_wells_per_plates_lines_total_sum) ++ "\n" ++
				      "  Materials to distribute:\n" ++
				      show(["    " ++ emptywells_controls_compounds_names_concentrations[1,i] ++ ": " ++ 
				            show(emptywells_controls_compounds_qtys_total[i]) ++ " required, " ++ 
				            show(sum(p in 1..num_plates_lines)(emptywells_controls_compounds_qtys[p,i])) ++ " allocated\n" 
				            | i in 1..emptywells_controls_num]) ++
				      "\nRECOMMENDATIONS:\n" ++
				      "  1. Use 'replicates_on_different_plates = true' instead (distributes materials evenly across plates)\n" ++
				      "  2. Manually balance compound distribution to create more uniform available spaces\n" ++
				      "  3. Use a Python preprocessing script to compute optimal empty wells/controls distribution\n"
				  );

array[int] of int: emptywells_controls_compounds_indices = [i | i in 1..len_emptywells_controls_compounds];
% emptywells_controls_compounds_indices - a short-hand index to make modeling easier

%%%%%%% Detecting some unsatisfiable cases & Data validation %%%%%%%%%%
% this one to be enabled later, if needed
%constraint assert(ceil(sum(compound_replicates)/numplates)*min(compound_concentrations++[infinity]) + sum([floor(control_concentrations[i]*control_replicates[i]/numplates) | i in 1..num_controls]) <= inner_plate_size, "Invalid data: the design is unsatisfiable. It is not possible to divide the compounds and controls evenly across the plates. (E01)");  
                                                                
% this one to be enabled later, if needed
%constraint assert((floor(sum(compound_replicates)/numplates)-1)*min(compound_concentrations++[infinity]) + max_compound_concentrations + sum([floor(control_concentrations[i]*control_replicates[i]/numplates) | i in 1..num_controls]) <= inner_plate_size, "Invalid data: the design is unsatisfiable. It is not possible to divide the compounds and controls evenly across the plates. (E02)");  
constraint forall(i in 1..num_plates_lines)
           (assert(sum(row(emptywells_controls_compounds_qtys,i)) = num_total_wells_line, "Model ERROR! Number of variables for wells in a line \(i) does not match!\n\(sum(row(emptywells_controls_compounds_qtys,i))) != \(num_total_wells_line)\n"));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Split the types of wells into criteria sets.
% A criteria set means that we would like to maximize distance between wells belonging to this set. The sets are:
% - empty wells only
% - control wells. If there are too many control wells (more than a half of designated line) then controls of different types are given their own set.
%   e.g. if we have too many positive and negative controls in total, then instead of trying to spread them as a whole (which would be impossible under given criteria) we will split them into two sets - only positive and only negative controls
% - compounds of the same type but with different concentrations. e.g. if we have "comp2" with concentrations "c21" and "c22" and "comp3" with concentrations "c31" and "c32" then we will have to sets: {"c21", "c22"} and {"c31", "c32"}
%   But if the total number of replicates of a given compound within a set is higher than the half of the designated plate space then we split each concentration into its own set, i.e. {"c21"}, {"c22"}, etc
array[1..num_plates_lines] of int: len_emptywells_controls_compounds_criteria_sets = 
 [length([1]
          ++
         if flag_consolidate_controls /\
            sum(i in 1..len_emptywells_controls_compounds where emptywells_controls_compounds_categories[i] = "ctrs")
               (emptywells_controls_compounds_qtys[p,i]) < num_total_wells_line div 2
         then  [1]
         else  [1 | i in 1..len_emptywells_controls_compounds where emptywells_controls_compounds_categories[i] = "ctrs"]
         endif
          ++
         [1
          | j in 1..compounds, k in 1..(if sum(l in 1..len_emptywells_controls_compounds
                                               where emptywells_controls_compounds_names_concentrations[1,l] =
                                                     compound_names[j])(emptywells_controls_compounds_qtys[p,l]) < num_total_wells_line div 2
                                      then 1
                                      else sum(l in 1..len_emptywells_controls_compounds
                                               where emptywells_controls_compounds_names_concentrations[1,l] = compound_names[j])(1)
                                           endif)])
         | p in 1..num_plates_lines];

int: num_criteria_sets = max(len_emptywells_controls_compounds_criteria_sets);

array[1..num_plates_lines, 1..num_criteria_sets] of set of int:
         emptywells_controls_compounds_criteria_sets =
           array2d(1..num_plates_lines, 1..num_criteria_sets,
                [([{i | i in 1..len_emptywells_controls_compounds where emptywells_controls_compounds_categories[i] = "ew"}]
                    ++
                   if flag_consolidate_controls /\
                      sum(i in 1..len_emptywells_controls_compounds where emptywells_controls_compounds_categories[i] = "ctrs")
                         (emptywells_controls_compounds_qtys[p,i]) < num_total_wells_line div 2
                   then  [{i | i in 1..len_emptywells_controls_compounds where emptywells_controls_compounds_categories[i] = "ctrs"}]
                   else  [{i} | i in 1..len_emptywells_controls_compounds where emptywells_controls_compounds_categories[i] = "ctrs"]
                   endif
                    ++
                   [if sum(l in 1..len_emptywells_controls_compounds where emptywells_controls_compounds_names_concentrations[1,l] =
                                                                           compound_names[j])(emptywells_controls_compounds_qtys[p,l]) < num_total_wells_line div 2 
                    then {i | i in 1..len_emptywells_controls_compounds where
                                      emptywells_controls_compounds_names_concentrations[1,i] = compound_names[j]}
                    else {min([l | l in 1..len_emptywells_controls_compounds where emptywells_controls_compounds_names_concentrations[1,l] =
                                                                           compound_names[j]]) + k - 1}
                   endif
                   | j in 1..compounds, k in 1..(if sum(l in 1..len_emptywells_controls_compounds
                                                        where emptywells_controls_compounds_names_concentrations[1,l] =
                                                           compound_names[j])(emptywells_controls_compounds_qtys[p,l]) < num_total_wells_line div 2
                                                           then 1
                                                 else sum(l in 1..len_emptywells_controls_compounds
                                                          where emptywells_controls_compounds_names_concentrations[1,l] =
                                                                      compound_names[j])(1)
                                                endif)]
                    ++
                   [{} | i in 1..(num_criteria_sets - len_emptywells_controls_compounds_criteria_sets[p])]
                   )[r]
                 | p in 1..num_plates_lines, r in 1..num_criteria_sets]);

array[1..num_plates_lines, 1..num_criteria_sets] of string:
         emptywells_controls_compounds_criteria_sets_categories =
           array2d(1..num_plates_lines, 1..num_criteria_sets,
                [(["ew"]
                    ++
                   if flag_consolidate_controls /\
                      sum(i in 1..len_emptywells_controls_compounds where emptywells_controls_compounds_categories[i] = "ctrs")
                         (emptywells_controls_compounds_qtys[p,i]) < num_total_wells_line div 2
                   then  ["ctrs"]
                   else  ["ctrs" | i in 1..len_emptywells_controls_compounds where emptywells_controls_compounds_categories[i] = "ctrs"]
                   endif
                    ++
                   ["comp"
                    | j in 1..compounds, k in 1..(if sum(l in 1..len_emptywells_controls_compounds
                                                         where emptywells_controls_compounds_names_concentrations[1,l] =
                                                               compound_names[j])(emptywells_controls_compounds_qtys[p,l]) < num_total_wells_line div 2
                                                then 1
                                                else sum(l in 1..len_emptywells_controls_compounds
                                                         where emptywells_controls_compounds_names_concentrations[1,l] = compound_names[j])(1)
                                                endif)]
                    ++
                   ["" | i in 1..(num_criteria_sets - len_emptywells_controls_compounds_criteria_sets[p])]
                   )[r]
                 | p in 1..num_plates_lines, r in 1..num_criteria_sets]);

array[1..num_plates_lines, 1..num_total_wells_line] of int: emptywells_controls_compounds_indices_ordered =
                             % Perform a test for non-negative values
                            if forall(p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds)
                                     (emptywells_controls_compounds_qtys[p,i] >= 0)
                             % If passed - calculate the matrix properly
                          then array2d(1..num_plates_lines, 1..num_total_wells_line,
                                       [emptywells_controls_compounds_indices[i] | p in 1..num_plates_lines,
                                                                                   i in 1..len_emptywells_controls_compounds,
                                                                                   j in 1..emptywells_controls_compounds_qtys[p,i]])
                             % If not passed, because of faulty pre-calculations, put dummy data and let an assert output the error
                          else array2d(1..num_plates_lines, 1..num_total_wells_line,
                                       [1 | p in 1..num_plates_lines, i in 1..num_total_wells_line])
                         endif;

% REMINDER: We assume that empty wells, controls and compounds are placed sequentially!
array[1..num_plates_lines, 1..len_emptywells_controls_compounds] of int:
     emptywells_controls_compounds_criteria_sets_transformed =
                               array2d(1..num_plates_lines, 1..len_emptywells_controls_compounds,
                                       [i | p in 1..num_plates_lines, i in 1..num_criteria_sets,
                                            j in emptywells_controls_compounds_criteria_sets[p,i]]);
array[1..num_plates_lines, 1..num_total_wells_line] of int: emptywells_controls_compounds_criteria_ordered = 
           array2d(1..num_plates_lines, 1..num_total_wells_line,
                   [emptywells_controls_compounds_criteria_sets_transformed[p,emptywells_controls_compounds_indices_ordered[p,i]]
                    | p in 1..num_plates_lines, i in 1..num_total_wells_line]);

% For quadrant distribution we use Deterministic round-robin allocation followed by sequential bin packing
% Imagine we have 5x5 grid. We have 4 materials (we call them 1, 2, 3 and 4, for simplicity) with quantities 10, 6, 4 and 5,
% that we will assign between 4 quadrants (I-IV). Since it's there are odd numbers of rows and columns,
% Q-I will take 9 replicates, Q-II and Q-III - 6 replicates each and Q-IV - 4 replicates.
% 
% So, what I do is to first list all the replicates in order:
% 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 4 4 4 4 4
% 
% I match them with mod 4 counts (to make them distinct I use i-iv):
% 
% 1  1   1  1 1  1   1  1 1  1   2  2 2  2   2  2 3  3   3  3 4  4   4  4 4
% i ii iii iv i ii iii iv i ii iii iv i ii iii iv i ii iii iv i ii iii iv i
% 
% I then list all the replicates associated with i, then with ii, iii and iv:
% 
% 1 1 1 2 3 4 4  1  1  1  2  3  4   1   1   2   2   3   4  1  1  2  2  3  4
% i i i i i i i ii ii ii ii ii ii iii iii iii iii iii iii iv iv iv iv iv iv
% 
% And lastly, I say that first 9 replicates in the list are Q-I, the next 6 - Q-II, etc:
% 
% 1 1 1 2 3 4 4  1  1  1  2  3  4   1   1   2   2   3   4   1   1  2  2  3  4
% i i i i i i i ii ii ii ii ii ii iii iii iii iii iii iii  iv  iv iv iv iv iv
% I I I I I I I  I  I II II II II  II  II III III III III III III IV IV IV IV
% 
% i.e:
% Q-I   will have 5, 1, 1, 2 replicates of materials 1-4 resp.
% Q-II  will have 3, 1, 1, 1 replicates of materials 1-4 resp.
% Q-III will have 2, 2, 1, 1 replicates of materials 1-4 resp.
% Q-IV  will have 0, 2, 1, 1 replicates of materials 1-4 resp.

int: base_quadrant_size = num_total_wells_line div 4;
int: remainder_wells = num_total_wells_line mod 4;
array[1..4] of int: quadrants_distribution = [base_quadrant_size + (remainder_wells >= i) | i in 1..4];
array[1..4] of int: quadrants_sizes = [floor(num_rows_line / 2) *  ceil(num_cols_line / 2) - num_forced_empty_wells_quadrant,
                                        ceil(num_rows_line / 2) *  ceil(num_cols_line / 2) - num_forced_empty_wells_quadrant,
                                       floor(num_rows_line / 2) * floor(num_cols_line / 2) - num_forced_empty_wells_quadrant,
                                        ceil(num_rows_line / 2) * floor(num_cols_line / 2) - num_forced_empty_wells_quadrant];

%%%%%%% Detecting some unsatisfiable cases & Data validation %%%%%%%%%%
% these two asserts should never fail, but they are kept here just in case
constraint assert(sum(quadrants_sizes) = num_total_wells_line, "Model ERROR! The total number of assigned wells per quadrant does not match:\nsum(\(emptywells_controls_compounds_qtys)) != \(num_total_wells_line)\n");
constraint assert(sum(quadrants_sizes) = sum(quadrants_distribution), "Model ERROR! The total number of assigned wells per quadrant does not match:\nsum(\(emptywells_controls_compounds_qtys)) != sum(\(quadrants_distribution))\n");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Splitting empty wells, controls, compounds by quadrants
% quadrant_shift_line rotates quadrant assignments across plates to avoid systematic bias
% Plate 1: Q-I first, Plate 2: Q-II first, Plate 3: Q-III first, Plate 4: Q-IV first, etc.
% PURPOSE: Prevents concentration from always being in same quadrant across plates
%          Distributes potential plate-specific effects more evenly
array[1..num_plates_lines] of int: quadrant_shift_line = [p - 1 | p in 1..num_plates_lines];

array[1..num_plates_lines, 1..num_total_wells_line] of 1..num_total_wells_line: quadrant_link_indices =
                          array2d(1..num_plates_lines, 1..num_total_wells_line,
                                  [([4 - ((3 + quadrant_shift_line[p]) mod 4) + 4 * (i - 1)
                                     | i in 1..quadrants_distribution[4 - ((3 + quadrant_shift_line[p]) mod 4)]] ++
                                    [4 - (     quadrant_shift_line[p]  mod 4) + 4 * (i - 1)
                                     | i in 1..quadrants_distribution[4 - (     quadrant_shift_line[p]  mod 4)]] ++
                                    [4 - ((1 + quadrant_shift_line[p]) mod 4) + 4 * (i - 1)
                                     | i in 1..quadrants_distribution[4 - ((1 + quadrant_shift_line[p]) mod 4)]] ++
                                    [4 - ((2 + quadrant_shift_line[p]) mod 4) + 4 * (i - 1)
                                     | i in 1..quadrants_distribution[4 - ((2 + quadrant_shift_line[p]) mod 4)]])[j]
                                   | p in 1..num_plates_lines, j in 1..num_total_wells_line]);

% When use_quadrant_distribution = false (i.e., only 1 row or 1 column), ALL wells go to quadrant 1
% This means symmetry breaking happens via global alldifferent instead of per-quadrant lex
bool: use_quadrant_distribution = (num_rows_line > 1) /\ (num_cols_line > 1);

array[1..num_plates_lines, 1..num_total_wells_line] of 1..4: quadrant_wells_assignment =
array2d(1..num_plates_lines, 1..num_total_wells_line,
    [if not use_quadrant_distribution
   then 1
 elseif has_element(i, [quadrant_link_indices[p,j] | j in 1..quadrants_sizes[1]]) 
   then 1
 elseif has_element(i, [quadrant_link_indices[p,j] | j in quadrants_sizes[1]..sum(k in 1..2)(quadrants_sizes[k])]) 
   then 2
 elseif has_element(i, [quadrant_link_indices[p,j] | j in sum(k in 1..2)(quadrants_sizes[k])+1..sum(k in 1..3)(quadrants_sizes[k])]) 
   then 3
   else 4 endif | p in 1..num_plates_lines, i in 1..num_total_wells_line]);

int: quadrant_1_rows = if use_quadrant_distribution then floor(num_rows_line / 2) else num_rows_line endif;
int: quadrant_1_cols = if use_quadrant_distribution then ceil(num_cols_line / 2) else num_cols_line endif;
                % if only 1 row/column per line no forced empty wells
constraint assert(not use_quadrant_distribution -> (num_corner_empty_wells = 0), "Invalid data: no forced empty corners when there is only 1 row and/or 1 column.");

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% OPTIONAL CONSTRAINT FLAGS %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% If concentrations_on_different_rows and/or concentrations_on_different_columns are true then:
% - we exclude compounds concentrations during application of the symmetry-breaking lex_chain constraint within its criteria set
% - we apply an all_different constraint on the relevant concentration later
%   to account for this we need to collect flags to denote which compound concentrations on which plate will be subjected to this treatment
% - we add an array - if the number of replicates of a concentration across all plates is less than the number of rows/columns
%   then the all_different constraint isn't just applied per single plate but across all plates (a stronger constraint)

array[int] of string: not_force_optional_ctrs = ["ew"];

% Changes from version 1.1:
% If a half (top/bottom/right/left) has leq nb of wells than its nb of rows/columns then apply all_different (same as before)
% If a half (top/bottom/right/left) has gt  nb of wells than its nb of rows/columns then apply gcc
%    to make sure that there are at least one well per row/column in this half.
%
% Thus we'll first split existing flags emptywells_controls_compounds_concentrations_on_different_* into:
%  - emptywells_controls_compounds_concentrations_on_different_rows_top     (top, all_different)
%  - emptywells_controls_compounds_concentrations_on_different_rows_top_gcc (top, gcc)
%  - emptywells_controls_compounds_concentrations_on_different_rows_btm     (bottom, all_different)
%  - emptywells_controls_compounds_concentrations_on_different_rows_btm_gcc (bottom, gcc)
%  - emptywells_controls_compounds_concentrations_on_different_columns_lft     (left, all_different)
%  - emptywells_controls_compounds_concentrations_on_different_columns_lft_gcc (left, gcc)
%  - emptywells_controls_compounds_concentrations_on_different_columns_rgt     (right, all_different)
%  - emptywells_controls_compounds_concentrations_on_different_columns_rgt_gcc (right, gcc)
%
% Same with flags emptywells_controls_compounds_concentrations_on_different_*_all_plates
% This behaviour should:
%  - reduce the number of situations when the model is unsatisfiable due to optional constraints
%  - ensure that every row/column is covered if there are too many wells with the criteria set (previously it was not enforced)
%  - the previous point, potentially, can lead to increased performance as the number of well placements is reduced (to be tested)
%
% Note that this requires tweaking the symmetry-breaking constraints to ensure that gcc flags are handled as well

% Helper functions to compute optional constraint flags:
function bool: optional_constraints_main_conditions(int: material, bool: global_flag) =
    % if there are only 1 row/column then we do not apply these constraints /\ the corresponding row/column flag must be enabled:
    use_quadrant_distribution /\ global_flag /\
	% no matter the flags, empty wells and controls are exempt from the treatment:
    not (emptywells_controls_compounds_categories[material] in not_force_optional_ctrs);
function bool: compute_optional_constraint_flag_alldifferent(int: plate, int: material, set of int: target_quadrants, bool: global_flag, int: threshold) =
    optional_constraints_main_conditions(material, global_flag) /\
	% must not exceed the threshold:
    sum(k in 1..num_total_wells_line where quadrant_wells_assignment[plate, k] in target_quadrants /\
                                           emptywells_controls_compounds_indices_ordered[plate, k] = material)(1) <= threshold;
function bool: compute_optional_constraint_flag_gcc(int: plate, int: material, set of int: target_quadrants, bool: global_flag, int: threshold) =
    optional_constraints_main_conditions(material, global_flag) /\
    % must exceed the threshold:
    sum(k in 1..num_total_wells_line where quadrant_wells_assignment[plate, k] in target_quadrants /\
                                           emptywells_controls_compounds_indices_ordered[plate, k] = material)(1) > threshold;
function bool: compute_optional_constraint_flag_allplates_alldifferent(int: material, set of int: target_quadrants, bool: global_flag, int: threshold) =
    % if there are only 1 row/column then we do not apply these constraints /\ the corresponding row/column flag must be enabled:
    use_quadrant_distribution /\ global_flag /\
	% no matter the flags, empty wells and controls are exempt from the treatment:
    not (emptywells_controls_compounds_categories[material] in not_force_optional_ctrs) /\
	% must not exceed the threshold:
    sum(p in 1..num_plates_lines, k in 1..num_total_wells_line where quadrant_wells_assignment[p, k] in target_quadrants /\
                                                                     emptywells_controls_compounds_indices_ordered[p, k] = material)(1) <= threshold;
function bool: compute_optional_constraint_flag_allplates_gcc(int: material, set of int: target_quadrants, bool: global_flag, int: threshold) =
    use_quadrant_distribution /\ global_flag /\
    not (emptywells_controls_compounds_categories[material] in not_force_optional_ctrs) /\
	% must exceed the threshold:
    sum(p in 1..num_plates_lines, k in 1..num_total_wells_line where quadrant_wells_assignment[p, k] in target_quadrants /\
                                                                     emptywells_controls_compounds_indices_ordered[p, k] = material)(1) > threshold;

% Flags for optional constraints (all_different):
array[1..num_plates_lines, 1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_rows_top =
       array2d(1..num_plates_lines, 1..len_emptywells_controls_compounds,
               [compute_optional_constraint_flag_alldifferent(p, i, {1,3}, concentrations_on_different_rows, floor(num_rows_line / 2))
                | p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds]);

array[1..num_plates_lines, 1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_rows_btm =
       array2d(1..num_plates_lines, 1..len_emptywells_controls_compounds,
               [compute_optional_constraint_flag_alldifferent(p, i, {2,4}, concentrations_on_different_rows, ceil(num_rows_line / 2))
                | p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds]);

array[1..num_plates_lines, 1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_columns_lft =
       array2d(1..num_plates_lines, 1..len_emptywells_controls_compounds,
               [compute_optional_constraint_flag_alldifferent(p, i, {1,2}, concentrations_on_different_columns, ceil(num_cols_line / 2))
                | p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds]);

array[1..num_plates_lines, 1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_columns_rgt =
       array2d(1..num_plates_lines, 1..len_emptywells_controls_compounds,
               [compute_optional_constraint_flag_alldifferent(p, i, {3,4}, concentrations_on_different_columns, floor(num_cols_line / 2))
                | p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds]);

% Flags for optional constraints (GCC):
array[1..num_plates_lines, 1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_rows_top_gcc =
       array2d(1..num_plates_lines, 1..len_emptywells_controls_compounds,
               [compute_optional_constraint_flag_gcc(p, i, {1,3}, concentrations_on_different_rows, floor(num_rows_line / 2))
                | p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds]);

array[1..num_plates_lines, 1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_rows_btm_gcc =
       array2d(1..num_plates_lines, 1..len_emptywells_controls_compounds,
               [compute_optional_constraint_flag_gcc(p, i, {2,4}, concentrations_on_different_rows, ceil(num_rows_line / 2))
                | p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds]);

array[1..num_plates_lines, 1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_columns_lft_gcc =
       array2d(1..num_plates_lines, 1..len_emptywells_controls_compounds,
               [compute_optional_constraint_flag_gcc(p, i, {1,2}, concentrations_on_different_columns, ceil(num_cols_line / 2))
                | p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds]);

array[1..num_plates_lines, 1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_columns_rgt_gcc =
       array2d(1..num_plates_lines, 1..len_emptywells_controls_compounds,
               [compute_optional_constraint_flag_gcc(p, i, {3,4}, concentrations_on_different_columns, floor(num_cols_line / 2))
                | p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds]);

% Across all plates flags (summing across all plates, alldifferent):
array[1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_rows_all_plates_top =
               [compute_optional_constraint_flag_allplates_alldifferent(i, {1,3}, concentrations_on_different_rows, floor(num_rows_line / 2))
                | i in 1..len_emptywells_controls_compounds];

array[1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_rows_all_plates_btm =
               [compute_optional_constraint_flag_allplates_alldifferent(i, {2,4}, concentrations_on_different_rows, ceil(num_rows_line / 2))
                | i in 1..len_emptywells_controls_compounds];

array[1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_columns_all_plates_lft =
               [compute_optional_constraint_flag_allplates_alldifferent(i, {1,2}, concentrations_on_different_columns, ceil(num_cols_line / 2))
                | i in 1..len_emptywells_controls_compounds];

array[1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_columns_all_plates_rgt =
               [compute_optional_constraint_flag_allplates_alldifferent(i, {3,4}, concentrations_on_different_columns, floor(num_cols_line / 2))
                | i in 1..len_emptywells_controls_compounds];

% Across all plates flags (gcc variant)
array[1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_rows_all_plates_top_gcc =
               [compute_optional_constraint_flag_allplates_gcc(i, {1,3}, concentrations_on_different_rows, floor(num_rows_line / 2))
                | i in 1..len_emptywells_controls_compounds];

array[1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_rows_all_plates_btm_gcc =
               [compute_optional_constraint_flag_allplates_gcc(i, {2,4}, concentrations_on_different_rows, ceil(num_rows_line / 2))
                | i in 1..len_emptywells_controls_compounds];

array[1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_columns_all_plates_lft_gcc =
               [compute_optional_constraint_flag_allplates_gcc(i, {1,2}, concentrations_on_different_columns, ceil(num_cols_line / 2))
                | i in 1..len_emptywells_controls_compounds];

array[1..len_emptywells_controls_compounds] of bool:
      emptywells_controls_compounds_concentrations_on_different_columns_all_plates_rgt_gcc =
               [compute_optional_constraint_flag_allplates_gcc(i, {3,4}, concentrations_on_different_columns, floor(num_cols_line / 2))
                | i in 1..len_emptywells_controls_compounds];


% Calculate flags for each critera set - use this set in the optimization criteria, use fake edges when calculating the criteria
% TODO: See if it is possible to rewrite it later into a smarter criteria which counts "white" and "black" squares and tests is their number is equal
array[1..num_plates_lines, 1..num_criteria_sets] of int: num_criteria_set_line =
                                    array2d(1..num_plates_lines, 1..num_criteria_sets,
                                            [sum(i in emptywells_controls_compounds_criteria_sets[p,j])
                                                (emptywells_controls_compounds_qtys[p,i])
                                             | p in 1..num_plates_lines, j in 1..num_criteria_sets]);
array[1..num_plates_lines, 1..num_criteria_sets] of bool: use_min_dist_criteria =
                                    array2d(1..num_plates_lines, 1..num_criteria_sets,
                                            [% Do not apply if the set is empty or contains only 1 well:
											 num_criteria_set_line[p,j] > 1 /\
											 % Do not apply if the set has too many wells:
                                             num_criteria_set_line[p,j] <= min(num_wells_set_criteria,
											 								   num_threshold_min_dist_criteria * num_total_wells_line)
                                             | p in 1..num_plates_lines, j in 1..num_criteria_sets]);

% - case 1: when the number of wells covers less than allotted amount of the plate, then       use fake edge wells
% - case 2: when the number of wells covers more than allotted amount of the plate, then don't use fake edge wells
array[1..num_plates_lines, 1..num_criteria_sets] of bool: use_fake_edge_wells =
                                  array2d(1..num_plates_lines, 1..num_criteria_sets,
                                          [num_criteria_set_line[p,j] > 0 /\
                                           emptywells_controls_compounds_criteria_sets_categories[p,j] != "ctrs" /\
                                           emptywells_controls_compounds_criteria_sets_categories[p,j] != "ew" /\
                                           num_criteria_set_line[p,j] <= num_threshold_use_fake_edge_wells
                                           | p in 1..num_plates_lines, j in 1..num_criteria_sets]);
int: num_edges = 16;
% Fake edges are always generated
array[1..num_edges,1..2] of int: edges = array2d(1..num_edges,1..2,
                                                 [0,                            0,
                                                  num_rows_line + 1,            0,
                                                  0,                            num_cols_line + 1,
                                                  num_rows_line + 1,            num_cols_line + 1,
                                                  
                                                  num_rows_line + 2,            floor(1 * num_cols_line / 4),
                                                  num_rows_line + 2,            floor(2 * num_cols_line / 4),
                                                  num_rows_line + 2,            floor(3 * num_cols_line / 4),
                                                  
                                                  -1,                           floor(1 * num_cols_line / 4),
                                                  -1,                           floor(2 * num_cols_line / 4),
                                                  -1,                           floor(3 * num_cols_line / 4),
                                                  
                                                  floor(1 * num_rows_line / 4), num_cols_line + 2,
                                                  floor(2 * num_rows_line / 4), num_cols_line + 2,
                                                  floor(3 * num_rows_line / 4), num_cols_line + 2,
                                                  
                                                  floor(1 * num_rows_line / 4), -1,
                                                  floor(2 * num_rows_line / 4), -1,
                                                  floor(3 * num_rows_line / 4), -1
                                                  ]);
int: min_dist_edges = max(1, min([num_rows_line, num_cols_line]) div 2);
                                      
% Calculate distance matrix edges and wells
int: num_edges_wells = num_edges + num_total_wells_line;

int: len_min_dist_criteria = sum(p in 1..num_plates_lines)(length([1 | i in 1..num_criteria_sets where use_min_dist_criteria[p,i]]));

% Counting tensors - count number of elements per plate, quadrant, criteria set / concentration
% Some are used for symmetry breaking constraints and some for the calculation of the distance criteria
function bool: well_excluded_by_optional_constraints(int: quadrant, int: plate, int: material) =
    (quadrant in [1,3] -> (not emptywells_controls_compounds_concentrations_on_different_rows_top[plate,material] /\
                           not emptywells_controls_compounds_concentrations_on_different_rows_top_gcc[plate,material])) /\
    (quadrant in [2,4] -> (not emptywells_controls_compounds_concentrations_on_different_rows_btm[plate,material] /\
                           not emptywells_controls_compounds_concentrations_on_different_rows_btm_gcc[plate,material])) /\
    (quadrant in [1,2] -> (not emptywells_controls_compounds_concentrations_on_different_columns_lft[plate,material] /\
                           not emptywells_controls_compounds_concentrations_on_different_columns_lft_gcc[plate,material])) /\
    (quadrant in [3,4] -> (not emptywells_controls_compounds_concentrations_on_different_columns_rgt[plate,material] /\
                           not emptywells_controls_compounds_concentrations_on_different_columns_rgt_gcc[plate,material]));

array[1..num_plates_lines, 1..4, 1..num_criteria_sets] of int: quadrant_criteria_counts =
    array3d(1..num_plates_lines, 1..4, 1..num_criteria_sets,
            [ sum([1 | k in 1..num_total_wells_line where emptywells_controls_compounds_criteria_ordered[p,k] = j /\
                                                          quadrant_wells_assignment[p,k] = i /\
                                                          % Changed in 1.2: the splitting of the flags necessitated the change in how to account for them
                                                          well_excluded_by_optional_constraints(i,p,emptywells_controls_compounds_indices_ordered[p,k])
                   ])
             | p in 1..num_plates_lines, i in 1..4, j in 1..num_criteria_sets]);
array[1..num_plates_lines, 1..4, 1..len_emptywells_controls_compounds] of int: quadrant_on_different_rows_columns =
    array3d(1..num_plates_lines, 1..4, 1..len_emptywells_controls_compounds,
            [ sum([1 | k in 1..num_total_wells_line where emptywells_controls_compounds_indices_ordered[p,k] = j /\
                                                          quadrant_wells_assignment[p,k] = i /\
                                                          % Changed in 1.2: the splitting of the flags necessitated the change in how to use them
                                                          not well_excluded_by_optional_constraints(i,p,j)
                   ])
             | p in 1..num_plates_lines, i in 1..4, j in 1..len_emptywells_controls_compounds]);

array[1..len_min_dist_criteria] of int: distances_edges_controls_compounds_plates_match =
                                        [p | p in 1..num_plates_lines, i in 1..num_criteria_sets where use_min_dist_criteria[p,i]];
array[1..len_min_dist_criteria] of int: distances_edges_controls_compounds_sets_match =
                                        [i | p in 1..num_plates_lines, i in 1..num_criteria_sets where use_min_dist_criteria[p,i]];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%-------------------CONSTRAINT MODEL--------------------------
% PLACEMENT VARIABLES:
array[1..num_plates_lines, 1..num_total_wells_line, 1..2] of var 1..max(num_rows_line, num_cols_line): emptywells_controls_compounds_coordinates::no_output;
array[1..num_plates_lines, 1..num_total_wells_line] of var 1..num_rows_line*num_cols_line: wells_line_lex_all_different_substitute::is_defined_var::no_output;

% if we have too many empty wells, so much so that we don't consider them to be a part of the criteria function then we set empty wells coordinates to a default value
predicate set_default_empty_well_coords(int: quadrant, int: default_row, int: default_col) =
    forall(p in 1..num_plates_lines, i in 1..num_total_wells_line 
           where quadrant_wells_assignment[p,i] = quadrant /\
                 emptywells_controls_compounds_indices_ordered[p,i] = 1 /\
                 not use_min_dist_criteria[p,emptywells_controls_compounds_criteria_ordered[p,i]])
          (emptywells_controls_compounds_coordinates[p, i, 1] = default_row /\
           emptywells_controls_compounds_coordinates[p, i, 2] = default_col);

constraint set_default_empty_well_coords(1, quadrant_1_rows, quadrant_1_cols);
constraint set_default_empty_well_coords(2, floor(num_rows_line / 2) + 1, ceil(num_cols_line / 2));
constraint set_default_empty_well_coords(3, floor(num_rows_line / 2), ceil(num_cols_line / 2) + 1);
constraint set_default_empty_well_coords(4, floor(num_rows_line / 2) + 1, ceil(num_cols_line / 2) + 1);

% domain restrictions on wells of each quadrant
predicate set_quadrant_domain(int: quadrant, int: dim, set of int: domain) =
    forall(p in 1..num_plates_lines, i in 1..num_total_wells_line where quadrant_wells_assignment[p,i] = quadrant)
          (emptywells_controls_compounds_coordinates[p, i, dim] in domain);
predicate enforce_corner_exclusion(int: quadrant, 
                                   int: row_threshold, bool: row_less_than,
                                   int: col_threshold, bool: col_less_than) =
    forall(p in 1..num_plates_lines, i in 1..num_total_wells_line 
           where quadrant_wells_assignment[p,i] = quadrant /\ num_corner_empty_wells > 0)
          ((if row_less_than 
            then emptywells_controls_compounds_coordinates[p, i, 1] <= row_threshold
            else emptywells_controls_compounds_coordinates[p, i, 1] >= row_threshold endif)
           ->
           (if col_less_than 
            then emptywells_controls_compounds_coordinates[p, i, 2] <= col_threshold
            else emptywells_controls_compounds_coordinates[p, i, 2] >= col_threshold endif));

constraint set_quadrant_domain(1, 1, 1..quadrant_1_rows);
constraint set_quadrant_domain(1, 2, 1..quadrant_1_cols);
constraint enforce_corner_exclusion(1, num_corner_empty_wells, true, num_corner_empty_wells + 1, false);

constraint set_quadrant_domain(2, 1, floor(num_rows_line / 2) + 1..num_rows_line);
constraint set_quadrant_domain(2, 2, 1..ceil(num_cols_line / 2));
constraint enforce_corner_exclusion(2, num_rows_line - num_corner_empty_wells + 1, false, num_corner_empty_wells + 1, false);

constraint set_quadrant_domain(3, 1, 1..floor(num_rows_line / 2));
constraint set_quadrant_domain(3, 2, ceil(num_cols_line / 2) + 1..num_cols_line);
constraint enforce_corner_exclusion(3, num_corner_empty_wells, true, num_cols_line - num_corner_empty_wells, true);

constraint set_quadrant_domain(4, 1, floor(num_rows_line / 2) + 1..num_rows_line);
constraint set_quadrant_domain(4, 2, ceil(num_cols_line / 2) + 1..num_cols_line);
constraint enforce_corner_exclusion(4, num_rows_line - num_corner_empty_wells + 1, false, num_cols_line - num_corner_empty_wells, true);


% a symmetry breaking ctr - with the caveat that it isn't the best one I can think of. It has a drawback (somewhat grouping one type of concentrations together, in some cases, - it is fine for now, but later I must think on how to restructure it
constraint forall(p in 1..num_plates_lines, i in 1..4, j in 1..num_criteria_sets where (j = 1 -> use_min_dist_criteria[p,j]) /\
                                                                                       quadrant_criteria_counts[p,i,j] >= 2) % Only for 2+ wells
  (% we first shuffle the list of replicates, to ensure that specific concentrations are not forced to be on certain ranges of rows
   let {
    % 1. Collect wells in original order
    array[int] of int: wells_original = 
      [k | k in 1..num_total_wells_line where emptywells_controls_compounds_criteria_ordered[p,k] = j /\
                                              quadrant_wells_assignment[p,k] = i /\
                                              well_excluded_by_optional_constraints(i,p,emptywells_controls_compounds_indices_ordered[p,k]) /\
                                              (j = 1 -> use_min_dist_criteria[p,j])], % exlude empty wells if not a part of the criteria],

    % 2. Generate random permutation
    array[int] of int: perm = random_permutation(length(wells_original)),
    
    % 3. Reorder wells
    array[int] of int: wells_permuted = 
      [wells_original[perm[idx]] | idx in 1..length(wells_original)]
  } in
  lex_chain_greater(array2d(1..2, 1..length(wells_permuted),
                            [emptywells_controls_compounds_coordinates[p, wells_permuted[idx], dim] 
                             | dim in 1..2, idx in 1..length(wells_permuted)])));


% symmetry breaking ctr - applied if concentrations_on_different_rows or concentrations_on_different_columns
constraint forall(p in 1..num_plates_lines, i in 1..4, j in 1..len_emptywells_controls_compounds where quadrant_on_different_rows_columns[p,i,j] >= 2)
                                                                                                       % Only for 2+ wells 
           (lex_chain_greater(array2d(1..2, 1..quadrant_on_different_rows_columns[p,i,j],
                                      [emptywells_controls_compounds_coordinates[p,k,l] | l in 1..2, k in 1..num_total_wells_line
                                       where emptywells_controls_compounds_indices_ordered[p,k] = j /\
                                             quadrant_wells_assignment[p,k] = i /\
                                             not well_excluded_by_optional_constraints(i,p,j)
                                       ])));

% Optional ctrs: apply all_different/gcc constraints when concentrations_on_different_rows and/or concentrations_on_different_columns WITHIN a plate
predicate apply_alldifferent_constraint(array[int, int] of bool: flag_array, int: coordinate_dim, set of int: target_quadrants) =
    forall(p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds where flag_array[p,i])
          (all_different([emptywells_controls_compounds_coordinates[p,j,coordinate_dim] | j in 1..num_total_wells_line
                          where emptywells_controls_compounds_indices_ordered[p,j] = i /\ quadrant_wells_assignment[p,j] in target_quadrants]));
% Applies Global Cardinality Constraint with lower bound = 1, upper bound = parameter
% DESIGN RATIONALE for upper bounds:
%   - Single-plate GCC: upper_bound = num_cols (or num_rows)
%         Each row can have at most num_cols wells ✓ (physical constraint)
%   - All-plates GCC: upper_bound = num_cols * num_plates
%         Across plates, each row position can have up to (num_cols × plates) wells
%         This is LOOSE (theoretically all could be on plate 1), but the critical constraint
%          is the LOWER BOUND = 1, ensuring every row/col is covered at least once
%         Tighter upper bound would be more complex and not practically beneficial

predicate apply_gcc_constraint(array[int, int] of bool: flag_array, int: coordinate_dim, set of int: target_quadrants, set of int: value_range, int: upper_bound) =
    forall(p in 1..num_plates_lines, i in 1..len_emptywells_controls_compounds where flag_array[p,i])
          (global_cardinality([emptywells_controls_compounds_coordinates[p,j,coordinate_dim] | j in 1..num_total_wells_line
                               where emptywells_controls_compounds_indices_ordered[p,j] = i /\ quadrant_wells_assignment[p,j] in target_quadrants],
                              [j | j in value_range],
                              [1 | j in value_range],
                              [upper_bound | j in value_range]));

constraint apply_alldifferent_constraint(emptywells_controls_compounds_concentrations_on_different_rows_top, 1, {1,3});
constraint apply_alldifferent_constraint(emptywells_controls_compounds_concentrations_on_different_rows_btm, 1, {2,4});
constraint apply_alldifferent_constraint(emptywells_controls_compounds_concentrations_on_different_columns_lft, 2, {1,2});
constraint apply_alldifferent_constraint(emptywells_controls_compounds_concentrations_on_different_columns_rgt, 2, {3,4});

constraint apply_gcc_constraint(emptywells_controls_compounds_concentrations_on_different_rows_top_gcc, 1, {1,3}, 1..floor(num_rows_line / 2), num_cols_line);
constraint apply_gcc_constraint(emptywells_controls_compounds_concentrations_on_different_rows_btm_gcc, 1, {2,4}, floor(num_rows_line / 2) + 1..num_rows_line, num_cols_line);
constraint apply_gcc_constraint(emptywells_controls_compounds_concentrations_on_different_columns_lft_gcc, 2, {1,2}, 1..ceil(num_cols_line / 2), num_rows_line);
constraint apply_gcc_constraint(emptywells_controls_compounds_concentrations_on_different_columns_rgt_gcc, 2, {3,4}, ceil(num_cols_line / 2) + 1..num_cols_line, num_rows_line);

% Optional ctrs: apply all_different/gcc constraints when concentrations_on_different_rows and/or concentrations_on_different_columns ACROSS all plates
predicate apply_alldifferent_constraint_allplates(array[int] of bool: flag_array, int: coordinate_dim, set of int: target_quadrants) =
    forall(i in 1..len_emptywells_controls_compounds where flag_array[i])
          (all_different([emptywells_controls_compounds_coordinates[p,j,coordinate_dim] | p in 1..num_plates_lines, j in 1..num_total_wells_line
                          where emptywells_controls_compounds_indices_ordered[p,j] = i /\ quadrant_wells_assignment[p,j] in target_quadrants]));
predicate apply_gcc_constraint_allplates(array[int] of bool: flag_array, int: coordinate_dim, set of int: target_quadrants, set of int: value_range, int: upper_bound) =
    forall(i in 1..len_emptywells_controls_compounds where flag_array[i])
          (global_cardinality([emptywells_controls_compounds_coordinates[p,j,coordinate_dim] | p in 1..num_plates_lines, j in 1..num_total_wells_line
                               where emptywells_controls_compounds_indices_ordered[p,j] = i /\ quadrant_wells_assignment[p,j] in target_quadrants],
                              [j | j in value_range],
                              [1 | j in value_range],
                              [upper_bound | j in value_range]));

constraint apply_alldifferent_constraint_allplates(emptywells_controls_compounds_concentrations_on_different_rows_all_plates_top, 1, {1,3});
constraint apply_alldifferent_constraint_allplates(emptywells_controls_compounds_concentrations_on_different_rows_all_plates_btm, 1, {2,4});
constraint apply_alldifferent_constraint_allplates(emptywells_controls_compounds_concentrations_on_different_columns_all_plates_lft, 2, {1,2});
constraint apply_alldifferent_constraint_allplates(emptywells_controls_compounds_concentrations_on_different_columns_all_plates_rgt, 2, {3,4});

constraint apply_gcc_constraint_allplates(emptywells_controls_compounds_concentrations_on_different_rows_all_plates_top_gcc,
                                          1, {1,3}, 1..floor(num_rows_line / 2), num_cols_line * num_plates_lines);
constraint apply_gcc_constraint_allplates(emptywells_controls_compounds_concentrations_on_different_rows_all_plates_btm_gcc,
                                          1, {2,4}, floor(num_rows_line / 2) + 1..num_rows_line, num_cols_line * num_plates_lines);
constraint apply_gcc_constraint_allplates(emptywells_controls_compounds_concentrations_on_different_columns_all_plates_lft_gcc,
                                          2, {1,2}, 1..ceil(num_cols_line / 2), num_rows_line * num_plates_lines);
constraint apply_gcc_constraint_allplates(emptywells_controls_compounds_concentrations_on_different_columns_all_plates_rgt_gcc,
                                          2, {3,4}, ceil(num_cols_line / 2) + 1..num_cols_line, num_rows_line * num_plates_lines);


% making sure that we don't place two or more materials in a same well with a lex_alldifferent constraint substitution placed upon replicate coordinates
constraint forall(p in 1..num_plates_lines, i in 1..num_total_wells_line)
                 (wells_line_lex_all_different_substitute[p,i] = 
                  (emptywells_controls_compounds_coordinates[p,i,1] - 1) * num_cols_line + emptywells_controls_compounds_coordinates[p,i,2]);

constraint forall(p in 1..num_plates_lines, q in 1..4)
           (all_different([wells_line_lex_all_different_substitute[p,i] | i in 1..num_total_wells_line where quadrant_wells_assignment[p,i] = q
                                                                           % exlude empty wells if not a part of the criteria
                                                                          /\ (emptywells_controls_compounds_indices_ordered[p,i] = 1 ->
                                                                              use_min_dist_criteria[p,emptywells_controls_compounds_criteria_ordered[p,i]])]));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% CALCULATING THE OPTIMIZATION CRITERIA - distance matrix, minimum distances within sets
% Notes:
% - distances between fake edges are not calculated
% - distances between a fake edge and a well are calculated if the well is under the group with true use_min_dist_criteria and use_fake_edge_wells
% - distances between two wells are calculated if they are only under the same group AND the group is under use_min_dist_criteria flag

array[1..len_min_dist_criteria] of var 1..max(2, num_rows_line + num_cols_line - 2): distances_edges_controls_compounds_minimum;
set of int: ctrs_criteria = {j | j in 1..len_min_dist_criteria where emptywells_controls_compounds_criteria_sets_categories[distances_edges_controls_compounds_plates_match[j],
                                                                                                                            distances_edges_controls_compounds_sets_match[j]] = "ctrs"};
var 1..max(2, num_rows_line + num_cols_line - 2): min_distance;
var 1..max(2, num_rows_line + num_cols_line - 2): min_distance_ctrs;

% Define min_distance as the minimum across all criteria sets
constraint if len_min_dist_criteria > 0 then
    % All individual criteria distances must be >= overall minimum
    forall(i in 1..len_min_dist_criteria)(
        distances_edges_controls_compounds_minimum[i] >= min_distance
    )
else
    % No criteria sets to optimize - set default value
    min_distance = 1
endif;

% Define min_distance_ctrs as the minimum across control criteria sets only
constraint if card(ctrs_criteria) > 0 then
    % All control criteria distances must be >= control-specific minimum
    forall(j in ctrs_criteria)(
        distances_edges_controls_compounds_minimum[j] >= min_distance_ctrs
    )
else
    % No control criteria sets - set default value
    min_distance_ctrs = 1
endif;
% Restricting minimum distances within sets
constraint forall(i in 1..len_min_dist_criteria
				          where use_min_dist_criteria[distances_edges_controls_compounds_plates_match[i],distances_edges_controls_compounds_sets_match[i]] /\
		        		          use_fake_edge_wells[  distances_edges_controls_compounds_plates_match[i],distances_edges_controls_compounds_sets_match[i]])
				        (forall(k in 1..num_total_wells_line, j in 1..num_edges where
				                (emptywells_controls_compounds_criteria_ordered[distances_edges_controls_compounds_plates_match[i],k] = distances_edges_controls_compounds_sets_match[i]))
                        (distance_manhattan(row(edges, j),
                                            [emptywells_controls_compounds_coordinates[distances_edges_controls_compounds_plates_match[i],k,1],
                                             emptywells_controls_compounds_coordinates[distances_edges_controls_compounds_plates_match[i],k,2]]) >= distances_edges_controls_compounds_minimum[i]
			      ));
constraint forall(i in 1..len_min_dist_criteria
				          where use_min_dist_criteria[distances_edges_controls_compounds_plates_match[i],distances_edges_controls_compounds_sets_match[i]])
				        (forall(j in 1..num_total_wells_line - 1, k in j + 1..num_total_wells_line where
				                (distances_edges_controls_compounds_sets_match[i] = emptywells_controls_compounds_criteria_ordered[distances_edges_controls_compounds_plates_match[i],j]) /\
				                (distances_edges_controls_compounds_sets_match[i] = emptywells_controls_compounds_criteria_ordered[distances_edges_controls_compounds_plates_match[i],k]))
                       (distance_manhattan([emptywells_controls_compounds_coordinates[distances_edges_controls_compounds_plates_match[i],j,1],
                                            emptywells_controls_compounds_coordinates[distances_edges_controls_compounds_plates_match[i],j,2]],
                                           [emptywells_controls_compounds_coordinates[distances_edges_controls_compounds_plates_match[i],k,1],
                                            emptywells_controls_compounds_coordinates[distances_edges_controls_compounds_plates_match[i],k,2]]) >= distances_edges_controls_compounds_minimum[i]
                  ));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

solve minimize - min_distance * 1000 - min_distance_ctrs * 100 - sum(distances_edges_controls_compounds_minimum);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Debugging output %%%
output [if debugging then "[num_rows, num_cols, num_corner_empty_wells, horizontal_cell_lines, vertical_cell_lines] = \([num_rows, num_cols, num_corner_empty_wells, horizontal_cell_lines, vertical_cell_lines])\n\n" else "" endif];

output [if debugging then "[num_rows_line, num_cols_line, num_total_wells_line, num_total_wells, num_plates_lines, num_empty_wells] = \([num_rows_line, num_cols_line, num_total_wells_line, num_total_wells, num_plates_lines, num_empty_wells])\n\n" else "" endif];

output [if debugging then "len_emptywells_controls_compounds = \(len_emptywells_controls_compounds)\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_categories =\n" else "" endif];
output [if debugging then show(emptywells_controls_compounds_categories) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_names_concentrations =\n" else "" endif];
output [if debugging then show2d(emptywells_controls_compounds_names_concentrations) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "[replicates_on_different_plates, replicates_on_same_plate] = \([replicates_on_different_plates, replicates_on_same_plate])\n\n" else "" endif];

output [if debugging then "[compounds_per_plate_min, compounds_per_plate_mod, emptywells_controls_num] = \([compounds_per_plate_min, compounds_per_plate_mod, emptywells_controls_num])\n\n" else "" endif];

output [if debugging then "compounds_on_plates_lines =\n" else "" endif];
output [if debugging then show(compounds_on_plates_lines) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "compounds_per_plates_lines_total =\n" else "" endif];
output [if debugging then show(compounds_per_plates_lines_total) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "available_wells_per_plates_lines_total =\n" else "" endif];
output [if debugging then show(available_wells_per_plates_lines_total) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_qtys_total =\n" else "" endif];
output [if debugging then show(emptywells_controls_compounds_qtys_total) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_qtys =\n" else "" endif];
output [if debugging then show2d(emptywells_controls_compounds_qtys) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_qtys (per replicate, for verification) =\n" else "" endif];
output [if debugging then show([sum(p in 1..num_plates_lines)(emptywells_controls_compounds_qtys[p,i])
        | i in 1..len_emptywells_controls_compounds]) else "" endif];
output [if debugging then "\n\n" else "" endif];
        
output [if debugging then "emptywells_controls_compounds_qtys (per plate line, for verification) =\n" else "" endif];
output [if debugging then show([num_total_wells_line] ++ [sum(i in 1..len_emptywells_controls_compounds)(emptywells_controls_compounds_qtys[p,i])
              | p in 1..num_plates_lines]) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_rows_top =\n" else "" endif];
output [if debugging then show2d(emptywells_controls_compounds_concentrations_on_different_rows_top) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_rows_top_gcc =\n" else "" endif];
output [if debugging then show2d(emptywells_controls_compounds_concentrations_on_different_rows_top_gcc) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_rows_btm =\n" else "" endif];
output [if debugging then show2d(emptywells_controls_compounds_concentrations_on_different_rows_btm) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_rows_btm_gcc =\n" else "" endif];
output [if debugging then show2d(emptywells_controls_compounds_concentrations_on_different_rows_btm_gcc) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_columns_lft =\n" else "" endif];
output [if debugging then show2d(emptywells_controls_compounds_concentrations_on_different_columns_lft) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_columns_lft_gcc =\n" else "" endif];
output [if debugging then show2d(emptywells_controls_compounds_concentrations_on_different_columns_lft_gcc) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_columns_rgt =\n" else "" endif];
output [if debugging then show2d(emptywells_controls_compounds_concentrations_on_different_columns_rgt) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_columns_rgt_gcc =\n" else "" endif];
output [if debugging then show2d(emptywells_controls_compounds_concentrations_on_different_columns_rgt_gcc) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_rows_all_plates_top_gcc =\n" else "" endif];
output [if debugging then show(emptywells_controls_compounds_concentrations_on_different_rows_all_plates_top_gcc) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_rows_all_plates_top_gcc =\n" else "" endif];
output [if debugging then show(emptywells_controls_compounds_concentrations_on_different_rows_all_plates_top_gcc) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_rows_all_plates_btm =\n" else "" endif];
output [if debugging then show(emptywells_controls_compounds_concentrations_on_different_rows_all_plates_btm) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_rows_all_plates_btm_gcc =\n" else "" endif];
output [if debugging then show(emptywells_controls_compounds_concentrations_on_different_rows_all_plates_btm_gcc) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_columns_all_plates_lft =\n" else "" endif];
output [if debugging then show(emptywells_controls_compounds_concentrations_on_different_columns_all_plates_lft) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_columns_all_plates_lft_gcc =\n" else "" endif];
output [if debugging then show(emptywells_controls_compounds_concentrations_on_different_columns_all_plates_lft_gcc) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_columns_all_plates_rgt =\n" else "" endif];
output [if debugging then show(emptywells_controls_compounds_concentrations_on_different_columns_all_plates_rgt) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_concentrations_on_different_columns_all_plates_rgt_gcc =\n" else "" endif];
output [if debugging then show(emptywells_controls_compounds_concentrations_on_different_columns_all_plates_rgt_gcc) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "num_criteria_set_line =\n" else "" endif];
output [if debugging then show2d(num_criteria_set_line) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "[num_criteria_sets, len_min_dist_criteria, num_edges_wells] = \([num_criteria_sets, len_min_dist_criteria, num_edges_wells])\n\n" else "" endif];

output [if debugging then "emptywells_controls_compounds_criteria_sets =\n" else "" endif];
output [if debugging then show2d(emptywells_controls_compounds_criteria_sets) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "use_min_dist_criteria =\n" else "" endif];
output [if debugging then show2d(use_min_dist_criteria) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "use_fake_edge_wells =\n" else "" endif];
output [if debugging then show2d(use_fake_edge_wells) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "quadrant_on_different_rows_columns =\n" else "" endif];
output [if debugging then show3d(quadrant_on_different_rows_columns) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "distances_edges_controls_compounds_plates_match =\n" else "" endif];
output [if debugging then show(distances_edges_controls_compounds_plates_match) else "" endif];
output [if debugging then "\n\n" else "" endif];

output [if debugging then "distances_edges_controls_compounds_sets_match =\n" else "" endif];
output [if debugging then show(distances_edges_controls_compounds_sets_match) else "" endif];
output [if debugging then "\n\n" else "" endif];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% CSV Ouput %%%

array[int] of string: letters_capital = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S", "T", "U", "V", "W", "X", "Y", "Z"];
array[int] of string: letters_inline  = ["a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l", "m", "n", "o", "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z"];
array[int] of string: letters = letters_capital ++ [letters_capital[i] ++ letters_inline[j] | i in 1..length(letters_capital), j in 1..length(letters_inline)];

output [if not debugging then "plateID,well,cmpdname,CONCuM,cmpdnum,VOLuL\n" else "" endif];

output [if not debugging then "plate_\(p),"
 ++ letters[fix(calculate_line_index(inner_empty_edge, size_empty_edge, num_rows_line, h, emptywells_controls_compounds_coordinates[p,k,1]))]
 ++ if fix(calculate_line_index(inner_empty_edge, size_empty_edge, num_cols_line, v, emptywells_controls_compounds_coordinates[p,k,2])) < 10
      then "0" else "" endif
 ++ "\(calculate_line_index(inner_empty_edge, size_empty_edge, num_cols_line, v, emptywells_controls_compounds_coordinates[p,k,2]))"
 ++ ","
 ++ emptywells_controls_compounds_names_concentrations[1,emptywells_controls_compounds_indices_ordered[p,k]]
 ++ ","
 ++ emptywells_controls_compounds_names_concentrations[2,emptywells_controls_compounds_indices_ordered[p,k]]
 ++ ","
 ++ emptywells_controls_compounds_names_concentrations[1,emptywells_controls_compounds_indices_ordered[p,k]]
 ++ "_"
 ++ emptywells_controls_compounds_names_concentrations[2,emptywells_controls_compounds_indices_ordered[p,k]]
 ++ "\n" else "" endif | p in 1..num_plates_lines, k in 1..num_total_wells_line,
                         h in 1..horizontal_cell_lines, v in 1..vertical_cell_lines
                         where fix(emptywells_controls_compounds_indices_ordered[p,k]) > 1];
                             % fix(emptywells_controls_compounds_indices_ordered[p,k]) = 1]; % output only empty wells


output ["criteria function = "];
output ["\(- min_distance * 1000 - min_distance_ctrs * 100 - sum(distances_edges_controls_compounds_minimum))\n"];
output [if debugging then show([min_distance] ++ distances_edges_controls_compounds_minimum) else "" endif];
output [if debugging then "\n" else "" endif];