From 3aa6f3238811690b9c2568bf6900c3b50a3de213 Mon Sep 17 00:00:00 2001 From: d33bs Date: Tue, 27 Oct 2026 12:56:07 -0600 Subject: [PATCH] add plaid platemapping skill --- .agents/skills/plaid-platemap/SKILL.md | 155 ++ .../skills/plaid-platemap/assets/README.md | 11 + .../assets/layout_predicates.mzn | 62 + .../plaid-platemap/assets/plate-design.mzn | 2016 +++++++++++++++++ .../plaid-platemap/scripts/plaid_platemap.py | 844 +++++++ .../tests/test_plaid_platemap.py | 302 +++ 6 files changed, 3390 insertions(+) create mode 100644 .agents/skills/plaid-platemap/SKILL.md create mode 100644 .agents/skills/plaid-platemap/assets/README.md create mode 100644 .agents/skills/plaid-platemap/assets/layout_predicates.mzn create mode 100644 .agents/skills/plaid-platemap/assets/plate-design.mzn create mode 100644 .agents/skills/plaid-platemap/scripts/plaid_platemap.py create mode 100644 .agents/skills/plaid-platemap/tests/test_plaid_platemap.py diff --git a/.agents/skills/plaid-platemap/SKILL.md b/.agents/skills/plaid-platemap/SKILL.md new file mode 100644 index 0000000..c32f9b7 --- /dev/null +++ b/.agents/skills/plaid-platemap/SKILL.md @@ -0,0 +1,155 @@ +--- +name: plaid-platemap +description: Use when designing or improving plate layouts with plaid (MiniZinc). +version: 0.1.0 +author: Dave Bunten (d33bs), Hermes Agent +license: MIT +platforms: [macos, linux] +metadata: + hermes: + tags: [plaid, minizinc, platemap, plate-layout, screening] + related_skills: [cytotable, petalibrary] +--- + +# Plaid Platemap Skill + +Improves hit-validation plate layouts with plaid (pharmbio/plaid, constraint +programming in MiniZinc) to reduce well-related cross-effects: same-compound +wells touching each other, edge clustering, and unbalanced row/column spread. +The companion script handles the non-standard workbook format of the Way Lab +hit-validation platemap (three sheets, fixed two-zone cell plating) that plaid +cannot consume directly, and translates plaid's solution back into the same +Excel format for collaborators. + +## When to Use + +- Asked to design, randomize, or improve a microplate layout, especially when + worried about proximity effects between wells. +- Asked to use plaid (or MiniZinc/Gecode) for plate design on this project. +- Reviewing a hit-validation workbook ("Hit Validation ... for Way Lab.xlsx") + for cross-effect weaknesses before an experiment is run. +- Don't use for analysis of already-run experiments (cytotable territory), for + Alpine/Slurm submission, or when the layout has no fixed-zone constraint and + plaid's web GUI alone is enough. + +## Prerequisites + +- MiniZinc with the Gecode solver: `terminal(command="brew install minizinc")` + on macOS (`apt`/`dnf` equivalents on Linux). Verify with + `terminal(command="minizinc --version")`; Gecode ships in the bundle. +- `uv` for the companion script (PEP 723 inline metadata pulls openpyxl): + `terminal(command="uv --version")`. +- The plaid model files are vendored in `assets/` (plate-design.mzn and + layout_predicates.mzn, pinned at pharmbio/plaid commit ecf4cfe). To refresh, + clone the plaid repo and diff; do not silently upgrade the model. +- The workbook under test: 'Platemap grid' is the authoritative per-well table + (WellRow/WellCol/well_position/heart_number/cell_type/heart_failure_type + plus one treatment_plate N column per plate); 'Visual Platemaps' holds the + rosters ('Randomized order' table) and sketches. Sketches in V3 are stylized + and contradict the grid - the grid wins, and apply corrects the sketches. + +## How to Run + +All commands from the repo root. The pipeline is five subcommands of one +script; each consumes the previous output. + +```sh +# 1. workbook -> canonical layout JSON (validates the V3 structure) +uv run --script .agents/skills/plaid-platemap/scripts/plaid_platemap.py extract \ + --input "path/to/Hit Validation V3 for Way Lab.xlsx" \ + --output layout.json + +# 2. layout JSON -> plaid .dzn (compounds-only model on the used plate area) +uv run --script .agents/skills/plaid-platemap/scripts/plaid_platemap.py dzn \ + --layout layout.json --output plate.dzn + +# 3. solve with MiniZinc/Gecode until a seed passes verification +uv run --script .agents/skills/plaid-platemap/scripts/plaid_platemap.py solve \ + --mzn .agents/skills/plaid-platemap/assets/plate-design.mzn \ + --dzn plate.dzn --layout layout.json --output solution.json + +# 4. write the redesigned workbook (original untouched) + QC sheet +uv run --script .agents/skills/plaid-platemap/scripts/plaid_platemap.py apply \ + --input "path/to/Hit Validation V3 for Way Lab.xlsx" \ + --output "path/to/Hit Validation V3 for Way Lab - plaid.xlsx" \ + --layout layout.json --solution solution.json + +# 5. independently re-verify the output workbook (exit non-zero on failure) +uv run --script .agents/skills/plaid-platemap/scripts/plaid_platemap.py verify \ + --input "path/to/Hit Validation V3 for Way Lab - plaid.xlsx" \ + --layout layout.json --reference layout.json +``` + +Put layout.json, plate.dzn, solution.json, and the output workbook in the +task's scratch/work directory (not the repo) unless asked to commit them. + +## Procedure + +1. **Extract and validate.** Run `extract`. It asserts the used region is the + 96-well inner area (rows B-G, cols 2-11), compounds sit only on failing + wells, per-plate compound rosters match the grid, and all plates share one + compound well pattern. Completion: extract prints zone sizes and the + before-metrics (V3: 17 same-compound adjacent pairs per plate). +2. **Generate the dzn.** Run `dzn`. V3 translates to 3 compounds x 8 replicates + on an 8x12 plate with a 1-well empty edge (inner area = the 60 used wells). + Completion: dzn written with the roster dose as the single concentration. +3. **Solve.** Run `solve`. It tries seeds 1..25, relocates plaid wells that + land on nonfailing wells onto free failing wells (king-distance >= 2 from + the compound's other wells, same row preferred), and checks the final + pattern: 24 compound wells, 8 per compound, all inside the failing zone, + same-compound king distance >= 2. First passing seed wins. + Completion: solution.json written; PASS line printed. +4. **Apply.** Run `apply` to write the new workbook: updated 'Platemap grid' + treatments, 'Visual Platemaps' sketches/per-well block corrected to match + the grid, plus a 'Platemap QC' sheet with before/after metrics. The + original file is never modified. Completion: output workbook path printed. +5. **Verify.** Run `verify` on the output workbook. It re-reads the grid + independently and checks every invariant per plate (rosters, replicate + counts, zones, adjacency, shared pattern) and prints before/after metrics. + Completion: "VERIFY OK" with exit code 0; any violation exits non-zero. + +## Pitfalls + +- **Color coding must follow values.** 'Visual Platemaps' color-codes by value + (V3: Failing 25 / Nonfailing 7 / Cpd 1-3 / DMSO fills). apply learns the + legend from the source workbook at extract time and recolors every cell it + rewrites; verify re-checks each sketch cell's fill against the legend and + fails on any mismatch ("color check: N mismatches" plus details on stderr). +- **DMSO cannot be a plaid control here.** 36 DMSO/60 wells exceeds the max + independent set (~30), so plaid's same-kind control separation would be + UNSAT. DMSO is the complement of the compound pattern instead. +- **Labels-only reshuffling cannot fix V3.** The V3 compound band + (C10-C11, D-E, F02-F03) cannot host 8 king-separated replicates per + compound (max ~6 in that shape); the well pattern itself must spread over + the full failing region, so apply also updates the Visual Platemaps + sketches to the new reality. +- **The V3 sketches are stale.** Hearts/compound sketches and the Platemap 4/5 + headers disagree with the grid; the grid is authoritative, apply rewrites + the sketches, and verify only trusts the grid. +- **plaid geometry vs physical plate.** The model uses the full 8x12 plate + with size_empty_edge=1, so inner coordinates map identically to physical + wells B-G x 2-11 - no row/col arithmetic in the relocation step. +- **Random seed matters.** Gecode returns the first satisfying layout for a + seed; different seeds give different layouts. The solve loop accepts the + first verified seed and records it in solution.json and the QC sheet. +- **MiniZinc timeouts.** `solve --timeout` bounds each seed (default 60 s); + raise it for bigger compound sets or denser plates. + +## Verification + +- `verify` passes on the output workbook (exit 0) with min same-compound + distance >= 2 on every plate and a clean color check (0 mismatches). +- 'Platemap grid' still has 60 wells; failing/nonfailing zones and heart + numbers unchanged; every plate keeps its 3 compounds x 8 replicates + 36 + DMSO (18 failing + 18 nonfailing). +- 'Visual Platemaps' fills match the learned legend exactly: sketch cells + carry the fill of their current value, not the fill of the previous layout's + occupant (a flipped well changes color with its content). +- 'Platemap QC' sheet records the before/after numbers (V3: 17 adjacent + same-compound pairs per plate -> 0). +- The original workbook is byte-identical after the run (apply refuses + --output equal to --input). +- Tests: `uv run --no-project --with pytest --with openpyxl pytest + .agents/skills/plaid-platemap/tests/test_plaid_platemap.py` exercises + extract/dzn/solve/apply/verify round trips on a small synthetic workbook + fixture (the solve test skips automatically when minizinc is missing). diff --git a/.agents/skills/plaid-platemap/assets/README.md b/.agents/skills/plaid-platemap/assets/README.md new file mode 100644 index 0000000..ae2be04 --- /dev/null +++ b/.agents/skills/plaid-platemap/assets/README.md @@ -0,0 +1,11 @@ +# Vendored plaid model files + +- `plate-design.mzn`, `layout_predicates.mzn`: copied verbatim from + at commit + `ecf4cfe6d8029dfadfceb79f5d6f2a1814d792c5` (Apache-2.0, see the LICENSE at + ). +- Do not edit these files; if you upgrade the pin, re-copy both files from + the new commit and note the new hash in `scripts/plaid_platemap.py` + (`PLAID_COMMIT`) and in the skill's Prerequisites section. +- The original repo's own examples (regression tests, JSON examples) are not + vendored; only the two model files the solver needs. diff --git a/.agents/skills/plaid-platemap/assets/layout_predicates.mzn b/.agents/skills/plaid-platemap/assets/layout_predicates.mzn new file mode 100644 index 0000000..f3e7fdf --- /dev/null +++ b/.agents/skills/plaid-platemap/assets/layout_predicates.mzn @@ -0,0 +1,62 @@ +include "globals.mzn"; + + +predicate sorted_except_0(array[int] of var int: x, int: num, int: min, int: max) = + regular(x, concat(["(0* \(i)){\(min),\(max)} " | i in 1..num]++["0*"]))::domain; + + +%% Cost implied constraints, similar to the ones derived for a group automaton constraint +predicate my_implied_cost(var int:cost, array[int] of var int: x) = + let { + int: L = length(x); + array[1..L] of var int: c_cost; + } in + c_cost[1] = x[1] /\ + c_cost[2] <= 1 /\ + c_cost[L] = cost /\ + forall(i in 1..L-1)(x[i+1] = c_cost[i+1] - c_cost[i] ) /\ + forall(i in 1..L-1)(c_cost[i+1] = x[i+1] + c_cost[i] ) /\ + forall(i in 1..L-1)(c_cost[i] <= c_cost[i+1]) /\ + forall(i in 2..L-1)(c_cost[i+1] <= 1 + c_cost[i-1]) /\ + forall(i in 2..L-1)(c_cost[i+1] - 1 <= c_cost[i-1]); + + +%% Cost implied constraints, similar to the ones derived for a group automaton constraint +predicate my_implied_cost_ext(var int:cost, array[int] of var int: x) = + let { + int: L = length(x); + array[1..L] of var int: c_cost; + } in + c_cost[1] = x[1] /\ + c_cost[2] <= 1 /\ + c_cost[L] = cost /\ + forall(i in 1..L-1)(x[i+1] = c_cost[i+1] - c_cost[i] ) /\ + forall(i in 1..L-1)(c_cost[i+1] = x[i+1] + c_cost[i] ) /\ + forall(i in 1..L-1)(c_cost[i] <= c_cost[i+1]) /\ + forall(i in 2..L-2)(c_cost[i+2] <= 1 + c_cost[i-1]) /\ + forall(i in 2..L-2)(c_cost[i+2] - 1 <= c_cost[i-1]); + + +function bool: force_spread_controls(int: numcols, int: numrows, int: total_controls, int: numplates) = + % 96-well plate (This is a proven upperbound) + if (numcols = 12 /\ numrows = 8) \/ (numcols = 8 /\ numrows = 12) then (11*numplates >= total_controls) + % 96-well plate with empty border (This is a proven upperbound) + elseif (numcols = 10 /\ numrows = 6) \/ (numcols = 6 /\ numrows = 10) then (7*numplates >= total_controls) + % 384-well plate with empty border (I have not proven yet that 33 is unsat) + % 32 is solvable, but makes the model too slow. + % TODO: Think about an intermediate constraint + elseif (numcols = 22 /\ numrows = 14) \/ (numcols = 14 /\ numrows = 22) then (31*numplates >= total_controls) + % 384-well plate (I have not proven yet that 38 is unsat) + elseif (numcols = 24 /\ numrows = 16) \/ (numcols = 16 /\ numrows = 24) then (37*numplates >= total_controls) + % Some generic attempt that + elseif numcols*numrows < 60 then (ceil(numcols*numrows/9)*numplates >= total_controls) + else false + endif; + + + +%force_spread_controls = if inner_plate_size <= 60 then (ceil(inner_plate_size/9)*numplates >= total_controls) else (floor(inner_plate_size/9)*numplates >= total_controls) endif; + +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; diff --git a/.agents/skills/plaid-platemap/assets/plate-design.mzn b/.agents/skills/plaid-platemap/assets/plate-design.mzn new file mode 100644 index 0000000..7c932d3 --- /dev/null +++ b/.agents/skills/plaid-platemap/assets/plate-design.mzn @@ -0,0 +1,2016 @@ +% Copyright 2021 PLAID 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 highly effective microplate layout designer +% +% Authors: Maria Andreina FRANCISCO RODRIGUEZ (maria.andreina.francisco@farmbio.uu.se) +% Version: 1.0 +% Last Revision: September 2022 +% + +include "globals.mzn"; +include "gecode.mzn"; +include "layout_predicates.mzn"; +%include "chuffed.mzn"; + + +%%%% 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; +bool: force_spread_controls; +bool: force_spread_concentrations; + +%% 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; +int: replicates = max(compound_replicates++[0]); +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) %% +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]); + +%% Potentially new parameters %% +bool: balance_controls_inside_plate = true; +bool: interconnected_plates = true; %(numplates==1); %Makes the problem much harder by connecting the plates +int: control_slack = 0; %Makes the problem harder/easier by constraining how equally should we distribute controls +bool: spread_controls = (ceil(inner_plate_size/2)*numplates >= total_controls); %Makes the problem harder when there are many controls + +%% Which individual controls can be spread more than others? +array[1..num_controls*max_control_concentrations] of bool: spread_control = [ ceil(inner_plate_size/4)*numplates >= (control_replicates[floor((i-1)/max_control_concentrations)+1]*(((i-1) mod max_control_concentrations)= 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(compounds==0 \/ min(compound_replicates) > 0,"Invalid datafile: Number of replicates should be larger than zero."); +constraint assert(compounds==0 \/ min(compound_concentrations) > 0,"Invalid datafile: Number of concentrations 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"); +constraint assert(not (num_controls == 1 /\ spread_control[1]==false),"Invalid datafile: There are too many controls of only one kind. This is not allowed at the moment. If you believe this is a mistake, please contact the developers."); +constraint assert(num_rows_line mod 2 == 0, "Invalid datafile: Currently we only support plate sizes that have an inner area with an even number of rows."); +constraint assert(num_cols_line mod 2 == 0, "Invalid datafile: Currently we only support plate sizes that have an inner area with an even number of columns." ); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +%% Number of wells needed. Note that plates might not be full +int: total_wells = sum([compound_concentrations[i]*compound_replicates[i] | i in 1..compounds]) + total_controls; + +set of int: Rows = 1..num_rows_line; +set of int: Columns = 1..num_cols_line; +set of int: Plates = 1..numplates; + +enum Vertical = {upper,lower}; +enum Horizontal = {left,right}; + +array [Plates,Rows] of var 0..num_cols_line: experiments_in_plate_row; +array [Plates,Columns] of var 0..num_rows_line: experiments_in_plate_column; + +bool: even_columns = ((num_cols_line mod 2) == 0); +bool: even_rows = ((num_rows_line mod 2) == 0); + +int: inner_plate_size = num_cols_line*num_rows_line; + + +%%%%% Data validation %%%%% +constraint assert(total_wells > 0, "Invalid data: the plates cannot be completely empty."); +constraint assert(inner_plate_size>0, "Invalid data: There are no wells on the plate."); +constraint assert(min(compound_concentrations++[0]) <= inner_plate_size, "Invalid data: Number of concentrations does not fit in one plate. If you think this is a mistake, please contact the development team."); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +%% Number of plates needed +%% max is used to avoid division-by-zero errors +int: numplates = max(ceil(total_wells/inner_plate_size),1); + +% TODO: temporal fix. This makes the problem too difficult for the borderline case in large plates. +% Probably because it tries to prove if that optimal solution exists. +%force_spread_controls = (ceil(inner_plate_size/9)*numplates >= total_controls); + +force_spread_controls = force_spread_controls(num_cols_line, num_rows_line, total_controls, numplates); + + +%force_spread_concentrations = if replicates_on_same_plate \/ replicates <= numplates then (ceil(inner_plate_size/9) >= replicates*max_compound_concentrations) else (ceil(inner_plate_size/9) >= max_compound_concentrations) endif ; + + +% TODO: There are many other cases here that could be considered, for example, individual compounds with fewer concentrations/replicates. +force_spread_concentrations = if replicates_on_different_plates /\ replicates <= numplates then ((inner_plate_size div 5) >= max_compound_concentrations) else (ceil(inner_plate_size/5) >= replicates*max_compound_concentrations) endif ; + + +%% All types of experiments (excluding controls): compounds and combinations +int: experiments = compounds*max_compound_concentrations*replicates; + +%% Used in redundant constraints +int: emptywells = numplates*inner_plate_size - total_wells; + + + +%%%%%%% Detecting some unsatisfiable cases & Data validation %%%%%%%%%% +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)"); + +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 assert(emptywells == 0 \/ allow_empty_wells,"Warning: There are \(emptywells) empty wells. You can allow the plates to have inner empty wells or modify other quantities to fill them up."); + +constraint assert(emptywells >= 0,"Model ERROR! Inner empty wells is negative. This should never happen!"); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%% Plates (Main model & Solution) %%%%%%%%%%% +array [Plates,Rows,Columns] of var 0..(experiments+num_controls*max_control_concentrations): plates; + +% 0 = empty wells +% 1... compounds*compound_concentrations*replicates = compounds +% compounds*compound_concentrations*replicates+1 ... experiments = combinations +% experiments+1... = controls -> control1_conc1, control1_conc2, ...control2_conc1, control2_conc2... + +%% Alternative (extra) model +array [1..experiments] of var {0} union Plates: experiment_plate; + +array [1..experiments] of var {0} union Rows: experiment_row; + +array [1..experiments] of var {0} union Columns: experiment_column; + + + + + + +%%%% CONSTRAINTS %%%% +%% Place only the exact amount of compounds, controls, and empty wells. +constraint global_cardinality(plates,[i | i in 0..(experiments + num_controls*max_control_concentrations)],[emptywells]++[(floor((i-1)/(compounds*max_compound_concentrations)) (experiment_plate[i] = 0 /\ experiment_row[i] = 0 /\ experiment_column[i] = 0)); + +constraint forall(i in 1..experiments)(((floor((i-1)/(compounds*max_compound_concentrations)) (experiment_plate[i] != 0 /\ experiment_row[i] != 0 /\ experiment_column[i] != 0)); + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Channelling constraints: an experiment/compound is located at a given plate + +constraint forall(l in 1..experiments, i in Plates)(count_eq(plates[i,..,..],l,(experiment_plate[l] == i))); + +constraint forall(l in 1..experiments)(count_eq(plates,l,(experiment_plate[l] != 0))); + + +%% For experiments that do not exist = plate 0 (redundant) +%constraint forall(l in 1..experiments)( + % count_eq(array1d(1..numplates*inner_plate_size, plates),l,0) <-> (experiment_plate[l] = 0 /\ experiment_row[l] = 0 /\ experiment_column[l] = 0) +%); + +%constraint forall(l in 1..experiments)(count_eq(array1d(1..numplates*inner_plate_size, plates),l,0) <-> (experiment_row[l] == 0)); +%constraint forall(l in 1..experiments)(count_eq(array1d(1..numplates*inner_plate_size, plates),l,0) <-> (experiment_column[l] == 0)); + + + +%% Channeling constraints between the two models + +constraint forall(l in 1..experiments, j in Rows)(count_eq(plates[..,j,..],l,(experiment_row[l] == j))); +constraint forall(l in 1..experiments, k in Columns)(count_eq(plates[..,..,k],l,(experiment_column[l] == k))); + +constraint forall(l in 1..experiments)(count_eq(plates,l,(experiment_row[l] != 0))); +constraint forall(l in 1..experiments)(count_eq(plates,l,(experiment_column[l] != 0))); + + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% All concentrations of a given compound replica must appear on a single plate + +constraint forall(l in 1..experiments where ((l mod max_compound_concentrations) == 1))(all_equal([experiment_plate[i] | i in l..(l+compound_concentrations[((((l-1) div max_compound_concentrations)) mod compounds)+1]-1)])); + +%% TODO: check if other compounds have fewer concentrations and can be spread around more +%% TODO: change to include replicates +%constraint if force_spread_concentrations then forall(i in Plates, j in 1..num_rows_line-1, k in 1..num_cols_line-1,l in 1..experiments where ((l mod max_compound_concentrations) == 1)) (1>=among([plates[i,j,k],plates[i,j+1,k],plates[i,j,k+1],plates[i,j+1,k+1]],(l..(l+compound_concentrations[(floor(((l-1)/max_compound_concentrations)) mod compounds)+1]-1)))):: domain endif; + + +% constraint if force_spread_concentrations then forall(i in Plates, j in 1..num_rows_line-1, k in 1..num_cols_line-1,l in 1..experiments where ((l mod max_compound_concentrations) == 1)) ( +% let{ +% var 0..1: num_concentrations; +% } in among(num_concentrations,[plates[i,j,k],plates[i,j+1,k],plates[i,j,k+1],plates[i,j+1,k+1]],(l..(l+compound_concentrations[(floor(((l-1)/max_compound_concentrations)) mod compounds)+1]-1)))::domain) endif; + +% Replaced by constraint below for clarity (less index math) +%constraint if force_spread_concentrations then forall(i in Plates, j in 1..num_rows_line-1, k in 1..num_cols_line-1,l in 1..experiments where ((l mod max_compound_concentrations) == 1)) ( +%2 > among([plates[i,j,k],plates[i,j+1,k],plates[i,j,k+1],plates[i,j+1,k+1]],(l..(l+compound_concentrations[(floor(((l-1)/max_compound_concentrations)) mod compounds)+1]-1)))::domain) endif; + + +% Replaced by constraint below to in order to include replicates +%constraint if force_spread_concentrations then forall(i in Plates, j in 1..num_rows_line-1, k in 1..num_cols_line-1, c in 1..compounds, r in 0..replicates-1) ( +%2 > among([plates[i,j,k],plates[i,j+1,k],plates[i,j,k+1],plates[i,j+1,k+1]],((c-1)*r*max_compound_concentrations+1..c*r*max_compound_concentrations))::domain) endif; + + +constraint if force_spread_concentrations then forall(i in Plates, j in 1..num_rows_line-1, k in 1..num_cols_line-1, c in 1..compounds) ( +2 > among([plates[i,j,k],plates[i,j+1,k],plates[i,j,k+1],plates[i,j+1,k+1]], +array_union([(c-1)*max_compound_concentrations+1 + r*compounds*max_compound_concentrations..c*max_compound_concentrations + r*compounds*max_compound_concentrations | r in 0..replicates-1]) +)::domain) endif; + + +%%array [1..experiments] of var {0} union Plates: experiment_plate; array_union + +array [1..replicates] of set of int: tt = [r*max_compound_concentrations+1..r*max_compound_concentrations | r in 0..replicates-1] ; + +%compound_replicates[(floor(((l-1)/max_compound_concentrations)) mod compounds)+1] + + +% constraint forall(i in Plates, j in Rows) ( +% let{ +% var int: exprInRow; +% constraint exprInRow = among(plates[i,j,..], 1..experiments)::domain; +% %sum([experiment_row[l] == j /\ experiment_plate[l] == i| l in 1..experiments]); +% } in exprInRow <= num_cols_line-among([plates[i,j,k] | k in Columns], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))); + + +constraint forall(i in Plates, j in Rows) ( +(among(plates[i,j,..], 1..experiments)::domain) <= num_cols_line-among([plates[i,j,k] | k in Columns], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Spreading concentrations of a compound across different rows and columns +%%% (consider balancing the number of rows and columns!) +% Consider making it different on different plates + + + +% Might be redundant != implied +constraint forall(i in Plates, j in Rows, k in Columns) ((0 < plates[i,j,k] /\ plates[i,j,k] <= experiments) <-> experiment_plate[plates[i,j,k]]==i); +constraint forall(i in Plates, j in Rows, k in Columns) ((0 < plates[i,j,k] /\ plates[i,j,k] <= experiments) <-> experiment_row[plates[i,j,k]]==j); +constraint forall(i in Plates, j in Rows, k in Columns) ((0 < plates[i,j,k] /\ plates[i,j,k] <= experiments) <-> experiment_column[plates[i,j,k]]==k); + + +% Spread the compounds across different rows and columns + +%% Spreads compounds over different rows. It's been replaced by the more restrictive alldifferent version below. +%constraint if concentrations_on_different_rows then (forall(l in 1..compounds) (forall(r in 0..(compound_replicates[l]-1)) (nvalue(min(compound_concentrations[l],num_rows_line-2*size_empty_edge), [experiment_row[r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i]| i in 1..compound_concentrations[l]])))) endif; + +constraint if concentrations_on_different_rows then (forall(l in 1..compounds) (forall(r in 0..(compound_replicates[l]-1)) ( +alldifferent([experiment_row[r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i]| i in 1..min(compound_concentrations[l],num_rows_line)]):: domain +))) endif; + +constraint if concentrations_on_different_rows then (forall(l in 1..compounds) (forall(r in 0..(compound_replicates[l]-1)) ( +alldifferent([experiment_row[r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i]| i in 1..min(compound_concentrations[l],num_rows_line)]):: domain +))) endif; + +%% HERE!!! NOW! +%% +%constraint if force_spread_concentrations then forall(i in Plates, j in 1..num_rows_line-1, k in 1..num_cols_line-1, c in 1..compounds) ( +%2 > among([plates[i,j,k],plates[i,j+1,k],plates[i,j,k+1],plates[i,j+1,k+1]], +%array_union([(c-1)*max_compound_concentrations+1 + r*compounds*max_compound_concentrations..c*max_compound_concentrations + r*compounds*max_compound_concentrations | r in 0..replicates-1]) +%)::domain) endif; + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% WORK AREA!!! +%balanced_alldifferent(array[int] of var int: x, int: mid_val, int: mid_count) + +%array [Vertical, 1..experiments] of var {0} union Rows: ul_experiment_row; + +%constraint forall(v in Vertical, e in 1..experiments) (if experiment_row[e] <= (num_rows_line div 2) then ul_experiment_row[upper,e]==experiment_row[e] /\ ul_experiment_row[lower,e]==0 else ul_experiment_row[lower,e]==experiment_row[e] /\ ul_experiment_row[upper,e]==0 endif); + +% constraint if concentrations_on_different_rows then (forall(l in 1..compounds) (if compound_concentrations[l]>1 then +% forall(v in Vertical, r in 0..(compound_replicates[l]-1)) ( +% let{ +% var int: compounds_in_row; +% constraint compounds_in_row <= ((compound_concentrations[l]+1) div 2)+1; +% constraint compounds_in_row >= (compound_concentrations[l] div 2)+1; +% } +% in nvalue(compounds_in_row,[ul_experiment_row[v,r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i]| i in 1..min(compound_concentrations[l],num_rows_line)]) +% )endif)) endif; + + + +%constraint if concentrations_on_different_rows then (forall(l in 1..compounds) (forall(r in 0..(compound_replicates[l]-1)) ( +%nvalue((compound_concentrations[l] div 2)+1,[ul_experiment_row[lower,r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i]| i in 1..min(compound_concentrations[l],num_rows_line)]) +%))) endif; + +%constraint if concentrations_on_different_rows then (forall(l in 1..compounds) (forall(r in 0..(compound_replicates[l]-1)) ( +%nvalue(ul_concentrations[r*compounds + l,lower]+1,[ul_experiment_row[lower,r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i]| i in %1..min(compound_concentrations[l],num_rows_line)]) +%))) endif; + +%% MAYBE XXXX +%constraint if concentrations_on_different_rows then (forall(l in 1..compounds) (forall(v in Vertical, r in 0..(compound_replicates[l]-1)) ( +%alldifferent_except_0([ul_experiment_row[v,r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i]| i in 1..min(compound_concentrations[l],num_rows_line)]) +%))) endif; + + +%constraint count(e in 1..experiments)(ul_experiment_row[upper,e]>0) = sum(ul_half_plates[..,upper]); +%constraint count(e in 1..experiments)(ul_experiment_row[lower,e]>(num_rows_line div 2)) = sum(ul_half_plates[..,lower]); +%%%%% END MAYBE + + +%% IDEA: +% Implied constraint: when there are more concentrations than rows, help the search to not fill in a row completely too soon by saving space in each row for the other compounds +% constraint if (concentrations_on_different_rows /\ (max_compound_concentrations >= num_rows_line)) then +% forall(i in Plates, j in Rows, r in 1..(compounds*replicates-1)) +% ((count(experiment_plate_row[i,1..max_compound_concentrations*r],j) + sum([ experiment_plate_row[i,max_compound_concentrations*(k-1)+1]>0 /\ (compound_concentrations[((k-1) mod compounds)+1]>=num_rows_line) | k in (r+1) .. compounds*replicates])) <= experiments_in_plate_row[i,j]) endif; + +%int: temp = (compounds*replicates-1); + +% NOT WORKING!! +% constraint if concentrations_on_different_rows then +% forall(i in Plates, j in Rows, r in temp..(compounds*replicates-1))( +% let{ +% array [Plates,Rows] of var 0..1: temp_array; +% %compound_concentrations[(r mod compounds)+1]; +% } in sum(temp_array[i,..])==1 /\ count(experiment_plate_row[i,1..max_compound_concentrations*r],j) + temp_array[i,j] == experiments_in_plate_row[i,j] +% ) endif; + +%% array [Plates,Rows] of var 0..num_cols_line: experiments_in_plate_row; +%% constraint forall(i in Plates, j in Rows)(count(experiment_plate_row[i,..],j,experiments_in_plate_row[i,j])); + + + %int: temp = (compounds*replicates-1); +%%% array [Plates,Rows] of var 0..1: temp_array2; +%%% constraint forall(i in Plates, j in Rows, r in temp..(compounds*replicates-1))( +%%% count(experiment_plate_row[i,1..max_compound_concentrations*r],j) + temp_array2[i,j] == experiments_in_plate_row[i,j] +%%% ); + +%%%constraint forall(i in Plates)(if experiment_plate_row[i,max_compound_concentrations*compounds*replicates] == 0 then sum(temp_array2[i,..])==0 else sum(temp_array2[i,..])==2 %%%endif); + +array [Plates,1..compounds*replicates,Rows] of var 0..1: compounds_in_plate_row; +array [Plates,1..compounds*replicates,Columns] of var 0..1: compounds_in_plate_column; + +%% Constraints for compounds_in_plate_row +constraint forall(j in Rows, e in 1..experiments)( + (experiment_row[e]==j) -> (compounds_in_plate_row[experiment_plate[e], (((e-1) div max_compound_concentrations)+1),j] == 1) +); + +constraint forall(i in Plates, j in Rows, e in 1..compounds*replicates)( + (experiment_plate[(e-1)*max_compound_concentrations+1] != i) -> (compounds_in_plate_row[i,e,j] == 0) +); + +constraint if concentrations_on_different_rows then forall(i in Plates, e in 1..compounds*replicates)( + if experiment_plate[(e-1)*max_compound_concentrations+1]==i then sum(compounds_in_plate_row[i,e,..]) == min(compound_concentrations[((e-1) mod compounds)+1],num_rows_line) + else sum(compounds_in_plate_row[i,e,..]) == 0 endif +) endif; + + +constraint if max_compound_concentrations<=num_rows_line /\ concentrations_on_different_rows then forall(i in Plates, j in Rows)( + sum(compounds_in_plate_row[i,..,j]) == experiments_in_plate_row[i,j] +) +else forall(i in Plates, j in Rows)( + sum(compounds_in_plate_row[i,..,j]) <= experiments_in_plate_row[i,j] +) +endif; + +constraint if max_compound_concentrations<=num_rows_line /\ concentrations_on_different_rows then + sum(compounds_in_plate_row[..,..,..]) == sum([compound_concentrations[i]*compound_replicates[i] | i in 1..compounds]) else + sum(compounds_in_plate_row[..,..,..]) <= sum([compound_concentrations[i]*compound_replicates[i] | i in 1..compounds]) +endif; + + +constraint if concentrations_on_different_rows /\ max_compound_concentrations<=num_rows_line /\ even_rows then (forall(l in 1..compounds,r in 0..(replicates-1)) ( + ul_concentrations[r*compounds + l,upper] == sum(compounds_in_plate_row[..,r*compounds + l, 1..(num_rows_line div 2)]) +)) endif; + +constraint if concentrations_on_different_rows /\ max_compound_concentrations<=num_rows_line /\ even_rows then (forall(l in 1..compounds,r in 0..(replicates-1)) ( + ul_concentrations[r*compounds + l,lower] == sum(compounds_in_plate_row[..,r*compounds + l, (num_rows_line div 2)+1..num_rows_line]) +)) endif; + +constraint if concentrations_on_different_rows /\ even_rows then + forall(i in Plates)(sum(experiments_in_plate_row[i,1..(num_rows_line div 2)]) == ul_half_plates[i, upper]) endif; + +constraint if concentrations_on_different_rows /\ even_rows then + forall(i in Plates)(sum(experiments_in_plate_row[i,(num_rows_line div 2)+1..num_rows_line]) == ul_half_plates[i, lower]) endif; + + +%% Constraints for compounds_in_plate_column +constraint forall(k in Columns, e in 1..experiments)( + (experiment_column[e]==k) -> (compounds_in_plate_column[experiment_plate[e], (((e-1) div max_compound_concentrations)+1),k] == 1) +); + +constraint forall(i in Plates, k in Columns, e in 1..compounds*replicates)( + (experiment_plate[(e-1)*max_compound_concentrations+1] != i) -> (compounds_in_plate_column[i,e,k] == 0) +); + +constraint if concentrations_on_different_columns then forall(i in Plates, e in 1..compounds*replicates)( + if experiment_plate[(e-1)*max_compound_concentrations+1]==i then sum(compounds_in_plate_column[i,e,..]) == min(compound_concentrations[((e-1) mod compounds)+1],num_cols_line) + else sum(compounds_in_plate_column[i,e,..]) == 0 endif +) endif; + + +constraint if max_compound_concentrations<=num_cols_line /\ concentrations_on_different_columns then forall(i in Plates, k in Columns)( + sum(compounds_in_plate_column[i,..,k]) == experiments_in_plate_column[i,k] +) +else forall(i in Plates, k in Columns)( + sum(compounds_in_plate_column[i,..,k]) <= experiments_in_plate_column[i,k] +) +endif; + +constraint if max_compound_concentrations<=num_cols_line /\ concentrations_on_different_columns then + sum(compounds_in_plate_column[..,..,..]) == sum([compound_concentrations[i]*compound_replicates[i] | i in 1..compounds]) else + sum(compounds_in_plate_column[..,..,..]) <= sum([compound_concentrations[i]*compound_replicates[i] | i in 1..compounds]) +endif; + +constraint if concentrations_on_different_columns /\ max_compound_concentrations<=num_cols_line /\ even_columns then (forall(l in 1..compounds,r in 0..(replicates-1)) ( + lr_concentrations[r*compounds + l,left] == sum(compounds_in_plate_column[..,r*compounds + l, 1..(num_cols_line div 2)]) +)) endif; + +constraint if concentrations_on_different_columns /\ max_compound_concentrations<=num_cols_line /\ even_columns then (forall(l in 1..compounds,r in 0..(replicates-1)) ( + lr_concentrations[r*compounds + l,right] == sum(compounds_in_plate_column[..,r*compounds + l, (num_cols_line div 2)+1..num_cols_line]) +)) endif; + +%% HERE +constraint if concentrations_on_different_columns /\ even_columns then + forall(i in Plates)(sum(experiments_in_plate_column[i,1..(num_cols_line div 2)]) == lr_half_plates[i, left]) endif; + +constraint if concentrations_on_different_columns /\ even_columns then + forall(i in Plates)(sum(experiments_in_plate_column[i,(num_cols_line div 2)+1..num_cols_line]) == lr_half_plates[i, right]) endif; + + +% array [Rows] of var 0..2: temp_array2; +% constraint forall(i in Plates, j in Rows, r in temp-1..(compounds*replicates-2))( +% count(experiment_plate_row[i,1..max_compound_concentrations*r],j) + temp_array2[j]<= experiments_in_plate_row[i,j] +% ); + +% constraint sum(temp_array2) ==20; + + +% array [Rows] of var 0..3: temp_array3; +% constraint forall(i in Plates, j in Rows, r in temp-2..(compounds*replicates-3))( +% count(experiment_plate_row[i,1..max_compound_concentrations*r],j) + temp_array3[j]<= experiments_in_plate_row[i,j] +% ); + +% constraint sum(temp_array3) ==30; + + +% array [Rows] of var 0..4: temp_array4; +% constraint forall(i in Plates, j in Rows, r in temp-3..(compounds*replicates-4))( +% count(experiment_plate_row[i,1..max_compound_concentrations*r],j) + temp_array4[j]<= experiments_in_plate_row[i,j] +% ); + +% constraint sum(temp_array4) ==40; + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Split concentrations between upper and lower halves of the plate + +array [1..compounds*replicates,Vertical] of var 0..max_compound_concentrations: ul_concentrations; + +%% For each replicate of each compound, count how many concentrations are in the upper/lower parts of the plate +%% Definition of ul_concentrations (upper) + +constraint if concentrations_on_different_rows then (forall(l in 1..compounds) (forall(r in 0..(replicates-1)) ( + among(ul_concentrations[r*compounds + l,upper],[experiment_row[r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i]| i in 1..compound_concentrations[l]], 1..(num_rows_line div 2))::domain +))) endif; + +%% Definition of ul_concentrations (lower) +constraint if concentrations_on_different_rows then (forall(l in 1..compounds) (forall(r in 0..(replicates-1)) ( + among(ul_concentrations[r*compounds + l,lower],[experiment_row[r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i]| i in 1..compound_concentrations[l]], (num_rows_line div 2)+1..num_rows_line)::domain +))) endif; + +%% Constraints to divide the concentrations into top and bottom +constraint if concentrations_on_different_rows then (forall(l in 1..compounds, v in Vertical) (forall(r in 0..(compound_replicates[l]-1)) ( + ul_concentrations[r*compounds + l,v] <= ((compound_concentrations[l]+1) div 2) +))) endif; + +constraint if concentrations_on_different_rows then (forall(l in 1..compounds, v in Vertical) (forall(r in 0..(compound_replicates[l]-1)) ( + ul_concentrations[r*compounds + l,v] >= (compound_concentrations[l] div 2) +))) endif; + +constraint forall(l in 1..compounds, v in Vertical) (forall(r in compound_replicates[l]..replicates-1) ( + ul_concentrations[r*compounds + l,v] == 0 +)); + +%% Implied constraints: +constraint if concentrations_on_different_rows then (sum(ul_concentrations)=sum([compound_concentrations[i]*compound_replicates[i] | i in 1..compounds])) endif; + +constraint if concentrations_on_different_rows then (forall(l in 1..compounds) (forall(r in 0..(compound_replicates[l]-1)) ( + ul_concentrations[r*compounds + l,upper] + ul_concentrations[r*compounds + l,lower] == compound_concentrations[l] +))) endif; + +constraint if concentrations_on_different_rows then forall(v in Vertical) (sum(ul_concentrations[..,v]) == sum(ul_half_plates[..,v]) ) endif; + +constraint if concentrations_on_different_rows then forall(i in Plates, v in Vertical) (sum(c in 0..compounds*replicates-1)(ul_concentrations[c+1,v]*(experiment_plate[c*max_compound_concentrations+1]==i)) == ul_half_plates[i,v] ) endif; + + +%% Defining ul_concentrations in terms of plates +constraint if concentrations_on_different_rows then (forall(l in 1..compounds) (forall(r in 0..(replicates-1)) ( + among(ul_concentrations[r*compounds + l,upper], plates[..,1..(num_rows_line div 2),..], {r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i| i in 1..compound_concentrations[l]})::domain +))) endif; + +constraint if concentrations_on_different_rows then (forall(l in 1..compounds) (forall(r in 0..(replicates-1)) ( + among(ul_concentrations[r*compounds + l,lower], plates[..,(num_rows_line div 2)+1..num_rows_line,..], {r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i| i in 1..compound_concentrations[l]})::domain +))) endif; + +%% Defining lr_concentrations in terms of plates +constraint if concentrations_on_different_columns then (forall(l in 1..compounds) (forall(r in 0..(replicates-1)) ( + among(lr_concentrations[r*compounds + l,left], plates[..,..,1..(num_cols_line div 2)], {r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i| i in 1..compound_concentrations[l]})::domain +))) endif; + +constraint if concentrations_on_different_columns then (forall(l in 1..compounds) (forall(r in 0..(replicates-1)) ( + among(lr_concentrations[r*compounds + l,right], plates[..,..,(num_cols_line div 2)+1..num_cols_line], {r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i| i in 1..compound_concentrations[l]})::domain +))) endif; + + +%% Implied constraint +constraint if concentrations_on_different_rows /\ even_rows then (forall(v in Vertical)( +sum(ul_half_plates[..,v]) == sum(ul_concentrations[..,v]) +))endif; + + +%% Definition of experiments_in_plate_row +constraint forall(i in Plates, j in Rows)(count(experiment_plate_row[i,..],j,experiments_in_plate_row[i,j])); + +constraint forall(i in Plates, j in Rows)(among(experiments_in_plate_row[i,j],plates[i,j,..],1..experiments)::domain); + +constraint forall(i in Plates, k in Columns)(among(experiments_in_plate_column[i,k],plates[i,..,k],1..experiments)::domain); + +%constraint forall(i in Plates, k in Columns)(experiments_in_plate_column[i,k] == sum(j in Rows)(0= num_rows_line)) then forall(i in Plates, j in Rows, r in 1..compounds*replicates)( +count(experiment_plate_row[i,1..experiments-max_compound_concentrations*r],j) + count(experiment_plate_row[i,experiments-max_compound_concentrations*r+1..experiments],j) == experiments_in_plate_row[i,j] +) endif; + +constraint if (concentrations_on_different_rows /\ (max_compound_concentrations >= num_rows_line)) then forall(i in Plates, j in Rows, r in 1..compounds*replicates-1)( +if(compound_concentrations[((r-1) mod compounds)+1]>=num_rows_line) then +(count(experiment_plate_row[i,1..max_compound_concentrations*r],j) + count(experiment_plate_row[i,max_compound_concentrations*(r+1)+1..experiments],j) + (experiment_plate_row[i,max_compound_concentrations*r+1]>0)) <= experiments_in_plate_row[i,j] +endif +) endif; + +% Implied constraint: when there are more concentrations than rows, help the search to not fill in a row completely too soon by saving space in each row for the other compounds +constraint if (concentrations_on_different_rows /\ (max_compound_concentrations >= num_rows_line)) then forall(i in Plates, j in Rows, r in 1..(compounds*replicates-1))((count(experiment_plate_row[i,1..max_compound_concentrations*r],j) + sum([ experiment_plate_row[i,max_compound_concentrations*(k-1)+1]>0 /\ (compound_concentrations[((k-1) mod compounds)+1]>=num_rows_line) | k in (r+1) .. compounds*replicates])) <= experiments_in_plate_row[i,j]) endif; + + + +%%%%%% +%% Spreads concentrations over different columns. It's been replaced by the more restrictive alldifferent version below. +%constraint if concentrations_on_different_columns then (forall(l in 1..compounds) (forall(r in 0..(compound_replicates[l]-1)) (nvalue(min(compound_concentrations[l],num_cols_line-2*size_empty_edge), [experiment_column[r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i]| i in 1..compound_concentrations[l]])))) endif; + +constraint if concentrations_on_different_columns then (forall(l in 1..compounds) (forall(r in 0..(compound_replicates[l]-1)) ( +alldifferent( [experiment_column[r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i]| i in 1..min(compound_concentrations[l],num_cols_line)]):: domain +))) endif; + +% Split concentrations between left-most and right-most halves of the plate + +array [1..compounds*replicates,Horizontal] of var 0..max_compound_concentrations: lr_concentrations; + +%% For each replicate of each compound, count how many concentrations there are in the left/right halves of the plate +%% Definition of lr_concentrations (left) +constraint if concentrations_on_different_columns /\ even_columns then (forall(l in 1..compounds) (forall(r in 0..(replicates-1)) ( + among(lr_concentrations[r*compounds + l,left],[experiment_column[r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i]| i in 1..compound_concentrations[l]], 1..(num_cols_line div 2))::domain +))) endif; + + +%% Definition of lr_concentrations (right) +constraint if concentrations_on_different_columns /\ even_columns then (forall(l in 1..compounds) (forall(r in 0..(replicates-1)) ( + among(lr_concentrations[r*compounds + l,right],[experiment_column[r*compounds*max_compound_concentrations+(l-1)*max_compound_concentrations+i]| i in 1..compound_concentrations[l]], (num_cols_line div 2)+1..num_cols_line)::domain +))) endif; + +%% Balancing left and right concentrations +constraint if concentrations_on_different_columns /\ even_columns then (forall(l in 1..compounds, h in Horizontal) (forall(r in 0..(compound_replicates[l]-1)) ( + lr_concentrations[r*compounds + l,h] <= ((compound_concentrations[l] + 1) div 2) +))) endif; + +constraint if concentrations_on_different_columns /\ even_columns then (forall(l in 1..compounds, h in Horizontal) (forall(r in 0..(compound_replicates[l]-1)) ( + lr_concentrations[r*compounds + l,h] >= (compound_concentrations[l] div 2) +))) endif; + +constraint forall(l in 1..compounds, h in Horizontal) (forall(r in compound_replicates[l]..replicates-1) ( + lr_concentrations[r*compounds + l,h] == 0 +)); + +%% Implied constraints: +constraint if concentrations_on_different_columns /\ even_columns then (sum(lr_concentrations) == sum([compound_concentrations[i]*compound_replicates[i] | i in 1..compounds])) endif; + +constraint if concentrations_on_different_columns /\ even_columns then (forall(l in 1..compounds) (forall(r in 0..(compound_replicates[l]-1)) ( + lr_concentrations[r*compounds + l,left] + lr_concentrations[r*compounds + l,right] == compound_concentrations[l] +))) endif; + +constraint if concentrations_on_different_columns /\ even_columns then (forall(h in Horizontal)( +sum(lr_half_plates[..,h]) == sum(lr_concentrations[..,h]) +))endif; + +constraint if concentrations_on_different_columns then forall(i in Plates, h in Horizontal) (sum(c in 0..compounds*replicates-1)(lr_concentrations[c+1,h]*(experiment_plate[c*max_compound_concentrations+1]==i)) == lr_half_plates[i,h] ) endif; + +%% Implied constraints +%% These constraints are too weak when the number of plates increases +%% Right now it looks like it's mostly taking time + +constraint global_cardinality_low_up(experiment_row, [i | i in Rows], [0| i in Rows], [numplates*num_cols_line| i in Rows]); +constraint global_cardinality_low_up(experiment_column, [i | i in Columns], [0| i in Columns], [numplates*num_rows_line| i in Columns]); + +array [Plates,Rows] of var 0..num_cols_line: controls_in_plate_row; +array [Plates,Columns] of var 0..num_rows_line: controls_in_plate_column; +array [Plates,Rows] of var 0..num_cols_line: controls_only_in_plate_row; +array [Plates,Columns] of var 0..num_rows_line: controls_only_in_plate_column; +array [Plates,Rows] of var 0..num_cols_line: emptywells_in_plate_row; +array [Plates,Columns] of var 0..num_rows_line: emptywells_in_plate_column; +array [Rows] of var 0..numplates*num_cols_line: controls_in_row; +array [Columns] of var 0..numplates*num_rows_line: controls_in_column; + +constraint (sum(controls_in_plate_row) == total_controls + emptywells); +constraint (sum(controls_in_plate_column) == total_controls + emptywells); + +constraint (sum(controls_only_in_plate_row) == total_controls); +constraint (sum(controls_only_in_plate_column) == total_controls); + +constraint (sum(controls_in_row) == total_controls + emptywells); +constraint (sum(controls_in_column) == total_controls + emptywells); + +constraint (sum(emptywells_in_plate_row) == emptywells); +constraint (sum(emptywells_in_plate_column) == emptywells); + +%% Redundant +%constraint (sum(controls_in_plate_row) == sum(controls_only_in_plate_row) + emptywells); +%constraint (sum(controls_in_plate_column) == sum(controls_only_in_plate_column) + emptywells); + +constraint forall(i in Plates, j in Rows)(count(plates[i,j,..],0,emptywells_in_plate_row[i,j])); +constraint forall(i in Plates, k in Columns)(count(plates[i,..,k],0,emptywells_in_plate_column[i,k])); + +%Among: Requires exactly n variables in x to take one of the values in v . +constraint forall(i in Plates, j in Rows) (among(controls_in_plate_row[i,j], plates[i,j,..], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); +constraint forall(i in Plates, j in Rows) (among(controls_only_in_plate_row[i,j], plates[i,j,..], (experiments+1..experiments+num_controls*max_control_concentrations))::domain); + + +constraint forall(i in Plates, k in Columns) (among(controls_in_plate_column[i,k], plates[i,..,k], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); +constraint forall(i in Plates, k in Columns) (among(controls_only_in_plate_column[i,k], plates[i,..,k], (experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(i in Plates, j in Rows) (among(emptywells_in_plate_row[i,j], plates[i,j,..], {0})::domain); +constraint forall(i in Plates, k in Columns) (among(emptywells_in_plate_column[i,k], plates[i,..,k], {0})::domain); + +constraint forall(j in Rows) (among(controls_in_row[j], plates[..,j,..], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(k in Columns) (among(controls_in_column[k], plates[..,..,k], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(j in Rows) (controls_in_row[j] == sum(controls_in_plate_row[..,j])); + +constraint forall(j in Rows) (controls_in_row[j] == sum(emptywells_in_plate_row[..,j]) + sum(controls_only_in_plate_row[..,j])); + +constraint forall(k in Columns) (controls_in_column[k] == sum(controls_in_plate_column[..,k])); + +constraint forall(k in Columns) (controls_in_column[k] == sum(emptywells_in_plate_column[..,k]) + sum(controls_only_in_plate_column[..,k])); + +constraint forall(i in Plates) ( sum(controls_ul_plates[i,upper,..]) + count(plates[i,1..(num_rows_line div 2),..],0) == sum(controls_in_plate_row[i,1..(num_rows_line div 2)])); + +constraint forall(i in Plates) ( sum(controls_ul_plates[i,lower,..]) + count(plates[i,(num_rows_line div 2)+1..num_rows_line,..],0) == sum(controls_in_plate_row[i,(num_rows_line div 2)+1..num_rows_line])); + +constraint forall(i in Plates) ( sum(controls_lr_plates[i,left,..]) + count(plates[i,..,1..(num_cols_line div 2)],0) == sum(controls_in_plate_column[i,1..(num_cols_line div 2)])); +constraint forall(i in Plates) ( sum(controls_lr_plates[i,right,..]) + count(plates[i,..,(num_cols_line div 2)+1..num_cols_line],0) == sum(controls_in_plate_column[i,(num_cols_line div 2)+1..num_cols_line])); + +constraint forall(i in Plates) ( sum_controls_ul_plates[i,upper] == sum(j in 1..(num_rows_line div 2))(controls_only_in_plate_row[i,j])); +constraint forall(i in Plates) ( sum_controls_ul_plates[i,lower] == sum(j in (num_rows_line div 2)+1..num_rows_line)(controls_only_in_plate_row[i,j])); + +constraint forall(i in Plates) ( sum_controls_lr_plates[i,left] == sum(j in 1..(num_cols_line div 2))(controls_only_in_plate_column[i,j])); +constraint forall(i in Plates) ( sum_controls_lr_plates[i,right] == sum(j in (num_cols_line div 2)+1..num_cols_line)(controls_only_in_plate_column[i,j])); + + +constraint sum(sum_controls_ul_plates) == total_controls; +constraint sum(sum_controls_lr_plates) == total_controls; + +constraint forall(v in Vertical) (sum(sum_controls_ul_plates[..,v]) >= (total_controls div 2)); +constraint forall(v in Vertical) (sum(sum_controls_ul_plates[..,v]) <= ((total_controls + 1) div 2)); + +constraint forall(h in Horizontal) (sum(sum_controls_lr_plates[..,h]) >= (total_controls div 2)); +constraint forall(h in Horizontal) (sum(sum_controls_lr_plates[..,h]) <= ((total_controls + 1) div 2)); + +%% Implied constraint: how many experiments can be in a row/column? It's too weak and needs to be done per plate. +constraint global_cardinality(experiment_row, [j | j in Rows], [(numplates*num_cols_line)-controls_in_row[j]| j in Rows]); +constraint global_cardinality(experiment_column, [k | k in Columns], [(numplates*num_rows_line)-controls_in_column[k]| k in Columns]); + + + + +%% Implied constraint: how many experiments can there be in a row? Per plate +% I think this still works... +%% TODO: Redundant +constraint forall(i in Plates, j in Rows) (sum([ experiment_row[l] == j /\ experiment_plate[l] == i| l in 1..experiments])== num_cols_line-among(plates[i,j,..], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(i in Plates, j in Rows) (sum([ experiment_row[l] == j /\ experiment_plate[l] == i| l in 1..experiments])== num_cols_line-controls_in_plate_row[i,j]); + +constraint forall(i in Plates, j in Rows) (sum([ experiment_row[l] == j /\ experiment_plate[l] == i| l in 1..experiments])== num_cols_line-emptywells_in_plate_row[i,j]-controls_only_in_plate_row[i,j]); + +constraint forall(i in Plates, k in Columns) (sum([ experiment_column[l] == k /\ experiment_plate[l] == i| l in 1..experiments])== num_rows_line-among(plates[i,..,k], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(i in Plates, k in Columns) (sum([ experiment_column[l] == k /\ experiment_plate[l] == i| l in 1..experiments])== num_rows_line-controls_in_plate_column[i,k]); + +constraint forall(i in Plates, k in Columns) (sum([ experiment_column[l] == k /\ experiment_plate[l] == i| l in 1..experiments])== num_rows_line-emptywells_in_plate_column[i,k]-controls_only_in_plate_column[i,k]); + + +constraint forall(i in Plates, j in Rows) (sum([ experiment_row[l] == j /\ experiment_plate[l] == i| l in 1..experiments]) == experiments_in_plate_row[i,j]); + +constraint forall(i in Plates, k in Columns) (sum([ experiment_column[l] == k /\ experiment_plate[l] == i| l in 1..experiments]) == experiments_in_plate_column[i,k]); + +%experiments_in_plate_row +constraint forall(i in Plates, j in Rows) (experiments_in_plate_row[i,j] == num_cols_line-among(plates[i,j,..], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(i in Plates, j in Rows) (experiments_in_plate_row[i,j] == num_cols_line-controls_in_plate_row[i,j]); + +constraint forall(i in Plates, j in Rows) (experiments_in_plate_row[i,j] == num_cols_line-emptywells_in_plate_row[i,j]-controls_only_in_plate_row[i,j]); + +constraint forall(i in Plates, k in Columns) (experiments_in_plate_column[i,k] == num_rows_line-among(plates[i,..,k], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(i in Plates, k in Columns) (experiments_in_plate_column[i,k] == num_rows_line-controls_in_plate_column[i,k]); + +constraint forall(i in Plates, k in Columns) (experiments_in_plate_column[i,k] == num_rows_line-emptywells_in_plate_column[i,k]-controls_only_in_plate_column[i,k]); + + +%% Below suggested by Gustav. It improves flattening but worsens the solving time +% constraint forall(i in Plates, j in Rows) ( +% let{ +% var int: exprInRow; +% constraint exprInRow = among(plates[i,j,..], 1..experiments)::domain; +% %sum([experiment_row[l] == j /\ experiment_plate[l] == i| l in 1..experiments]); +% } in exprInRow <= num_cols_line-among([plates[i,j,k] | k in Columns], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))); + + +%% Replaced by constraint below (suggested by Gustav) +%constraint forall(i in Plates, k in Columns) (sum([ experiment_column[l] == k /\ experiment_plate[l] == i| l in 1..experiments])<= num_rows_line-among(plates[i,..,k], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))); + +%% Gustav: +constraint forall(i in Plates, k in Columns) ( +% sum([ experiment_column[l] == k /\ experiment_plate[l] == i| l in 1..experiments]) +(among(plates[i,..,k], 1..experiments)::domain) <= num_rows_line-among([plates[i,j,k] | j in Rows], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Different compound replicas must appear on different plates (if possible) +int: min_plates = min(numplates,replicates); + +constraint if replicates_on_different_plates then forall(l in 1..compounds) (nvalue(min(numplates,compound_replicates[l]), [experiment_plate[(l-1)*max_compound_concentrations + i*compounds*max_compound_concentrations + 1] | i in 0..(compound_replicates[l]-1)])) endif; + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Different compound replicas must appear on the same plate +constraint if replicates_on_same_plate then forall(l in 1..compounds) (all_equal([experiment_plate[(l-1)*max_compound_concentrations + i*compounds*max_compound_concentrations + 1] | i in 0..(compound_replicates[l]-1)])) endif; + + +constraint if concentrations_on_different_rows then forall(l in 1..compounds) (forall(conc in 0..compound_concentrations[l]-1)(alldifferent([experiment_row[(l-1)*max_compound_concentrations + i*compounds*max_compound_concentrations + conc + 1] | i in 0..(min(compound_replicates[l],num_rows_line)-1)])::domain)) endif; + +constraint if concentrations_on_different_columns then forall(l in 1..compounds) (forall(conc in 0..compound_concentrations[l]-1)(alldifferent([experiment_column[(l-1)*max_compound_concentrations + i*compounds*max_compound_concentrations + conc + 1] | i in 0..(min(compound_replicates[l],num_cols_line)-1)])::domain)) endif; + + +%% Balancing compounds between plates: distribute the compounds equitatively among all the plates + + +%% TODO: could this lower bound be stronger? +%% FIX ME! +%int: at_least = max_compound_concentrations*floor((compounds*replicates)/numplates) + combination_concentrations*floor((combinations*replicates)/numplates); + +%int: at_least = min(compound_concentrations)*floor(sum(compound_replicates)/numplates) + combination_concentrations*floor((combinations*replicates)/numplates); + +%% FIX ME! +%% FIX ME!! +%int: at_most = min(max_compound_concentrations*ceil(sum(compound_replicates)/numplates) + combination_concentrations*ceil((combinations*replicates)/numplates), + % inner_plate_size -(sum([floor(control_replicates[i]*control_concentrations[i]/numplates) | i in 1..num_controls]))); + + + +%% This implied constraint is very important for propagation! + +%constraint global_cardinality_low_up(experiment_plate,[ i | i in Plates],[ at_least | i in Plates],[ at_most | i in Plates]); + +% I could make this stronger, but I would need to assume that experiments can indeed be distrubuted in a balanced way. +% I would like to do something that allows example-02 and example-03 to work +% This constraint assumes that there is no plate with only controls, but that seems fair enough + + +int: min_compounds_plate = max(floor((sum([compound_concentrations[i]*compound_replicates[i] | i in 1..compounds]))/inner_plate_size), min(compound_concentrations++[0]) ); + +constraint global_cardinality_low_up(experiment_plate,[ i | i in Plates],[ min_compounds_plate | i in Plates],[ inner_plate_size | i in Plates]); + +%% TODO: Think more about this constraint and what happens when there are replicates, etc. +%% This option is possibly going to make some plans unfeasible (when there are different numbers of replicates, etc). +constraint if sorted_compounds == true /\ max(compound_replicates)<=1 then sorted_except_0(experiment_plate,numplates,at_least_compounds,at_most_compounds) endif; + +%% WARNING! This option assumes that compounds can be placed on plates in order! This works, for example, when all compounds have the same +% number of replicates and concentrations and one would like to generate several layouts with the compounds on the exact same plates. +% Otherwise, experiment_plate needs to be given as input +constraint if sorted_compounds == true /\ replicates_on_same_plate == true then sorted_except_0(experiment_plate[1..compounds*max(compound_concentrations)],numplates,at_least_compounds div numplates,at_most_compounds div numplates) endif; + + + +% Balancing compounds between plates (and knowing that there must be some sort of balance!) +%% FIX ME! +%% I want to remove these constraints but I need to make other stronger first +int: at_least_compounds = min(compound_concentrations++[infinity])*(sum(compound_replicates) div numplates); +int: at_most_compounds = max_compound_concentrations*ceil(sum(compound_replicates)/numplates); + +%% TODO: False only if there are no compounds +constraint if at_least_compounds <= at_most_compounds then global_cardinality_low_up([experiment_plate[i] | i in 1..experiments] ,[ i | i in Plates],[ at_least_compounds | i in Plates],[ at_most_compounds | i in Plates]) endif; + + + + +%% Implied constraints +%% These constraints are too weak when the number of plates increases +constraint if concentrations_on_different_rows then global_cardinality_low_up(experiment_row, [i | i in Rows], [numplates*(at_least_compounds div num_rows_line)| i in Rows], [numplates*num_cols_line| i in Rows]) endif; + +constraint if concentrations_on_different_columns then global_cardinality_low_up(experiment_column, [i | i in Columns], [numplates*floor(at_least_compounds/num_cols_line)| i in Columns], [numplates*num_rows_line| i in Columns]) endif; + +constraint if concentrations_on_different_rows then forall(j in Rows)(sum(experiments_in_plate_row[..,j])>=numplates*(at_least_compounds div num_rows_line)) endif; +constraint if concentrations_on_different_columns then forall(k in Columns)(sum(experiments_in_plate_column[..,k])>=numplates*(at_least_compounds div num_cols_line)) endif; + +constraint if concentrations_on_different_rows then forall(i in Plates, j in Rows)(experiments_in_plate_row[i,j]>=(at_least_compounds div num_rows_line)) endif; +constraint if concentrations_on_different_columns then forall(i in Plates, k in Columns)(experiments_in_plate_column[i,k]>=(at_least_compounds div num_cols_line)) endif; + +constraint if concentrations_on_different_rows then forall(i in Plates, j in Rows)(experiments_in_plate_row[i,j]<=ceil(at_most_compounds / num_rows_line)) endif; +constraint if concentrations_on_different_columns then forall(i in Plates, k in Columns)(experiments_in_plate_column[i,k]<=ceil(at_most_compounds / num_cols_line)) endif; + +array [Plates,Vertical] of var 0..num_cols_line*(num_rows_line div 2): ul_half_plates; +array [Plates,Horizontal] of var 0..num_rows_line*(num_cols_line div 2): lr_half_plates; + +%%% Counting the number of experiments in the upper and lower half-plates +constraint forall(i in Plates) (among(ul_half_plates[i,upper],plates[i,1..num_rows_line div 2,..], 1..experiments)::domain); +constraint forall(i in Plates) (among(ul_half_plates[i,lower],plates[i,(num_rows_line div 2)+1..num_rows_line,..], 1..experiments)::domain); + +constraint forall(i in Plates) (sum(e in 1..experiments)(experiment_plate[e] == i /\ experiment_row[e] <= (num_rows_line div 2)) == ul_half_plates[i,upper]); +constraint forall(i in Plates) (sum(e in 1..experiments)(experiment_plate[e] == i /\ experiment_row[e] > (num_rows_line div 2)) == ul_half_plates[i,lower]); + +constraint (sum(e in 1..experiments)(experiment_row[e] in 1..(num_rows_line div 2)) == sum(i in Plates)(ul_half_plates[i,upper])); +constraint (sum(e in 1..experiments)(experiment_row[e] in (num_rows_line div 2)+1..num_rows_line) == sum(i in Plates)(ul_half_plates[i,lower])); + +constraint if concentrations_on_different_rows /\ even_rows then forall(i in Plates)(sum(experiments_in_plate_row[i,1..(num_rows_line div 2)]) == ul_half_plates[i,upper]) endif; +constraint if concentrations_on_different_rows /\ even_rows then forall(i in Plates)(sum(experiments_in_plate_row[i,(num_rows_line div 2)+1..num_rows_line]) == ul_half_plates[i,lower]) endif; + +%% The following constraints should be redundant when using the two constraints above. +constraint (among(sum(i in Plates)(ul_half_plates[i,upper]),experiment_row, 1..(num_rows_line div 2))::domain); +constraint (among(sum(i in Plates)(ul_half_plates[i,lower]),experiment_row, (num_rows_line div 2)+1..num_rows_line)::domain); + +constraint (among(sum(i in Plates)(lr_half_plates[i,left]),experiment_column, 1..(num_cols_line div 2))::domain); +constraint (among(sum(i in Plates)(lr_half_plates[i,right]),experiment_column, (num_cols_line div 2)+1..num_cols_line)::domain); + + +constraint if even_rows then (forall(i in Plates)(abs(ul_half_plates[i,upper]-ul_half_plates[i,lower])<=1)) endif; + +%%% Counting the number of experiments in the left and right half-plates +constraint forall(i in Plates) (among(lr_half_plates[i,left],plates[i,..,1..num_cols_line div 2], 1..experiments)::domain); +constraint forall(i in Plates) (among(lr_half_plates[i,right],plates[i,..,(num_cols_line div 2)+1..num_cols_line], 1..experiments)::domain); + +constraint forall(i in Plates) (sum(e in 1..experiments)(experiment_plate[e] == i /\ experiment_column[e] <= (num_cols_line div 2)) == lr_half_plates[i,left]); +constraint forall(i in Plates) (sum(e in 1..experiments)(experiment_plate[e] == i /\ experiment_column[e] > (num_cols_line div 2)) == lr_half_plates[i,right]); + +constraint if even_columns then (sum(e in 1..experiments)(experiment_column[e] in 1..(num_cols_line div 2)) == sum(i in Plates)(lr_half_plates[i,left])) endif; +constraint if even_columns then (sum(e in 1..experiments)(experiment_column[e] in (num_cols_line div 2)+1..num_cols_line) == sum(i in Plates)(lr_half_plates[i,right])) endif; + +constraint if even_columns then (forall(i in Plates)(abs(lr_half_plates[i,left]-lr_half_plates[i,right])<=1)) endif; + + + + +%%% Counting the number of experiments in left-most and right-most half-plates +%constraint forall(i in Plates) (among(lr_half_plates[i,left],[plates[i,j,k] | j in Rows, k in 1..ceil(num_cols_line/2)], 1..experiments)); %%::domain +%constraint forall(i in Plates) (among(lr_half_plates[i,right],[plates[i,j,k] | j in Rows, k in floor(num_cols_line/2)+1..num_cols_line], 1..experiments)); + +%constraint forall(i in Plates) (sum(e in 1..experiments, k in 1..ceil(num_cols_line/2))(experiment_plate[e] == i /\ experiment_column[e] == k) == lr_half_plates[i,left]); +%constraint forall(i in Plates) (sum(e in 1..experiments, k in ceil(num_cols_line/2)+1..num_cols_line)(experiment_plate[e] == i /\ experiment_column[e] == k) == lr_half_plates[i,right]); + + +%array [Plates,Vertical,Horizontal] of var floor((at_least_compounds+at_least_combinations)/4)..ceil((at_most_compounds+at_most_combinations)/4): quarter_plates; + +% Balancing experiments in quarters +%constraint forall(i in Plates) (among(quarter_plates[i,upper,left],[plates[i,j,k] | j in 1..ceil(num_rows_line/2), k in 1..ceil(num_cols_line/2)], 1..experiments)); +%constraint forall(i in Plates) (among(quarter_plates[i,upper,right],[plates[i,j,k] | j in 1..ceil(num_rows_line/2), k in floor(num_cols_line/2)+1..num_cols_line], 1..experiments)); + +%constraint forall(i in Plates) (among(quarter_plates[i,lower,left],[plates[i,j,k] | j in floor(num_rows_line/2)+1..num_rows_line, k in 1..ceil(num_cols_line/2)], 1..experiments)); +%constraint forall(i in Plates) (among(quarter_plates[i,lower,right],[plates[i,j,k] | j in floor(num_rows_line/2)+1..num_rows_line, k in floor(num_cols_line/2)+1..num_cols_line], 1..experiments)); + + +%%% NEW %%% +%% Implied constraints %% +constraint sum(ul_half_plates) == sum([compound_concentrations[i]*compound_replicates[i] | i in 1..compounds]); +constraint sum(lr_half_plates) == sum([compound_concentrations[i]*compound_replicates[i] | i in 1..compounds]); + +constraint forall(v in Vertical)(sum(ul_half_plates[..,v]) <= ceil(sum([compound_concentrations[i]*compound_replicates[i] | i in 1..compounds])/2)); +constraint forall(h in Horizontal)(sum(lr_half_plates[..,h]) <= ceil(sum([compound_concentrations[i]*compound_replicates[i] | i in 1..compounds])/2)); + +constraint forall(i in Plates, v in Vertical) ( ul_half_plates[i,v] <= ((sum([ (experiment_plate[c] == i) | c in 1..experiments]) + 1) div 2)); +constraint forall(i in Plates, h in Horizontal) ( lr_half_plates[i,h] <= ((sum([ (experiment_plate[c] == i) | c in 1..experiments]) + 1) div 2)); + +constraint forall(i in Plates, v in Vertical)(((floor(sum(compound_replicates)/numplates)-1)*min(compound_concentrations++[infinity]) + max_compound_concentrations) div 2 <= ul_half_plates[i,v]); + +constraint forall(i in Plates, h in Horizontal)(((floor(sum(compound_replicates)/numplates)-1)*min(compound_concentrations++[infinity]) + max_compound_concentrations) div 2 <= lr_half_plates[i,h]); + +%constraint forall(i in Plates)(sum(v in Vertical)(ul_half_plates[i,v])<=at_most); +%constraint forall(i in Plates)(sum(h in Horizontal)(lr_half_plates[i,h])<=at_most); + +%constraint forall(i in Plates)(sum(v in Vertical)(ul_half_plates[i,v])>=at_least); +%constraint forall(i in Plates)(sum(h in Horizontal)(lr_half_plates[i,h])>=at_least); + + +%constraint forall(i in Plates)(sum(v in Vertical)(ul_half_plates[i,v]) == sum(h in Horizontal)(lr_half_plates[i,h])); + +%constraint forall(i in Plates)(sum(v in Vertical)(ul_half_plates[i,v]) == sum(v in Vertical, h in Horizontal)(quarter_plates[i,v,h])); + +%% Implied constraints +constraint forall(i in Plates)(count_eq(experiment_plate, i, sum(h in Horizontal)(lr_half_plates[i,h]))); +constraint forall(i in Plates)(count_eq(experiment_plate, i, sum(v in Vertical)(ul_half_plates[i,v]))); + +%constraint forall(i in Plates)(count_eq(experiment_plate, i, sum(v in Vertical, h in Horizontal)(quarter_plates[i,v,h]))); + + + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% Balancing controls between plates. +%% Asumes controls can be balanced this way, that is, that the designer expects about the same number of controls on each plate and not a plate full of controls. + +array[int] of float: controls_per_plate = [control_replicates[floor((i-1)/max_control_concentrations)+1]*(((i-1) mod max_control_concentrations)=floor(min_controls_per_plate[d]/2)); + +constraint forall(i in Plates, d in 1..num_controls*max_control_concentrations)(abs(controls_ul_plates[i,upper,d] - controls_ul_plates[i,lower,d]) <= 1); + +constraint if even_columns then forall(i in Plates, v in Vertical, d in 1..num_controls*max_control_concentrations)(abs(controls_quarter_plates[i,v,left,d] - controls_quarter_plates[i,v,right,d]) <= 1) endif; + +constraint if even_columns then forall(i in Plates, h in Horizontal, d in 1..num_controls*max_control_concentrations)(abs(controls_quarter_plates[i,upper,h,d] - controls_quarter_plates[i,lower,h,d]) <= 1) endif; + + +constraint if even_columns then forall(i in Plates, v in Vertical, d in 1..num_controls*max_control_concentrations)(controls_quarter_plates[i,v,left,d] + controls_quarter_plates[i,v,right,d] == controls_ul_plates[i,v,d]) endif; + +constraint if even_columns then forall(i in Plates, h in Horizontal, d in 1..num_controls*max_control_concentrations)(controls_quarter_plates[i,upper,h,d] + controls_quarter_plates[i,lower,h,d] == controls_lr_plates[i,h,d]) endif; + +constraint if even_columns then forall(i in Plates, v in Vertical, d in 1..num_controls*max_control_concentrations)(controls_quarter_plates[i,v,left,d] >= (controls_ul_plates[i,v,d] div 2)) endif; + +constraint if even_columns then forall(i in Plates, v in Vertical, d in 1..num_controls*max_control_concentrations)(controls_quarter_plates[i,v,right,d] >= (controls_ul_plates[i,v,d] div 2)) endif; + +constraint if even_columns then forall(i in Plates, v in Vertical, d in 1..num_controls*max_control_concentrations)(controls_quarter_plates[i,v,left,d] <= ((controls_ul_plates[i,v,d]+1) div 2)) endif; + +constraint if even_columns then forall(i in Plates, v in Vertical, d in 1..num_controls*max_control_concentrations)(controls_quarter_plates[i,v,right,d] <= ((controls_ul_plates[i,v,d]+1) div 2)) endif; + +constraint if even_rows then forall(i in Plates, h in Horizontal, d in 1..num_controls*max_control_concentrations)(controls_quarter_plates[i,upper,h,d] >= (controls_lr_plates[i,h,d] div 2)) endif; + +constraint if even_rows then forall(i in Plates, h in Horizontal, d in 1..num_controls*max_control_concentrations)(controls_quarter_plates[i,lower,h,d] >= (controls_lr_plates[i,h,d] div 2)) endif; + +constraint if even_rows then forall(i in Plates, h in Horizontal, d in 1..num_controls*max_control_concentrations)(controls_quarter_plates[i,upper,h,d] <= ((controls_lr_plates[i,h,d]+1) div 2)) endif; + +constraint if even_rows then forall(i in Plates, h in Horizontal, d in 1..num_controls*max_control_concentrations)(controls_quarter_plates[i,lower,h,d] <= ((controls_lr_plates[i,h,d]+1) div 2)) endif; + +constraint forall(i in Plates, d in 1..num_controls*max_control_concentrations)( sum(v in Vertical)(controls_ul_plates[i,v,d]) <= max_controls_per_plate[d]); + +constraint forall(i in 1..num_controls*max_control_concentrations)(sum([controls_ul_plates[p,v,i] | v in Vertical, p in Plates]) == control_replicates[floor((i-1)/max_control_concentrations)+1]*(((i-1) mod max_control_concentrations)= (min_controls_per_plate[d] div 2)) endif; + +constraint if even_columns then forall(i in Plates, d in 1..num_controls*max_control_concentrations)(abs(controls_lr_plates[i,left,d] - controls_lr_plates[i,right,d]) <= 1) endif; + +constraint forall(i in Plates, d in 1..num_controls*max_control_concentrations)(controls_lr_plates[i,left,d] + controls_lr_plates[i,right,d] <= max_controls_per_plate[d]); + +constraint forall(i in 1..num_controls*max_control_concentrations)(sum([controls_lr_plates[p,h,i] | h in Horizontal, p in Plates]) == control_replicates[floor((i-1)/max_control_concentrations)+1]*(((i-1) mod max_control_concentrations)constraint forall(i in Plates)(ul_half_plates[i,upper] + ul_half_plates[i,lower] + sum(d in 1..num_controls*max_control_concentrations,v in Vertical)(controls_ul_plates[i,v,d]) + count([plates[i,j,k] | j in 1+size_empty_edge..num_rows_line-size_empty_edge, k in 1+size_empty_edge..num_cols_line-size_empty_edge], 0) == inner_plate_size); + + +%% Is this helping? I don't think so! But it might... +%%->constraint (sum(r in 1..floor(num_rows_line/2))(controls_in_row[r]) == sum(i in Plates, d in 1..num_controls*max_control_concentrations)(controls_ul_plates[i,upper,d]) + count([plates[i,j,k] | i in Plates, j in 1..floor(num_rows_line/2), k in Columns], 0)); + +%%->constraint (sum(r in floor(num_rows_line/2)+1..num_rows_line)(controls_in_row[r]) == sum(i in Plates, d in 1..num_controls*max_control_concentrations)(controls_ul_plates[i,lower,d]) + count([plates[i,j,k] | i in Plates, j in floor(num_rows_line/2)+1..num_rows_line, k in Columns], 0)); + + +%%% Controls per row across plates%%% + +%constraint if temp_bool /\ interconnected_plates then forall(j in 1+size_empty_edge..num_rows_line-size_empty_edge) (controls_in_row[j] <= ceil((sum(l in 1..num_controls)(control_replicates[l]*control_concentrations[l])+emptywells)/(num_rows_line-2*size_empty_edge)) + 2*size_empty_edge*numplates) endif; + +%constraint if temp_bool /\ interconnected_plates then forall(j in 1+size_empty_edge..num_rows_line-size_empty_edge) (controls_in_row[j] >= floor((sum(l in 1..num_controls)(control_replicates[l]*control_concentrations[l])+emptywells)/(num_rows_line-2*size_empty_edge)) + 2*size_empty_edge*numplates) endif; + +constraint if interconnected_plates then forall(j in Rows) (controls_in_row[j] <= ceil((total_controls+emptywells)/num_rows_line) ) endif; + +constraint if interconnected_plates then forall(j in Rows) (controls_in_row[j] >= floor((total_controls+emptywells)/num_rows_line) ) endif; + +constraint forall(i in Plates, j in Rows) (controls_only_in_plate_row[i,j] <= ceil(ceil(total_controls/numplates)/num_rows_line) ) ; + +constraint forall(i in Plates, j in Rows) (controls_only_in_plate_row[i,j] >= ((total_controls div numplates) div num_rows_line) ) ; + + +%%% Controls per row per plate %%% +%% Relaxing this constraint a bit... +%constraint forall(i in Plates, j in 1+size_empty_edge..num_rows_line-size_empty_edge) (controls_in_plate_row[i,j] <= control_slack + ceil(ceil((sum(l in 1..num_controls)(control_replicates[l]*control_concentrations[l])+emptywells)/(num_rows_line-2*size_empty_edge))/numplates) + 2*size_empty_edge); + +%constraint forall(i in Plates, j in 1+size_empty_edge..num_rows_line-size_empty_edge) (controls_in_plate_row[i,j] >= floor(floor((sum(l in 1..num_controls)(control_replicates[l]*control_concentrations[l])+emptywells)/(num_rows_line-2*size_empty_edge))/numplates) + 2*size_empty_edge); + + +constraint forall(i in Plates, j in Rows) (controls_in_plate_row[i,j] <= control_slack + ceil(ceil((total_controls+emptywells)/numplates)/num_rows_line) ); + +constraint forall(i in Plates, j in Rows) (controls_in_plate_row[i,j] >= floor(floor((total_controls+emptywells)/numplates)/num_rows_line) ); + +%constraint if emptywells == 0 then forall(i in Plates, j in Rows) (controls_only_in_plate_row[i,j] == controls_in_plate_row[i,j] ) else forall(i in Plates, j in Rows) (controls_only_in_plate_row[i,j] <= controls_in_plate_row[i,j] ) endif; + +constraint forall(i in Plates, j in Rows)(controls_in_plate_row[i,j] == controls_only_in_plate_row[i,j] + emptywells_in_plate_row[i,j]); + +constraint forall(i in Plates, k in Columns) (controls_in_plate_column[i,k] <= control_slack + ceil(ceil((total_controls+emptywells)/numplates)/num_cols_line) ); + +constraint forall(i in Plates, k in Columns) (controls_in_plate_column[i,k] >= floor(floor((total_controls+emptywells)/numplates)/num_cols_line) ); + +constraint forall(i in Plates, k in Columns) (controls_only_in_plate_column[i,k] + emptywells_in_plate_column[i,k] == controls_in_plate_column[i,k] ); + +constraint if interconnected_plates then forall(i in Plates, j in Rows) (controls_in_plate_row[i,j] >= (controls_in_row[j] div numplates)) endif; + +constraint if interconnected_plates then forall(i in Plates, j in Rows) (controls_in_plate_row[i,j] <= ((controls_in_row[j]+numplates-1)div numplates)) endif; + + +%%% Controls per column across plates%%% + +constraint if interconnected_plates then forall(k in Columns) (controls_in_column[k] <= ceil((total_controls+emptywells)/num_cols_line) ) endif; + +constraint if interconnected_plates then forall(k in Columns) (controls_in_column[k] >= floor((total_controls+emptywells)/num_cols_line) ) endif; + +constraint forall(i in Plates, k in Columns) (controls_only_in_plate_column[i,k] <= ceil(ceil(total_controls / numplates) / num_cols_line) ); + +constraint forall(i in Plates, k in Columns) (controls_only_in_plate_column[i,k] >= ((total_controls div numplates) div num_cols_line) ); + + +%%% Attempting to balance controls inside each row + +array [Rows,Horizontal] of var 0..numplates*num_cols_line: controls_in_half_row; + +array [Rows,Horizontal] of var 0..numplates*num_cols_line: controls_only_in_half_row; %Does NOT include empty wells + +array [Plates,Rows,Horizontal] of var 0..ceil(num_cols_line/2): controls_in_half_row_plate; %includes empty wells + +array [Plates,Rows,Horizontal] of var 0..ceil(num_cols_line/2): controls_only_in_half_row_plate; %Does NOT include empty wells + +constraint forall(j in Rows) ( among(controls_in_half_row[j,left],[plates[i,j,k] | i in Plates, k in 1..floor(num_cols_line/2)], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(j in Rows) ( among(controls_in_half_row[j,right],[plates[i,j,k] | i in Plates, k in floor(num_cols_line/2)+1..num_cols_line], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(i in Plates, j in Rows) ( among(controls_in_half_row_plate[i,j,left],[plates[i,j,k] | k in 1..floor(num_cols_line/2)], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(i in Plates, j in Rows) ( among(controls_in_half_row_plate[i,j,right],[plates[i,j,k] | k in floor(num_cols_line/2)+1..num_cols_line], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(i in Plates, j in Rows) ( among(controls_only_in_half_row_plate[i,j,left],[plates[i,j,k] | k in 1..(num_cols_line div 2)], ( experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(i in Plates, j in Rows) ( among(controls_only_in_half_row_plate[i,j,right],[plates[i,j,k] | k in (num_cols_line div 2)+1..num_cols_line], (experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +% controls_only_in_half_row +constraint forall(j in Rows) ( among(controls_only_in_half_row[j,left],[plates[i,j,k] | i in Plates, k in 1..(num_cols_line div 2)], ( experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(j in Rows) ( among(controls_only_in_half_row[j,right],[plates[i,j,k] | i in Plates, k in (num_cols_line div 2)+1..num_cols_line], ( experiments+1..experiments+num_controls*max_control_concentrations))::domain); + + +constraint forall(j in Rows, h in Horizontal) (controls_in_half_row[j,h] == sum(i in Plates)(controls_in_half_row_plate[i,j,h])); + +constraint forall(i in Plates, j in Rows) (controls_only_in_plate_row[i,j] == sum(h in Horizontal)(controls_only_in_half_row_plate[i,j,h])); + +constraint forall(i in Plates, j in Rows, h in Horizontal) (controls_only_in_plate_row[i,j] div 2 <= controls_only_in_half_row_plate[i,j,h]); + +constraint forall(i in Plates, j in Rows, h in Horizontal) ((controls_only_in_plate_row[i,j]+1) div 2 >= controls_only_in_half_row_plate[i,j,h]); + +constraint forall(i in Plates, k in Columns, v in Vertical) (controls_only_in_plate_column[i,k] div 2 <= controls_only_in_half_column_plate[i,k,v]); + +constraint forall(i in Plates, k in Columns, v in Vertical) ((controls_only_in_plate_column[i,k]+1) div 2 >= controls_only_in_half_column_plate[i,k,v]); + + +%constraint forall(j in 1+size_empty_edge..num_rows_line-size_empty_edge) ( ceil((ceil((sum(l in 1..num_controls)(control_replicates[l]*control_concentrations[l])+emptywells)/(num_rows_line-2*size_empty_edge)) + 2*size_empty_edge*numplates)/2) >= controls_in_half_row[j,left]); + +%constraint forall(j in 1+size_empty_edge..num_rows_line-size_empty_edge) ( ceil((ceil((sum(l in 1..num_controls)(control_replicates[l]*control_concentrations[l])+emptywells)/(num_cols_line-2*size_empty_edge)) + 2*size_empty_edge*numplates)/2) >= controls_in_half_row[j,right]); + +constraint if interconnected_plates then forall(j in Rows) ( (floor((total_controls+emptywells)/num_rows_line) div 2) <= controls_in_half_row[j,left]) endif; + +constraint if interconnected_plates then forall(j in Rows) ( (floor((total_controls+emptywells)/num_rows_line) div 2) <= controls_in_half_row[j,right]) endif; + +constraint if interconnected_plates then forall(j in Rows, h in Horizontal) ( (controls_in_row[j] div 2) <= controls_in_half_row[j,h] ) endif; +% constraint if interconnected_plates then forall(j in Rows, h in Horizontal) ( ((controls_in_row[j]+1) div 2) >= controls_in_half_row[j,h] ) endif; + +constraint if interconnected_plates then forall(j in Rows) (controls_in_row[j] == (controls_in_half_row[j,left] + controls_in_half_row[j,right])) endif; + + +%constraint if temp_bool then forall(i in Plates, h in Horizontal, j in 1+size_empty_edge..num_rows_line-size_empty_edge) ( floor(floor((floor((sum(l in 1..num_controls)(control_replicates[l]*control_concentrations[l])+emptywells)/numplates)/(num_rows_line-2*size_empty_edge))/2) + size_empty_edge) <= controls_in_half_row_plate[i,j,h]) endif; + +constraint forall(i in Plates, h in Horizontal, j in Rows) ( (controls_in_plate_row[i,j] div 2) <= controls_in_half_row_plate[i,j,h]); + +constraint forall(i in Plates, h in Horizontal, j in Rows) ( ( (controls_in_plate_row[i,j]+1) div 2) >= controls_in_half_row_plate[i,j,h]); + + +%constraint if temp_bool then forall(i in Plates, h in Horizontal, j in 1+size_empty_edge..num_rows_line-size_empty_edge) ( ceil(ceil((ceil((sum(l in 1..num_controls)(control_replicates[l]*control_concentrations[l])+emptywells)/numplates)/(num_rows_line-2*size_empty_edge))/2) + size_empty_edge) >= controls_in_half_row_plate[i,j,h]) endif; + + +constraint forall(i in Plates, j in Rows) (controls_in_plate_row[i,j] == controls_in_half_row_plate[i,j,left] + controls_in_half_row_plate[i,j,right]); + + +%%% Attempting to balance controls inside each column + +array [Columns,Vertical] of var 0..numplates*num_rows_line: controls_in_half_column; + +array [Columns,Vertical] of var 0..numplates*num_rows_line: controls_only_in_half_column; %Does NOT include empty wells + +array [Plates,Columns,Vertical] of var 0..(num_rows_line div 2): controls_in_half_column_plate; %includes empty wells + +array [Plates,Columns,Vertical] of var 0..(num_rows_line div 2): controls_only_in_half_column_plate; %Does NOT include empty wells + +array [Plates,Rows,Horizontal] of var 0..ceil(num_cols_line/2): experiments_in_half_row_plate; + +array [Plates,Columns,Vertical] of var 0..ceil(num_rows_line/2): experiments_in_half_column_plate; + +constraint forall(k in Columns) ( among(controls_in_half_column[k,upper],[plates[i,j,k] | i in Plates, j in 1..(num_rows_line div 2)], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(k in Columns) ( among(controls_in_half_column[k,lower],[plates[i,j,k] | i in Plates, j in (num_rows_line div 2)+1..num_rows_line], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(i in Plates, k in Columns) ( among(controls_in_half_column_plate[i,k,upper],[plates[i,j,k] | j in 1..(num_rows_line div 2)], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(i in Plates, k in Columns) ( among(controls_in_half_column_plate[i,k,lower],[plates[i,j,k] | j in (num_rows_line div 2)+1..num_rows_line], ({0} union experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(i in Plates, k in Columns) ( among(controls_only_in_half_column_plate[i,k,upper],[plates[i,j,k] | j in 1..(num_rows_line div 2)], ( experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(i in Plates, k in Columns) ( among(controls_only_in_half_column_plate[i,k,lower],[plates[i,j,k] | j in (num_rows_line div 2)+1..num_rows_line], ( experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +%% controls_only_in_half_column +constraint forall(k in Columns) ( among(controls_only_in_half_column[k,upper],[plates[i,j,k] | i in Plates, j in 1..(num_rows_line div 2)], ( experiments+1..experiments+num_controls*max_control_concentrations))::domain); + +constraint forall(k in Columns) ( among(controls_only_in_half_column[k,lower],[plates[i,j,k] | i in Plates, j in (num_rows_line div 2)+1..num_rows_line], ( experiments+1..experiments+num_controls*max_control_concentrations))::domain); + + +constraint forall(k in Columns, v in Vertical) (controls_in_half_column[k,v] == sum(i in Plates)(controls_in_half_column_plate[i,k,v])); + +constraint forall(i in Plates, k in Columns) ( controls_in_plate_column[i,k] == sum(v in Vertical)(controls_in_half_column_plate[i,k,v])); + +constraint forall(i in Plates, v in Vertical, k in Columns) ( controls_in_plate_column[i,k] <= 2*controls_in_half_column_plate[i,k,v] + 1); + +constraint forall(i in Plates, v in Vertical, k in Columns) ( controls_in_plate_column[i,k] >= 2*controls_in_half_column_plate[i,k,v] - 1); + +constraint forall(i in Plates, v in Vertical, k in Columns) ( (controls_in_plate_column[i,k] div 2) <= controls_in_half_column_plate[i,k,v]); + +constraint forall(i in Plates, v in Vertical, k in Columns) ( ( (controls_in_plate_column[i,k]+1) div 2) >= controls_in_half_column_plate[i,k,v]); + +constraint forall(i in Plates, k in Columns) ( controls_only_in_plate_column[i,k] == sum(v in Vertical)(controls_only_in_half_column_plate[i,k,v])); + +%constraint forall(j in 1+size_empty_edge..num_rows_line-size_empty_edge) ( ceil((ceil((sum(l in 1..num_controls)(control_replicates[l]*control_concentrations[l])+emptywells)/(num_rows_line-2*size_empty_edge)) + 2*size_empty_edge*numplates)/2) >= controls_in_half_row[j,left]); + +%constraint forall(j in 1+size_empty_edge..num_rows_line-size_empty_edge) ( ceil((ceil((sum(l in 1..num_controls)(control_replicates[l]*control_concentrations[l])+emptywells)/(num_cols_line-2*size_empty_edge)) + 2*size_empty_edge*numplates)/2) >= controls_in_half_row[j,right]); + +constraint if interconnected_plates then forall(k in Columns, v in Vertical) ( (floor((total_controls+emptywells)/num_cols_line) div 2) <= controls_in_half_column[k,v]) endif; + +constraint if interconnected_plates then forall(k in Columns) (controls_in_column[k] == (controls_in_half_column[k,upper] + controls_in_half_column[k,lower])) endif; + +%%% Defining experiments_in_half_row_plate /\ experiments_in_half_column_plate + +constraint forall(i in Plates, j in Rows) ( among(experiments_in_half_row_plate[i,j,left], [plates[i,j,k] | k in 1..(num_cols_line div 2)], (1..experiments))::domain); +constraint forall(i in Plates, j in Rows) ( among(experiments_in_half_row_plate[i,j,right],[plates[i,j,k] | k in (num_cols_line div 2)+1..num_cols_line], (1..experiments))::domain); + +constraint forall(i in Plates, k in Columns) ( among(experiments_in_half_column_plate[i,k,upper],[plates[i,j,k] | j in 1..(num_rows_line div 2)], (1..experiments))::domain); +constraint forall(i in Plates, k in Columns) ( among(experiments_in_half_column_plate[i,k,lower],[plates[i,j,k] | j in (num_rows_line div 2)+1..num_rows_line], (1..experiments))::domain); + +constraint forall(i in Plates, j in Rows) ( sum(h in Horizontal)(experiments_in_half_row_plate[i,j,h]) == experiments_in_plate_row[i,j]); +constraint forall(i in Plates, k in Columns) ( sum(v in Vertical)(experiments_in_half_column_plate[i,k,v]) == experiments_in_plate_column[i,k]); + +constraint forall(i in Plates, j in Rows, h in Horizontal) ( experiments_in_half_row_plate[i,j,h] >= (experiments_in_plate_row[i,j] div 2)); %%WARNING +constraint forall(i in Plates, j in Rows, h in Horizontal) ( experiments_in_half_row_plate[i,j,h] <= ((experiments_in_plate_row[i,j] + 1) div 2)); %%WARNING + +constraint forall(i in Plates, k in Columns, v in Vertical) (experiments_in_half_column_plate[i,k,v] >= (experiments_in_plate_column[i,k] div 2)); %%WARNING +constraint forall(i in Plates, k in Columns, v in Vertical) (experiments_in_half_column_plate[i,k,v] <= ((experiments_in_plate_column[i,k] + 1) div 2)); %%WARNING + +constraint forall(i in Plates, h in Horizontal, j in Rows) ( (num_cols_line div 2) == controls_in_half_row_plate[i,j,h] + experiments_in_half_row_plate[i,j,h]); +constraint forall(i in Plates, v in Vertical, k in Columns) ( (num_rows_line div 2) == controls_in_half_column_plate[i,k,v] + experiments_in_half_column_plate[i,k,v]); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Balancing each type of control per row + +constraint forall(i in Plates, j in Rows, c in 1..num_controls*max_control_concentrations) (count_eq(plates[i,j,..],experiments+c,controls_row_plates[i,j,c])); + +constraint forall(i in Plates, j in Rows) (count_eq(plates[i,j,..],0,controls_row_plates[i,j,0])); + +constraint forall(i in Plates, j in Rows, c in 1..num_controls*max_control_concentrations) (controls_row_plates[i,j,c] <= ceil(max_controls_per_plate[c]/num_rows_line)); + +constraint forall(i in Plates, j in Rows, c in 1..num_controls*max_control_concentrations) (controls_row_plates[i,j,c] >= (min_controls_per_plate[c] div num_rows_line)); + +constraint forall(i in Plates, j in Rows) (sum(controls_row_plates[i,j,..]) + experiments_in_plate_row[i,j] == num_cols_line); + +constraint sum(controls_row_plates[..,..,0]) == emptywells; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Balancing each type of control per column + +constraint forall(i in Plates, k in Columns, c in 1..num_controls*max_control_concentrations) (count_eq(plates[i,..,k],experiments+c,controls_column_plates[i,k,c])); + +constraint forall(i in Plates, k in Columns) (count_eq(plates[i,..,k],0,controls_column_plates[i,k,0])); + +constraint forall(i in Plates, k in Columns, c in 1..num_controls*max_control_concentrations) (controls_column_plates[i,k,c] <= ceil(max_controls_per_plate[c]/num_cols_line)); + +constraint forall(i in Plates, k in Columns, c in 1..num_controls*max_control_concentrations) (controls_column_plates[i,k,c] >= (min_controls_per_plate[c] div num_cols_line)); + +constraint forall(i in Plates, k in Columns) (sum(controls_column_plates[i,k,..]) + experiments_in_plate_column[i,k] == num_rows_line); + +constraint sum(controls_column_plates[..,..,0]) == emptywells; + +%% Implied all_different constraint: all experiments/compounds are different +constraint alldifferent_except(plates, {0} union experiments+1..experiments+num_controls*max_control_concentrations):: domain; + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% Distributing controls + +%% Controls of the same type are not located next to each other + +%% Implied constraint +constraint if (min(spread_control++[1])>0) then forall(i in Plates, j in 1..num_rows_line-1, k in 1..num_cols_line-1) (alldifferent_except_0([plates[i,j,k],plates[i,j+1,k],plates[i,j,k+1],plates[i,j+1,k+1]])):: domain endif; + +array [Plates,Rows,Columns] of var 0..1: controls_layout; + +constraint forall(i in Plates, j in Rows, k in Columns)(controls_layout[i,j,k] == ((plates[i,j,k] > experiments)) ); + + +constraint if spread_controls then forall(i in Plates, j in Rows)(regular(controls_layout[i,j,..], "(0|(1 0))* (1?)")::domain) endif; + +constraint if spread_controls then forall(i in Plates, k in Columns)(regular(controls_layout[i,..,k], "(0|(1 0))* (1?)")::domain) endif; + +constraint if force_spread_controls then forall(i in Plates, j in Rows)(regular(controls_layout[i,j,..], "(0|(1 0 0))* (1|(1 0))?")::domain) endif; + +constraint if force_spread_controls then forall(i in Plates, k in Columns)(regular(controls_layout[i,..,k], "(0|(1 0 0))* (1|(1 0))?")::domain) endif; + + +constraint if spread_controls then forall(i in Plates, j in Rows)(my_implied_cost(controls_only_in_plate_row[i,j],controls_layout[i,j,..])::domain) endif; + +constraint if spread_controls /\ num_cols_line > 2 then forall(i in Plates, j in Rows)(my_implied_cost(controls_only_in_half_row_plate[i,j,left],controls_layout[i,j,1..(num_cols_line div 2)])::domain) endif; + +constraint if spread_controls /\ num_cols_line > 2 then forall(i in Plates, j in Rows)(my_implied_cost(controls_only_in_half_row_plate[i,j,right],[controls_layout[i,j,k] | k in ((num_cols_line div 2)+1)..num_cols_line])::domain) endif; + +constraint if spread_controls /\ num_cols_line > 2 then forall(i in Plates, j in Rows)(my_implied_cost(controls_only_in_half_row_plate[i,j,right],reverse([controls_layout[i,j,k] | k in ((num_cols_line div 2)+1)..num_cols_line]))::domain) endif; + +constraint if spread_controls then forall(i in Plates, k in Columns)(my_implied_cost(controls_only_in_plate_column[i,k],controls_layout[i,..,k])::domain) endif; + +constraint if spread_controls /\ num_rows_line > 2 then forall(i in Plates, k in Columns)(my_implied_cost(controls_only_in_half_column_plate[i,k,upper],[controls_layout[i,j,k] | j in 1..(num_rows_line div 2)])::domain) endif; + +constraint if spread_controls /\ num_rows_line > 2 then forall(i in Plates, k in Columns)(my_implied_cost(controls_only_in_half_column_plate[i,k,lower],[controls_layout[i,j,k] | j in (num_rows_line div 2)+1..num_rows_line])::domain) endif; + +constraint if force_spread_controls then forall(i in Plates, j in Rows)(my_implied_cost_ext(controls_only_in_plate_row[i,j],controls_layout[i,j,..])::domain) endif; + +constraint if force_spread_controls then forall(i in Plates, k in Columns)(my_implied_cost_ext(controls_only_in_plate_column[i,k],controls_layout[i,..,k])::domain) endif; + + +constraint forall(i in Plates) (sum_controls_ul_plates[i,upper] == sum(controls_layout[i,1..(num_rows_line div 2),..])); +constraint forall(i in Plates) (sum_controls_ul_plates[i,lower] == sum(controls_layout[i,(num_rows_line div 2)+1..num_rows_line,..])); + +constraint forall(i in Plates) (sum_controls_lr_plates[i,left] == sum(controls_layout[i,..,1..(num_cols_line div 2)])); +constraint forall(i in Plates) (sum_controls_lr_plates[i,right] == sum(controls_layout[i,..,(num_cols_line div 2)+1..num_cols_line])); + + + +constraint forall(i in Plates, j in Rows) (global_cardinality(controls_layout[i,j,..], + [0,1], + [num_cols_line-controls_only_in_plate_row[i,j],controls_only_in_plate_row[i,j]])); + +constraint forall(i in Plates, k in Columns) (global_cardinality(controls_layout[i,..,k], + [0,1], + [num_rows_line-controls_only_in_plate_column[i,k],controls_only_in_plate_column[i,k]])); + + + +constraint forall(ctr in 1..num_controls*max_control_concentrations)( if spread_control[ctr] then + forall(i in Plates, j in 2..num_rows_line-1, k in 2..num_cols_line-1)(plates[i,j,k] == experiments+ctr -> + plates[i,j-1,k-1]!= plates[i,j,k] /\ + plates[i,j-1,k] != plates[i,j,k] /\ + plates[i,j-1,k+1]!= plates[i,j,k] /\ + plates[i,j,k-1] != plates[i,j,k] /\ + plates[i,j,k+1] != plates[i,j,k] /\ + plates[i,j+1,k-1]!= plates[i,j,k] /\ + plates[i,j+1,k] != plates[i,j,k] /\ + plates[i,j+1,k+1]!= plates[i,j,k] + ) + %% This might be leading to unsat when there are too many of the same control +% else forall(i in Plates, j in 2..num_rows_line-1, k in 2..num_cols_line-1)(plates[i,j,k] == experiments+ctr -> + % plates[i,j-1,k] != plates[i,j,k] /\ + % plates[i,j,k-1] != plates[i,j,k] /\ + % plates[i,j,k+1] != plates[i,j,k] /\ + % plates[i,j+1,k] != plates[i,j,k] + % ) + endif ); + +%% Upper row j=1 +constraint forall(ctr in 1..num_controls*max_control_concentrations)( if spread_control[ctr] then + forall(i in Plates, k in 2..num_cols_line-1)(plates[i,1,k] == experiments+ctr -> + plates[i,1,k-1] != plates[i,1,k] /\ + plates[i,1,k+1] != plates[i,1,k] /\ + plates[i,2,k-1]!= plates[i,1,k] /\ + plates[i,2,k] != plates[i,1,k] /\ + plates[i,2,k+1]!= plates[i,1,k] + ) + %% This might be leading to unsat when there are too many of the same control +% else forall(i in Plates, k in 2..num_cols_line-1)(plates[i,1,k] == experiments+ctr -> + % plates[i,1,k-1] != plates[i,1,k] /\ + % plates[i,1,k+1] != plates[i,1,k] /\ + % plates[i,2,k] != plates[i,1,k] + % ) + endif ); + +%% Lower row j = num_rows_line +constraint forall(ctr in 1..num_controls*max_control_concentrations)( if spread_control[ctr] then + forall(i in Plates, k in 2..num_cols_line-1)(plates[i,num_rows_line,k] == experiments+ctr -> + plates[i,num_rows_line-1,k-1]!= plates[i,num_rows_line,k] /\ + plates[i,num_rows_line-1,k] != plates[i,num_rows_line,k] /\ + plates[i,num_rows_line-1,k+1]!= plates[i,num_rows_line,k] /\ + plates[i,num_rows_line,k-1] != plates[i,num_rows_line,k] /\ + plates[i,num_rows_line,k+1] != plates[i,num_rows_line,k] + ) + %% This might be leading to unsat when there are too many of the same control +% else forall(i in Plates, k in 2..num_cols_line-1)(plates[i,num_rows_line,k] == experiments+ctr -> + % plates[i,num_rows_line-1,k] != plates[i,num_rows_line,k] /\ + % plates[i,num_rows_line,k-1] != plates[i,num_rows_line,k] /\ + % plates[i,num_rows_line,k+1] != plates[i,num_rows_line,k] + % ) + endif ); + +%% Left column k=1 +constraint forall(ctr in 1..num_controls*max_control_concentrations)( if spread_control[ctr] then + forall(i in Plates, j in 2..num_rows_line-1)(plates[i,j,1] == experiments+ctr -> + plates[i,j-1,1] != plates[i,j,1] /\ + plates[i,j-1,2]!= plates[i,j,1] /\ + plates[i,j,2] != plates[i,j,1] /\ + plates[i,j+1,1] != plates[i,j,1] /\ + plates[i,j+1,2]!= plates[i,j,1] + ) + %% This might be leading to unsat when there are too many of the same control + % else forall(i in Plates, j in 2..num_rows_line-1)(plates[i,j,1] == experiments+ctr -> + % plates[i,j-1,1] != plates[i,j,1] /\ + % plates[i,j,2] != plates[i,j,1] /\ + % plates[i,j+1,1] != plates[i,j,1] + % ) + endif ); + +%% Right column k=num_cols_line +constraint forall(ctr in 1..num_controls*max_control_concentrations)( if spread_control[ctr] then + forall(i in Plates, j in 2..num_rows_line-1)(plates[i,j,num_cols_line] == experiments+ctr -> + plates[i,j-1,num_cols_line-1]!= plates[i,j,num_cols_line] /\ + plates[i,j-1,num_cols_line] != plates[i,j,num_cols_line] /\ + plates[i,j,num_cols_line-1] != plates[i,j,num_cols_line] /\ + plates[i,j+1,num_cols_line-1]!= plates[i,j,num_cols_line] /\ + plates[i,j+1,num_cols_line] != plates[i,j,num_cols_line] + ) + %% This might be leading to unsat when there are too many of the same control +% else forall(i in Plates, j in 2..num_rows_line-1)(plates[i,j,num_cols_line] == experiments+ctr -> + % plates[i,j-1,num_cols_line] != plates[i,j,num_cols_line] /\ + % plates[i,j,num_cols_line-1] != plates[i,j,num_cols_line] /\ + % plates[i,j+1,num_cols_line] != plates[i,j,num_cols_line] + % ) + endif ); + +%% Upper left corner j=1 k=1 +constraint forall(ctr in 1..num_controls*max_control_concentrations)( if spread_control[ctr] then + forall(i in Plates)(plates[i,1,1] == experiments+ctr -> + plates[i,1,2] != plates[i,1,1] /\ + plates[i,2,1] != plates[i,1,1] /\ + plates[i,2,2] != plates[i,1,1] + ) + %% This might be leading to unsat when there are too many of the same control +% else forall(i in Plates)(plates[i,1,1] == experiments+ctr -> + % plates[i,1,2] != plates[i,1,1] /\ + % plates[i,2,1] != plates[i,1,1] + % ) + endif ); + + +%% Upper right corner j=1 k=num_cols_line +constraint forall(ctr in 1..num_controls*max_control_concentrations)( if spread_control[ctr] then + forall(i in Plates)(plates[i,1,num_cols_line] == experiments+ctr -> + plates[i,1,num_cols_line-1] != plates[i,1,num_cols_line] /\ + plates[i,2,num_cols_line-1] != plates[i,1,num_cols_line] /\ + plates[i,2,num_cols_line] != plates[i,1,num_cols_line] + ) + %% This might be leading to unsat when there are too many of the same control +% else forall(i in Plates)(plates[i,1,num_cols_line] == experiments+ctr -> + % plates[i,1,num_cols_line-1] != plates[i,1,num_cols_line] /\ + % plates[i,2,num_cols_line] != plates[i,1,num_cols_line] + % ) + endif ); + + +%% Lower left corner j = num_rows_line k=1 +constraint forall(ctr in 1..num_controls*max_control_concentrations)( if spread_control[ctr] then + forall(i in Plates)(plates[i,num_rows_line,1] == experiments+ctr -> + plates[i,num_rows_line-1,1] != plates[i,num_rows_line,1] /\ + plates[i,num_rows_line-1,2]!= plates[i,num_rows_line,1] /\ + plates[i,num_rows_line,2] != plates[i,num_rows_line,1] + ) + %% This might be leading to unsat when there are too many of the same control +% else forall(i in Plates)(plates[i,num_rows_line,1] == experiments+ctr -> + % plates[i,num_rows_line-1,1] != plates[i,num_rows_line,1] /\ + % plates[i,num_rows_line,2] != plates[i,num_rows_line,1] + % ) + endif ); + +%% Lower right corner j = num_rows_line k = num_cols_line +constraint forall(ctr in 1..num_controls*max_control_concentrations)( if spread_control[ctr] then + forall(i in Plates, k in 2..num_cols_line-1)(plates[i,num_rows_line,num_cols_line] == experiments+ctr -> + plates[i,num_rows_line-1,num_cols_line-1]!= plates[i,num_rows_line,num_cols_line] /\ + plates[i,num_rows_line-1,num_cols_line] != plates[i,num_rows_line,num_cols_line] /\ + plates[i,num_rows_line,num_cols_line-1] != plates[i,num_rows_line,num_cols_line] + ) + %% This might be leading to unsat when there are too many of the same control +% else forall(i in Plates)(plates[i,num_rows_line,num_cols_line] == experiments+ctr -> + % plates[i,num_rows_line-1,num_cols_line] != plates[i,num_rows_line,num_cols_line] /\ + % plates[i,num_rows_line,num_cols_line-1] != plates[i,num_rows_line,num_cols_line] + % ) + endif ); + + +constraint if spread_controls then forall(i in Plates, j in 2..num_rows_line-1, k in 2..num_cols_line-1)(plates[i,j,k] > experiments -> + %plates[i,j-1,k-1] <= experiments /\ + plates[i,j-1,k] <= experiments /\ + %plates[i,j-1,k+1] <= experiments /\ + plates[i,j,k-1] <= experiments /\ + plates[i,j,k+1] <= experiments /\ + %plates[i,j+1,k-1] <= experiments /\ + plates[i,j+1,k] <= experiments %/\ + %plates[i,j+1,k+1] <= experiments + ) endif; + +constraint if spread_controls then forall(i in Plates, j in 2..num_rows_line-1)(plates[i,j,1] > experiments -> + plates[i,j-1,1] <= experiments /\ + plates[i,j,2] <= experiments /\ + plates[i,j+1,1] <= experiments + ) endif; + +constraint if spread_controls then forall(i in Plates, j in 2..num_rows_line-1)(plates[i,j,num_cols_line] > experiments -> + plates[i,j-1,num_cols_line] <= experiments /\ + plates[i,j,num_cols_line-1] <= experiments /\ + plates[i,j+1,num_cols_line] <= experiments + ) endif; + + +constraint if spread_controls then forall(i in Plates, k in 2..num_cols_line-1)(plates[i,1,k] > experiments -> + plates[i,1,k-1] <= experiments /\ + plates[i,1,k+1] <= experiments /\ + plates[i,2,k] <= experiments + ) endif; + +constraint if spread_controls then forall(i in Plates, k in 2..num_cols_line-1)(plates[i,num_rows_line,k] > experiments -> + plates[i,num_rows_line-1,k] <= experiments /\ + plates[i,num_rows_line,k-1] <= experiments /\ + plates[i,num_rows_line,k+1] <= experiments + ) endif; + + +constraint if spread_controls then forall(i in Plates)(plates[i,1,1] > experiments -> + plates[i,1,2] <= experiments /\ + plates[i,2,1] <= experiments + ) endif; + +constraint if spread_controls then forall(i in Plates)(plates[i,1,num_cols_line] > experiments -> + plates[i,1,num_cols_line-1] <= experiments /\ + plates[i,2,num_cols_line] <= experiments + ) endif; + +constraint if spread_controls then forall(i in Plates)(plates[i,num_rows_line,1] > experiments -> + plates[i,num_rows_line-1,1] <= experiments /\ + plates[i,num_rows_line,2] <= experiments + ) endif; + +constraint if spread_controls then forall(i in Plates)(plates[i,num_rows_line,num_cols_line] > experiments -> + plates[i,num_rows_line-1,num_cols_line] <= experiments /\ + plates[i,num_rows_line,num_cols_line-1] <= experiments + ) endif; +%% TODO: Add spread_controls contraint in borders! + + +% j=1 +constraint if spread_controls then forall(i in Plates, k in 2..num_cols_line-1)(plates[i,1,k] > experiments -> + plates[i,1,k-1] != plates[i,1,k] /\ + plates[i,1,k+1] != plates[i,1,k] /\ + plates[i,2,k-1]!= plates[i,1,k] /\ + plates[i,2,k] != plates[i,1,k] /\ + plates[i,2,k+1]!= plates[i,1,k] + ) endif; + +% j = num_rows_line +constraint if spread_controls then forall(i in Plates, k in 2..num_cols_line-1)(plates[i,num_rows_line,k] > experiments -> + plates[i,num_rows_line-1,k-1]!= plates[i,num_rows_line,k] /\ + plates[i,num_rows_line-1,k] != plates[i,num_rows_line,k] /\ + plates[i,num_rows_line-1,k+1]!= plates[i,num_rows_line,k] /\ + plates[i,num_rows_line,k-1] != plates[i,num_rows_line,k] /\ + plates[i,num_rows_line,k+1] != plates[i,num_rows_line,k] + ) endif; + +% First column: k = 1 +constraint if spread_controls then forall(i in Plates, j in 2..num_rows_line-1)(plates[i,j,1] > experiments -> + plates[i,j-1,1] != plates[i,j,1] /\ + plates[i,j-1,2]!= plates[i,j,1] /\ + plates[i,j,2] != plates[i,j,1] /\ + plates[i,j+1,1] != plates[i,j,1] /\ + plates[i,j+1,2]!= plates[i,j,1] + ) endif; + +% Last column: k = num_cols_line +constraint if spread_controls then forall(i in Plates, j in 2..num_rows_line-1)(plates[i,j,num_cols_line] > experiments -> + plates[i,j-1,num_cols_line-1]!= plates[i,j,num_cols_line] /\ + plates[i,j-1,num_cols_line] != plates[i,j,num_cols_line] /\ + plates[i,j,num_cols_line-1] != plates[i,j,num_cols_line] /\ + plates[i,j+1,num_cols_line-1]!= plates[i,j,num_cols_line] /\ + plates[i,j+1,num_cols_line] != plates[i,j,num_cols_line] + ) endif; + + +%% All especial cases for the corners are now covered by the alldifferent_except_0 constraints + + +constraint if force_spread_controls then (forall(i in Plates, j in 2..num_rows_line-1, k in 2..num_cols_line-1)(plates[i,j,k] > experiments -> + plates[i,j-1,k] <= experiments /\ + plates[i,j,k-1] <= experiments /\ + plates[i,j,k+1] <= experiments /\ + plates[i,j+1,k] <= experiments + ) + ) endif ; + +constraint if force_spread_controls then forall(i in Plates, j in 2..num_rows_line-1, k in 2..num_cols_line-1)(plates[i,j,k] > experiments -> + plates[i,j-1,k-1] <= experiments /\ + plates[i,j-1,k] <= experiments /\ + plates[i,j-1,k+1] <= experiments /\ + plates[i,j,k-1] <= experiments /\ + plates[i,j,k+1] <= experiments /\ + plates[i,j+1,k-1] <= experiments /\ + plates[i,j+1,k] <= experiments /\ + plates[i,j+1,k+1] <= experiments + ) endif; + +constraint if force_spread_controls then forall(i in Plates, j in 3..num_rows_line-2, k in 3..num_cols_line-2)(plates[i,j,k] > experiments -> + plates[i,j-2,k-2] <= experiments /\ + plates[i,j-2,k-1] <= experiments /\ + plates[i,j-2,k] <= experiments /\ + plates[i,j-2,k+1] <= experiments /\ + plates[i,j-2,k+2] <= experiments /\ + plates[i,j-1,k-2] <= experiments /\ + plates[i,j-1,k+2] <= experiments /\ + plates[i,j,k-2] <= experiments /\ + plates[i,j,k+2] <= experiments /\ + plates[i,j+1,k-2] <= experiments /\ + plates[i,j+1,k+2] <= experiments /\ + plates[i,j+2,k-2] <= experiments /\ + plates[i,j+2,k-1] <= experiments /\ + plates[i,j+2,k] <= experiments /\ + plates[i,j+2,k+1] <= experiments /\ + plates[i,j+2,k+2] <= experiments + ) endif; + + +constraint if force_spread_controls then forall(i in Plates, j = 2, k in 3..num_cols_line-2)(plates[i,j,k] > experiments -> + plates[i,j,k-2] <= experiments /\ + plates[i,j,k+2] <= experiments /\ + plates[i,j-1,k-2] <= experiments /\ + plates[i,j-1,k+2] <= experiments /\ + plates[i,j+1,k-2] <= experiments /\ + plates[i,j+1,k+2] <= experiments /\ + plates[i,j+2,k-2] <= experiments /\ + plates[i,j+2,k-1] <= experiments /\ + plates[i,j+2,k] <= experiments /\ + plates[i,j+2,k+1] <= experiments /\ + plates[i,j+2,k+2] <= experiments + ) endif; + + +constraint if force_spread_controls then forall(i in Plates, j = 1, k in 3..num_cols_line-2)(plates[i,j,k] > experiments -> + plates[i,j,k-2] <= experiments /\ + plates[i,j,k+2] <= experiments /\ + plates[i,j+1,k-2] <= experiments /\ + plates[i,j+1,k+2] <= experiments /\ + plates[i,j+2,k-2] <= experiments /\ + plates[i,j+2,k-1] <= experiments /\ + plates[i,j+2,k] <= experiments /\ + plates[i,j+2,k+1] <= experiments /\ + plates[i,j+2,k+2] <= experiments + ) endif; + + +constraint if force_spread_controls then forall(i in Plates, j = num_rows_line-1, k in 3..num_cols_line-2)(plates[i,j,k] > experiments -> + plates[i,j-2,k-2] <= experiments /\ + plates[i,j-2,k-1] <= experiments /\ + plates[i,j-2,k] <= experiments /\ + plates[i,j-2,k+1] <= experiments /\ + plates[i,j-2,k+2] <= experiments /\ + plates[i,j,k-2] <= experiments /\ + plates[i,j,k+2] <= experiments /\ + plates[i,j+1,k-2] <= experiments /\ + plates[i,j+1,k+2] <= experiments + ) endif; + +constraint if force_spread_controls then forall(i in Plates, j = num_rows_line, k in 3..num_cols_line-2)(plates[i,j,k] > experiments -> + plates[i,j-2,k-2] <= experiments /\ + plates[i,j-2,k-1] <= experiments /\ + plates[i,j-2,k] <= experiments /\ + plates[i,j-2,k+1] <= experiments /\ + plates[i,j-2,k+2] <= experiments /\ + plates[i,j,k-2] <= experiments /\ + plates[i,j,k+2] <= experiments + ) endif; + + +%% Border cases for k: k=1 +constraint if force_spread_controls then forall(i in Plates, j in 3..num_rows_line-2, k=1)(plates[i,j,1] > experiments -> + plates[i,j-2,k] <= experiments /\ + plates[i,j-2,k+1] <= experiments /\ + plates[i,j-2,k+2] <= experiments /\ + plates[i,j-1,k+2] <= experiments /\ + plates[i,j,k+2] <= experiments /\ + plates[i,j+1,k+2] <= experiments /\ + plates[i,j+2,k] <= experiments /\ + plates[i,j+2,k+1] <= experiments /\ + plates[i,j+2,k+2] <= experiments + ) endif; + + +%% Border cases for k +constraint if force_spread_controls then forall(i in Plates, j in 3..num_rows_line-2, k = 2)(plates[i,j,k] > experiments -> + plates[i,j-2,k-1] <= experiments /\ + plates[i,j-2,k] <= experiments /\ + plates[i,j-2,k+1] <= experiments /\ + plates[i,j-2,k+2] <= experiments /\ + plates[i,j-1,k+2] <= experiments /\ + plates[i,j,k+2] <= experiments /\ + plates[i,j+1,k+2] <= experiments /\ + plates[i,j+2,k-1] <= experiments /\ + plates[i,j+2,k] <= experiments /\ + plates[i,j+2,k+1] <= experiments /\ + plates[i,j+2,k+2] <= experiments + ) endif; + +%% Border cases for k +constraint if force_spread_controls then forall(i in Plates, j in 3..num_rows_line-2, k = num_cols_line-1)(plates[i,j,k] > experiments -> + plates[i,j-2,k-2] <= experiments /\ + plates[i,j-2,k-1] <= experiments /\ + plates[i,j-2,k] <= experiments /\ + plates[i,j-2,k+1] <= experiments /\ + plates[i,j-1,k-2] <= experiments /\ + plates[i,j,k-2] <= experiments /\ + plates[i,j+1,k-2] <= experiments /\ + plates[i,j+2,k-2] <= experiments /\ + plates[i,j+2,k-1] <= experiments /\ + plates[i,j+2,k] <= experiments /\ + plates[i,j+2,k+1] <= experiments + ) endif; + +%% Border cases for k +constraint if force_spread_controls then forall(i in Plates, j in 3..num_rows_line-2, k = num_cols_line)(plates[i,j,k] > experiments -> + plates[i,j-2,k-2] <= experiments /\ + plates[i,j-2,k-1] <= experiments /\ + plates[i,j-2,k] <= experiments /\ + plates[i,j-1,k-2] <= experiments /\ + plates[i,j,k-2] <= experiments /\ + plates[i,j+1,k-2] <= experiments /\ + plates[i,j+2,k-2] <= experiments /\ + plates[i,j+2,k-1] <= experiments /\ + plates[i,j+2,k] <= experiments + ) endif; + + +%%%% Find a solution! %%%% + +array [Plates,1..experiments] of var {0} union Rows: experiment_plate_row; + +constraint forall(l in 1..experiments)((experiment_plate[l] > 0) -> (experiment_plate_row[experiment_plate[l],l] == experiment_row[l])); + +constraint forall(i in Plates, l in 1..experiments)((experiment_plate[l] != i) -> (experiment_plate_row[i,l] == 0)); + + + +%array [1..numplates*2] of ann: search_order = [if (i mod 2)==1 then int_search(experiment_plate_row[(i+1)div 2,..], first_fail, indomain_random) else int_search(plates[i div 2,..,..], random, indomain_max) endif | i in 1..numplates*2]; + +%array [1..numplates*2] of ann: search_order = [if (i mod 2)==1 then int_search(experiment_plate_row[(i+1)div 2,..], first_fail, indomain_random) else int_search(plates[i div 2,..,..], random, indomain_max) endif | i in 1..numplates*2]; + +%% Strategy 1: good for example04 and example11 BACKUP, previous version! +%array [1..numplates*3] of ann: search_order_a = [if (i mod 3)==1 then int_search(experiment_plate_row[(i+2) div 3,..], first_fail, indomain_random) elseif (i mod 3) == 2 then int_search(controls_layout[(i+2) div 3,..,..], random, indomain_max) else int_search(plates[i div 3,..,..], first_fail, indomain_random) endif | i in 1..numplates*3]; + +array [1..numplates*3] of ann: search_order_a = [if (i mod 3)==2 then int_search(experiment_plate_row[(i+2) div 3,..], first_fail, indomain_random) elseif (i mod 3) == 1 then int_search(controls_layout[(i+2) div 3,..,..], random, indomain_max) else int_search(plates[i div 3,..,..], first_fail, indomain_random) endif | i in 1..numplates*3]; + + +%% The one I normally use?? +%array [1..numplates*3] of ann: search_order = [if (i mod 3)==2 then int_search(experiment_plate_row[(i+2)div 3,..], first_fail, indomain_random) elseif (i mod 3) == 1 then int_search(controls_layout[(i+2)div 3,..,..], random, indomain_max) else int_search(plates[(i+2) div 3,..,..], first_fail, indomain_random) endif | i in 1..numplates*3]; + +%% Strategy 2: +array [1..numplates*2] of ann: search_order_b = [if (i mod 2) == 1 then int_search(controls_layout[(i+1) div 2,..,..], random, indomain_max) else int_search(plates[(i+1) div 2,..,..], first_fail, indomain_random) endif | i in 1..numplates*2]; + + +%solve::seq_search([int_search(experiment_plate, first_fail, indomain_random)]++ +solve::seq_search((if sorted_compounds == true then [int_search(experiment_plate, first_fail, indomain_min)] else [int_search(experiment_plate, first_fail, indomain_random)] endif)++ +(if (swap_search mod 2 == 0) then search_order_b else search_order_a endif)) +%[int_search(controls_layout, random, indomain_max), %int_search(experiment_plate_row, first_fail, indomain_random), +%int_search(plates, first_fail, indomain_random)]) + :: restart_geometric(1.1,5*experiments) + satisfy ; + +%% Trying to use a more deterministic search strategy for testing +% solve:: if (false) then seq_search([int_search(experiment_plate, first_fail, indomain_min)]++search_order) +% else (seq_search([int_search(experiment_plate, first_fail, indomain_random)]++search_order) +% :: restart_geometric(1.1,5*experiments)) endif +% satisfy ; + + +%solve::seq_search([int_search(experiment_plate, first_fail, indomain_random)]++ +%search_order) +%[int_search(experiment_plate_row, first_fail, indomain_random)]++[int_search(plates[i,..,..], random, indomain_max) | i in Plates]) +%[int_search(plates[i,..,..], random, indomain_max) | i in Plates ]) +% :: restart_geometric(1.1,1000) + % :: if testing==true then restart_none else restart_geometric(1.1,5*experiments) endif + %:: restart_linear(1000) + %:: restart_constant(10*experiments) + % satisfy ; + + + +%% Chuffed +% solve::seq_search([int_search(experiment_plate, random_order, indomain_min), int_search(experiment_row, random_order, indomain_min), int_search(plates, random_order, indomain_max)]) +% :: restart_geometric(1.5,100) +% satisfy; + + +%solve::seq_search([int_search(experiment_plate, random_order, indomain_min), priority_search(experiment_plate,[int_search([experiment_row[i]],input_order,indomain_split_random)| i in 1..experiments],smallest,complete)]) satisfy; + + +%solve::seq_search([int_search(experiment_plate, first_fail, indomain_random), int_search(experiment_row, random, indomain_split), int_search(plates, first_fail, indomain_random)]) satisfy; + +%solve satisfy; + + + +%solve::seq_search([int_search(experiment_plate, first_fail, indomain_random), int_search(experiment_row, first_fail, indomain_random), int_search(plates, random, indomain_max)]) satisfy; + + + + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% %%% Pretty printing %%% + +output [if testing \/ debugging then "\(numplates) plates\n" else "" endif]; + +output [if debugging /\ (emptywells>0) then "warning: there are \(emptywells) empty wells. Consider using them as controls.\n\n" else "" endif]; + +output [if debugging then "Spreading controls... \(spread_control) ... \(min(spread_control++[1])>0)\n" else "" endif]; + +output [if debugging then "Force spread controls... \(force_spread_controls)\n" else "" endif]; + +output [if debugging then "Random value... \(normal(0,1))\n" else "" endif]; + +output [if debugging then "Plate 1:\n" else "" endif]; + +%%TODO: change output so symbols are not hardcoded +output [if debugging then if fix(plates[i,j,k]) == 0 then " ." + elseif fix(plates[i,j,k]) == (experiments+1) then " +" %% Some control + elseif fix(plates[i,j,k]) == (experiments+2) then " -" %% Other control... + elseif fix(plates[i,j,k]) == (experiments+3) then " x" %% ... + elseif fix(plates[i,j,k]) == (experiments+4) then " ~" %% ... + elseif fix(plates[i,j,k]) == (experiments+5) then " *" %% other + elseif fix(plates[i,j,k]) == (experiments+6) then " #" %% other + elseif fix(plates[i,j,k]) == (experiments+7) then " &" %% other + elseif fix(plates[i,j,k]) == (experiments+8) then " @" %% other + elseif fix(plates[i,j,k]) == (experiments+9) then " <" %% other + elseif fix(plates[i,j,k]) == (experiments+10) then " >" %% other + elseif fix(plates[i,j,k]) == (experiments+11) then " $" %% other + elseif fix(plates[i,j,k]) == (experiments+12) then " X" %% other + elseif fix(plates[i,j,k]) > (experiments+12) then " Y" %% other + else " " endif ++ + %% Adding lines and headings %% + if j== num_rows_line /\ k == num_cols_line /\ i 0 then "plate_\(i)," ++ letters[fix(calculate_line_index(inner_empty_edge, size_empty_edge, num_rows_line, h, j))] ++ if fix(calculate_line_index(inner_empty_edge, size_empty_edge, num_cols_line, v, k)) < 10 then "0" else "" endif ++ "\(fix(calculate_line_index(inner_empty_edge, size_empty_edge, num_cols_line, v, k)))," ++ + +%%%% Compounds %%%% + +if fix(plates[i,j,k]) <= experiments then compound_names[(((fix(plates[i,j,k])-1) mod (compounds*max_compound_concentrations)) div max_compound_concentrations)+1] +++ "," ++ compound_concentration_names[(((fix(plates[i,j,k])-1) mod (compounds*max_compound_concentrations)) div max_compound_concentrations)+1,(((fix(plates[i,j,k])-1) mod (compounds*max_compound_concentrations)) mod max_compound_concentrations)+1] ++ "," ++ + +compound_names[(floor(((fix(plates[i,j,k])-1)/max_compound_concentrations)) mod compounds)+1] ++ "_" ++ compound_concentration_names[(floor(((fix(plates[i,j,k])-1)/max_compound_concentrations)) mod compounds)+1,((floor((fix(plates[i,j,k])-1))) mod max_compound_concentrations)+1] + + +%%%% Controls %%%% + +else +%control_names[((fix(plates[i,j,k])-experiments-1) mod num_controls)+1] ++ "," ++ +control_names[floor((fix(plates[i,j,k])-experiments-1)/max_control_concentrations)+1] ++ "," ++ +%% control concentration + +%control_concentration_names[((fix(plates[i,j,k])-experiments-1) mod num_controls)+1,floor((fix(plates[i,j,k])-experiments-1)/num_controls)+1] +control_concentration_names[floor((fix(plates[i,j,k])-experiments-1)/max_control_concentrations)+1,((fix(plates[i,j,k])-experiments-1) mod max_control_concentrations)+1] + +++ "," ++ +%% Latex name of the control. Right now it's the same as the name +control_names[floor((fix(plates[i,j,k])-experiments-1)/max_control_concentrations)+1] ++ "_" ++ control_concentration_names[floor((fix(plates[i,j,k])-experiments-1)/max_control_concentrations)+1,((fix(plates[i,j,k])-experiments-1) mod max_control_concentrations)+1] + +endif + +++ "\n" else "" endif| + i in Plates, j in Rows, k in Columns, v in 1..vertical_cell_lines, h in 1..horizontal_cell_lines]; + + + diff --git a/.agents/skills/plaid-platemap/scripts/plaid_platemap.py b/.agents/skills/plaid-platemap/scripts/plaid_platemap.py new file mode 100644 index 0000000..cf1795b --- /dev/null +++ b/.agents/skills/plaid-platemap/scripts/plaid_platemap.py @@ -0,0 +1,844 @@ +#!/usr/bin/env python3 +# /// script +# requires-python = ">=3.11" +# dependencies = ["openpyxl>=3.1,<4"] +# /// +"""Translate a hit-validation platemap workbook to/from a plaid model. + +The experiment workbook fixes a cell-plating layout with two irregular zones +(failing vs nonfailing wells). plaid (pharmbio/plaid) cannot express irregular +zones or 36-well DMSO density, so this tool models the compound-placement +problem, lets plaid solve it, relocates any out-of-zone placements, and +retranslates the result into a workbook in the original format. + +Subcommands: + extract workbook -> canonical layout JSON (zones, rosters, before-metrics) + dzn layout JSON -> plaid .dzn (compounds-only model) + solve run MiniZinc/Gecode per seed until a verified layout passes + apply write redesigned workbook (styles preserved) + QC sheet + verify independently re-check an output workbook against invariants + +Run with: uv run --script plaid_platemap.py --help +""" +from __future__ import annotations + +import argparse +import itertools +import json +import re +import subprocess +import sys +from collections import Counter +from dataclasses import dataclass, field +from datetime import date +from pathlib import Path + +import openpyxl +from openpyxl.cell import Cell +from openpyxl.styles import Font, PatternFill +from openpyxl.styles.colors import Color + +PLAID_COMMIT = "ecf4cfe6d8029dfadfceb79f5d6f2a1814d792c5" # pharmbio/plaid pin +GRID_SHEET = "Platemap grid" +VISUAL_SHEET = "Visual Platemaps" +QC_SHEET = "Platemap QC" +ROSTER_HEADER = "Randomized order" +TREATMENT_PREFIX = "treatment_plate " +ROW_LETTERS = "ABCDEFGH" +DESCRIPTION = "Translate platemap workbooks to/from plaid constraint models." + + +# ---------------------------------------------------------------- well utils + +def parse_well(well: str) -> tuple[str, int]: + m = re.fullmatch(r"([A-H])(\d{2})", str(well).strip()) + if not m: + raise ValueError(f"bad well label: {well!r}") + return m.group(1), int(m.group(2)) + + +def well_sort_key(well: str) -> tuple[str, int]: + return parse_well(well) + + +def well_label(row: str, col: int) -> str: + return f"{row}{col:02d}" + + +def king_distance(a: str, b: str) -> int: + (ra, ca), (rb, cb) = parse_well(a), parse_well(b) + return max(abs(ord(ra) - ord(rb)), abs(ca - cb)) + + +def set_cell(ws, row: int, col: int, value) -> None: + cell = ws.cell(row=row, column=col) + if not isinstance(cell, Cell): + raise RuntimeError(f"unexpected non-writable cell at {ws.title} r{row}c{col}") + cell.value = value + + +# ---------------------------------------------------------------- layout + +@dataclass +class Layout: + wells: list[dict] = field(default_factory=list) # pos/row/col/heart/cell_type/hf_type/treatments + plates: list[dict] = field(default_factory=list) # index/compounds/dose per plate + problems: list[str] = field(default_factory=list) + color_legend: dict | None = None # value -> fill spec, learned from the source workbook + + @property + def n_plates(self) -> int: + return len(self.plates) + + def by_pos(self) -> dict[str, dict]: + return {w["pos"]: w for w in self.wells} + + def failing(self) -> set[str]: + return {w["pos"] for w in self.wells if w["cell_type"] == "failing"} + + def nonfailing(self) -> set[str]: + return {w["pos"] for w in self.wells if w["cell_type"] != "failing"} + + def compound_wells(self, plate_idx: int) -> set[str]: + return {w["pos"] for w in self.wells if w["treatments"][plate_idx] != "DMSO"} + + +def load_layout(path: str) -> Layout: + payload = json.loads(Path(path).read_text()) + return Layout(wells=payload["wells"], plates=payload["plates"], + color_legend=payload.get("color_legend")) + + +def fill_spec(fill) -> dict | None: + """Serialize a cell fill (pattern + foreground color) for JSON storage.""" + if fill is None or fill.patternType is None: + return None + fg = fill.fgColor + spec: dict = {"patternType": fill.patternType, "fg": {"type": fg.type, "tint": fg.tint}} + if fg.type == "rgb" and isinstance(fg.rgb, str): + spec["fg"]["rgb"] = fg.rgb + if fg.type == "theme" and isinstance(fg.theme, int): + spec["fg"]["theme"] = fg.theme + if fg.type == "indexed" and isinstance(fg.indexed, int): + spec["fg"]["indexed"] = fg.indexed + return spec + + +def fill_from_spec(spec: dict | None) -> PatternFill: + """Rebuild a fill from a stored spec; empty fill for None.""" + if spec is None: + return PatternFill(fill_type=None) + fg = spec.get("fg", {}) + color = None + if fg.get("type") == "theme": + color = Color(theme=fg["theme"], tint=fg.get("tint", 0.0)) + elif fg.get("type") == "rgb": + color = Color(rgb=fg.get("rgb", "FF000000")) + if color is None: + return PatternFill(fill_type=None) + return PatternFill(patternType=spec.get("patternType", "solid"), fgColor=color) + + +def capture_color_legend(ws) -> tuple[dict, list[str]]: + """Learn the value -> fill legend from the source Visual Platemaps sheet. + + Keys: 'Failing', 'Nonfailing' (from ' ' sketch cells), + 'Cpd N' labels, and 'DMSO'. Cells without a fill do not define the legend + (plaid-era formula cells were uncolored); conflicting fills for the same + key are reported so validation can fail loudly. + """ + legend: dict[str, dict] = {} + conflicts: list[str] = [] + for row in ws.iter_rows(): + for c in row: + v = c.value + if not isinstance(v, str): + continue + s = v.strip() + key = None + m = re.fullmatch(r"(Failing|Nonfailing) \d+", s) + if m: + key = m.group(1) + elif re.fullmatch(r"Cpd \d+", s): + key = s + elif s == "DMSO": + key = "DMSO" + if key is None: + continue + spec = fill_spec(c.fill) + if spec is None: + continue + if key in legend and legend[key] != spec: + conflicts.append(f"{key}: {legend[key]} vs {spec}") + else: + legend[key] = spec + return legend, conflicts + + +def legend_fill(legend: dict | None, key: str | None) -> PatternFill: + """Fill for a legend key; empty (cleared) fill for None or unknown keys.""" + if legend and key and key in legend: + return fill_from_spec(legend[key]) + return PatternFill(fill_type=None) + + +def find_anchor(ws, label: str, col: int = 1) -> int | None: + """Row of the first cell in a column whose value matches label (case-insensitive).""" + for r in range(1, ws.max_row + 1): + v = ws.cell(row=r, column=col).value + if isinstance(v, str) and v.strip().lower() == label: + return r + return None + + +def read_grid_sheet(ws) -> tuple[list[dict], list[str], int]: + rows = list(ws.iter_rows(values_only=True)) + header = list(rows[0]) + n_plates = sum(1 for h in header if str(h).startswith(TREATMENT_PREFIX)) + expected = ["WellRow", "WellCol", "well_position", "heart_number", + "cell_type", "heart_failure_type"] + [f"{TREATMENT_PREFIX}{i}" for i in range(1, n_plates + 1)] + if header[: len(expected)] != expected: + raise ValueError(f"unexpected {GRID_SHEET} header: {header}") + wells, problems = [], [] + for r in rows[1:]: + if r[2] is None: + continue + pos = str(r[2]).strip() + row_letter, col = parse_well(pos) + if str(r[0]).strip() != row_letter or int(r[1]) != col: + problems.append(f"{pos}: WellRow/WellCol disagree with well_position") + wells.append({ + "pos": pos, "row": row_letter, "col": col, + "heart": int(r[3]), "cell_type": str(r[4]).strip(), + "hf_type": str(r[5]).strip(), + "treatments": [str(r[6 + i]).strip() for i in range(n_plates)], + }) + return wells, problems, n_plates + + +def read_rosters(ws, n_plates: int) -> tuple[list[dict], list[str]]: + """Read the 'Randomized order' table: New identifier / plate / dose columns.""" + header_row = None + for row in ws.iter_rows(values_only=True): + if any(str(v).strip() == ROSTER_HEADER for v in row if v is not None): + header_row = row + break + if header_row is None: + raise ValueError(f"'{ROSTER_HEADER}' table not found in {VISUAL_SHEET}") + col = {str(v).strip(): i for i, v in enumerate(header_row) if v is not None} + need = ["New identifier", "Validation plate map #", "Dose (uM)"] + missing = [k for k in need if k not in col] + if missing: + raise ValueError(f"roster table missing columns {missing}: {col}") + id_c, plate_c, dose_c = col["New identifier"], col["Validation plate map #"], col["Dose (uM)"] + plates: dict[int, dict] = {} + problems = [] + for row in ws.iter_rows(min_row=2, values_only=True): + label = row[plate_c] if plate_c < len(row) else None + if label is None or not str(label).strip().lower().startswith("plate"): + continue # not a roster row (skips headers, mini-table rows, blanks) + m = re.search(r"(\d+)", str(label)) + if not m: + problems.append(f"roster row {row[id_c]!r}: bad plate label {label!r}") + continue + plate_idx = int(m.group(1)) + ident = row[id_c] if id_c < len(row) else None + if not isinstance(ident, str) or not ident.strip(): + problems.append(f"plate {plate_idx} roster row has no identifier: {row[:6]!r}") + continue + plates.setdefault(plate_idx, {"compounds": [], "dose": None}) + plates[plate_idx]["compounds"].append(ident.strip()) + plates[plate_idx]["dose"] = row[dose_c] if dose_c < len(row) else None + ordered = [] + for i in range(1, n_plates + 1): + if i not in plates: + raise ValueError(f"no roster entries for plate {i}") + ordered.append({"index": i, "compounds": plates[i]["compounds"], "dose": plates[i]["dose"]}) + return ordered, problems + + +def layout_metrics(layout: Layout, plate_idx: int) -> dict: + """Per-plate cross-effect metrics for the assignments in layout.""" + by_cpd: dict[str, list[str]] = {} + for w in layout.wells: + t = w["treatments"][plate_idx] + if t != "DMSO": + by_cpd.setdefault(t, []).append(w["pos"]) + adjacent, min_dist, rows_used = [], {}, {} + for cpd, ws in by_cpd.items(): + pairs = [(king_distance(a, b), a, b) for a, b in itertools.combinations(sorted(ws), 2)] + adjacent += [(cpd, a, b) for d, a, b in pairs if d == 1] + min_dist[cpd] = min((d for d, _, _ in pairs), default=None) + rows_used[cpd] = "".join(sorted({parse_well(w)[0] for w in ws})) + counts = {cpd: len(ws) for cpd, ws in by_cpd.items()} + return { + "plate": plate_idx + 1, + "compound_counts": dict(sorted(counts.items())), + "same_compound_adjacent_pairs": len(adjacent), + "same_compound_adjacent_pairs_detail": sorted(adjacent), + "min_same_compound_distance": dict(sorted(min_dist.items())), + "rows_used_per_compound": dict(sorted(rows_used.items())), + } + + +def validate_layout(layout: Layout) -> None: + p = list(layout.problems) + positions = [w["pos"] for w in layout.wells] + if len(set(positions)) != len(positions): + p.append("duplicate well_position values") + rows = {w["row"] for w in layout.wells} + cols = {w["col"] for w in layout.wells} + combos = {(w["row"], w["col"]) for w in layout.wells} + if len(combos) != len(rows) * len(cols): + p.append(f"used region is not a full rectangle: {len(rows)} rows x {len(cols)} cols " + f"but {len(combos)} distinct wells; plaid models rectangular inner areas only") + failing = layout.failing() + for i in range(layout.n_plates): + cw = layout.compound_wells(i) + off_zone = sorted(cw - failing) + if off_zone: + p.append(f"plate {i + 1}: compounds on nonfailing wells {off_zone}") + counts = Counter(w["treatments"][i] for w in layout.wells if w["treatments"][i] != "DMSO") + roster = layout.plates[i]["compounds"] + if sorted(counts) != sorted(roster): + p.append(f"plate {i + 1}: grid compounds {sorted(counts)} != roster {sorted(roster)}") + if len(set(counts.values())) > 1: + p.append(f"plate {i + 1}: uneven replicate counts {dict(counts)}") + patterns = [layout.compound_wells(i) for i in range(layout.n_plates)] + if all(x == patterns[0] for x in patterns): + print(f"note: all {layout.n_plates} plates share one compound well pattern " + f"({len(patterns[0])} wells); the redesign keeps this property") + else: + p.append("plates use different compound well patterns; solve each plate separately") + if p: + raise SystemExit("layout validation failed:\n - " + "\n - ".join(p)) + + +def cmd_extract(args: argparse.Namespace) -> None: + wb = openpyxl.load_workbook(args.input, data_only=True) + if GRID_SHEET not in wb.sheetnames or VISUAL_SHEET not in wb.sheetnames: + raise SystemExit(f"workbook must contain '{GRID_SHEET}' and '{VISUAL_SHEET}' sheets") + wells, problems, n_plates = read_grid_sheet(wb[GRID_SHEET]) + rosters, roster_problems = read_rosters(wb[VISUAL_SHEET], n_plates) + color_legend, legend_conflicts = capture_color_legend(wb[VISUAL_SHEET]) + layout = Layout(wells=wells, plates=rosters, + problems=problems + roster_problems + legend_conflicts, color_legend=color_legend) + validate_layout(layout) + payload = { + "source": str(Path(args.input).resolve()), + "n_plates": n_plates, + "wells": layout.wells, + "plates": layout.plates, + "color_legend": layout.color_legend, + "before_metrics": [layout_metrics(layout, i) for i in range(n_plates)], + } + Path(args.output).write_text(json.dumps(payload, indent=2)) + m = payload["before_metrics"][0] + print(f"extracted {len(wells)} wells, {n_plates} plates -> {args.output}") + print(f"zone sizes: failing={len(layout.failing())} nonfailing={len(layout.nonfailing())}") + print(f"color legend keys: {sorted(color_legend)}") + print(f"plate 1 before: {m['same_compound_adjacent_pairs']} same-compound adjacent pairs, " + f"min distances {m['min_same_compound_distance']}") + + +# ---------------------------------------------------------------- dzn + +def cmd_dzn(args: argparse.Namespace) -> None: + layout = load_layout(args.layout) + # derive plate geometry from the used region and the requested empty edge + rows = sorted({w["row"] for w in layout.wells}) + cols = sorted({w["col"] for w in layout.wells}) + e = args.size_empty_edge + num_rows = ord(rows[-1]) - ord(rows[0]) + 1 + 2 * e + num_cols = cols[-1] - cols[0] + 1 + 2 * e + inner_rows = ord(rows[-1]) - ord(rows[0]) + 1 + inner_cols = cols[-1] - cols[0] + 1 + if inner_rows % 2 != 0 or inner_cols % 2 != 0: + raise SystemExit(f"plaid requires even inner dimensions; used region is " + f"{inner_rows} x {inner_cols}; choose another --size-empty-edge " + f"or adjust the layout region") + reps = Counter() + for i in range(layout.n_plates): + for w in layout.wells: + if w["treatments"][i] != "DMSO": + reps[w["treatments"][i]] += 1 + per_plate = [len(layout.plates[i]["compounds"]) for i in range(layout.n_plates)] + if len(set(per_plate)) != 1: + raise SystemExit(f"unequal compounds per plate {per_plate}; single-pattern dzn unsupported") + n_cpd = per_plate[0] + if len(set(reps.values())) != 1: + raise SystemExit(f"unequal replicate counts across compounds: {dict(reps)}") + seen: dict[str, set[int]] = {} + for i in range(layout.n_plates): + for c in layout.plates[i]["compounds"]: + seen.setdefault(c, set()).add(i) + multi = {c: sorted(ps) for c, ps in seen.items() if len(ps) > 1} + if multi: + raise SystemExit(f"compounds appearing on multiple plates are unsupported by the " + f"single-pattern model: {multi}") + reps_per = next(iter(reps.values())) + dose = layout.plates[0]["dose"] + conc = f"{dose} uM" if dose is not None else "1x" + names = ", ".join(f'"Cpd {i}"' for i in range(1, n_cpd + 1)) + conc_names = " ".join(f'| "{conc}"' for _ in range(n_cpd)) + text = f"""%% Generated by plaid_platemap.py dzn - compound placement on the used plate area. +%% Physical plate derived from the used region {rows[0]}{cols[0]:02d}-{rows[-1]}{cols[-1]:02d} +%% plus a {e}-well empty edge on each side ({num_rows} x {num_cols} total). +%% Compounds only: the DMSO controls are too dense for plaid's same-kind +%% separation and are assigned as the complement of the compound wells by this +%% tool. +num_rows = {num_rows}; +num_cols = {num_cols}; + +horizontal_cell_lines = 1; +vertical_cell_lines = 1; +allow_empty_wells = true; +size_empty_edge = {args.size_empty_edge}; + +concentrations_on_different_rows = true; +concentrations_on_different_columns = true; +replicates_on_different_plates = false; +replicates_on_same_plate = true; + +%%% Compounds %%% +compounds = {n_cpd}; +compound_names = [{names}]; +compound_replicates = [{", ".join(str(reps_per) for _ in range(n_cpd))}]; +compound_concentrations = [{", ".join("1" for _ in range(n_cpd))}]; +compound_concentration_names = [{conc_names} |]; +compound_concentration_indicators = [""]; + +%%% Combinations %%% +combinations = 0; +combination_concentrations = 0; +combination_names = []; +combination_concentration_names = []; + +%%% Controls %%% +num_controls = 0; +control_concentrations = []; +control_replicates = []; +control_names = []; +control_concentration_names = []; +""" + Path(args.output).write_text(text) + print(f"wrote {args.output}: {n_cpd} compounds x {reps_per} replicates, concentration {conc!r}") + + +# ---------------------------------------------------------------- solve + +def run_plaid(mzn: Path, dzn: Path, seed: int, timeout: int) -> dict[str, str]: + cmd = ["minizinc", "--solver", "gecode", "-r", str(seed), + "--time-limit", str(timeout * 1000), str(mzn), str(dzn)] + try: + out = subprocess.run(cmd, capture_output=True, text=True, + timeout=timeout + 5, check=True) + except subprocess.TimeoutExpired: + raise RuntimeError(f"seed {seed}: minizinc timed out after {timeout + 5}s") + except subprocess.CalledProcessError as e: + raise RuntimeError(f"seed {seed}: minizinc failed: {e.stderr[-500:]}") + assign = {} + for line in out.stdout.splitlines(): + m = re.match(r"^plate_1,([A-H])(\d{2}),(.+?),", line.strip()) + if m: + assign[well_label(m.group(1), int(m.group(2)))] = m.group(3) + if not assign: + raise RuntimeError(f"seed {seed}: no layout rows in minizinc output") + return assign + + +def relocate(assign: dict[str, str], layout: Layout, min_dist: int) -> dict[str, str]: + """Move compound wells that fall on out-of-zone wells onto free in-zone wells.""" + failing = layout.failing() + placed = dict(assign) + displaced = sorted(w for w in placed if w not in failing) + free = sorted(failing - set(placed)) + for w in displaced: + cpd = placed.pop(w) + same = [x for x, c in placed.items() if c == cpd] + cands = [f for f in free if all(king_distance(f, s) >= min_dist for s in same)] + if not cands: + raise RuntimeError(f"no in-zone relocation candidate for {w} ({cpd})") + + def score(f: str) -> tuple: + row_count = Counter(parse_well(x)[0] for x in placed) + return (0 if parse_well(f)[0] == parse_well(w)[0] else 1, + -min(king_distance(f, s) for s in same), + row_count.get(parse_well(f)[0], 0), + well_sort_key(f)) + + best = min(cands, key=score) + placed[best] = cpd + free.remove(best) + return placed + + +def check_pattern(assign: dict[str, str], layout: Layout, n_cpd: int, reps: int, + min_dist: int) -> tuple[list[str], dict]: + problems, by_cpd = [], {} + for w, c in assign.items(): + by_cpd.setdefault(c, []).append(w) + if sorted(by_cpd) != [f"Cpd {i}" for i in range(1, n_cpd + 1)]: + problems.append(f"compound labels {sorted(by_cpd)} != expected {n_cpd} compounds") + for c, ws in sorted(by_cpd.items()): + if len(ws) != reps: + problems.append(f"{c}: {len(ws)} wells, expected {reps}") + off_zone = sorted(set(assign) - layout.failing()) + if off_zone: + problems.append(f"wells outside failing zone: {off_zone}") + for c, ws in sorted(by_cpd.items()): + for a, b in itertools.combinations(sorted(ws), 2): + if king_distance(a, b) < min_dist: + problems.append(f"{c}: {a}/{b} closer than {min_dist} (king)") + metrics = { + "min_same_compound_distance": { + c: min(king_distance(a, b) for a, b in itertools.combinations(sorted(ws), 2)) + for c, ws in sorted(by_cpd.items())}, + "rows_used_per_compound": { + c: "".join(sorted({parse_well(w)[0] for w in ws})) for c, ws in sorted(by_cpd.items())}, + "compound_wells_per_row": dict(sorted(Counter(parse_well(w)[0] for w in assign).items())), + } + return problems, metrics + + +def cmd_solve(args: argparse.Namespace) -> None: + layout = load_layout(args.layout) + n_cpd = len(layout.plates[0]["compounds"]) + reps = len(layout.compound_wells(0)) // n_cpd + seeds = list(range(args.first_seed, args.first_seed + args.max_seeds)) + failures = [] + for seed in seeds: + try: + assign = run_plaid(Path(args.mzn), Path(args.dzn), seed, args.timeout) + placed = relocate(assign, layout, args.min_dist) + except RuntimeError as e: + failures.append(f"seed {seed}: {e}") + continue + problems, metrics = check_pattern(placed, layout, n_cpd, reps, args.min_dist) + if problems: + failures.append(f"seed {seed}: " + "; ".join(problems[:4])) + continue + n_displaced = sum(1 for w in assign if w not in layout.failing()) + solution = { + "seed": seed, + "min_dist": args.min_dist, + "displaced_wells_relocated": n_displaced, + "plaid_commit": PLAID_COMMIT, + "mzn": str(Path(args.mzn).resolve()), + "dzn": str(Path(args.dzn).resolve()), + "n_compounds": n_cpd, + "replicates": reps, + "pattern": placed, + "metrics": metrics, + } + Path(args.output).write_text(json.dumps(solution, indent=2)) + print(f"seed {seed}: PASS (relocated {n_displaced} out-of-zone wells)") + print(f" min same-compound distance: {metrics['min_same_compound_distance']}") + print(f" compound wells per row: {metrics['compound_wells_per_row']}") + print(f"solution -> {args.output}") + return + raise SystemExit(f"no seed passed out of {len(seeds)} tried:\n - " + "\n - ".join(failures[:10])) + + +# ---------------------------------------------------------------- apply + +def pattern_treatments(layout: Layout, pattern: dict[str, str]) -> list[dict]: + """Materialize per-well treatments for every plate from a well->'Cpd N' pattern.""" + new_wells = [] + for w in layout.wells: + cpd = pattern.get(w["pos"]) + treatments = ["DMSO" if cpd is None + else layout.plates[p]["compounds"][int(cpd.split()[1]) - 1] + for p in range(layout.n_plates)] + new_wells.append({**w, "treatments": treatments}) + return new_wells + + +def cmd_apply(args: argparse.Namespace) -> None: + src, dst = Path(args.input).resolve(), Path(args.output).resolve() + if src == dst: + raise SystemExit("refusing to overwrite the input workbook; choose a different --output") + layout = load_layout(args.layout) + solution = json.loads(Path(args.solution).read_text()) + pattern: dict[str, str] = solution["pattern"] + legend = layout.color_legend + wb = openpyxl.load_workbook(args.input) + + def set_value_fill(ws, row: int, col: int, value, fill_key: str | None) -> None: + set_cell(ws, row, col, value) + cell = ws.cell(row=row, column=col) + if not isinstance(cell, Cell): + raise RuntimeError(f"unexpected non-writable cell at {ws.title} r{row}c{col}") + cell.fill = legend_fill(legend, fill_key) + + # 1. Platemap grid: new per-plate treatments (grid carries no fill scheme + # in V3, so no recoloring here; the Visual Platemaps sheet holds the colors) + ws = wb[GRID_SHEET] + for row_cells in ws.iter_rows(min_row=2, min_col=1, max_col=6 + layout.n_plates): + pos = row_cells[2].value + if pos is None: + continue + pos = str(pos).strip() + for plate_idx in range(layout.n_plates): + cpd = pattern.get(pos) + value = "DMSO" if cpd is None else layout.plates[plate_idx]["compounds"][int(cpd.split()[1]) - 1] + set_cell(ws, row_cells[0].row or 0, 7 + plate_idx, value) + + # 2. Visual Platemaps: hearts + compounds sketches now mirror the grid, + # with fills recolored from the learned legend so flipped wells match + ws = wb[VISUAL_SHEET] + by_pos = layout.by_pos() + rows_in_region = sorted({w["row"] for w in layout.wells}) + cols_in_region = sorted({w["col"] for w in layout.wells}) + + hearts_anchor = find_anchor(ws, "hearts") + if hearts_anchor is not None: + for row_letter in rows_in_region: + r = hearts_anchor + 1 + (ord(row_letter) - ord("A")) + for c in cols_in_region: + w = by_pos[well_label(row_letter, c)] + if w["cell_type"] == "failing": + label, key = f"Failing {w['heart']}", "Failing" + else: + label, key = f"Nonfailing {w['heart']}", "Nonfailing" + set_value_fill(ws, r, c + 1, label, key) + compounds_anchor = find_anchor(ws, "compounds") + if compounds_anchor is not None: + for row_letter in rows_in_region: + r = compounds_anchor + 1 + (ord(row_letter) - ord("A")) + for c in cols_in_region: + cpd = pattern.get(well_label(row_letter, c)) + set_value_fill(ws, r, c + 1, cpd if cpd is not None else "DMSO", cpd or "DMSO") + # per-well block: located by the 'Initial list' header in column C; rows + # below it list well index 1..N in raster order. C holds the "Initial list" + # (sequential grouping: reps wells per compound in raster order), D the + # "Randomized" assignment from the new pattern, E.. per-plate K-IDs as + # values (V3 kept formulas here; values survive programmatic re-reads). + # Fix the stale plate headers (V3 mislabels plates 4/5) when present. + positions = sorted(by_pos, key=well_sort_key) + n_cpd = len(layout.plates[0]["compounds"]) + reps = sum(1 for pos in positions if pos in pattern) // n_cpd if n_cpd else 0 + block_header = None + for r in range(1, ws.max_row + 1): + v = ws.cell(row=r, column=3).value + if isinstance(v, str) and v.strip().lower() == "initial list": + block_header = r + break + if block_header is not None: + for i in range(1, layout.n_plates + 1): + set_cell(ws, block_header, 4 + i, f"Platemap {i}") + initial: dict[str, str] = {} + for k, pos in enumerate([p for p in positions if p in pattern], start=1): + initial[pos] = f"Cpd {(k - 1) // reps + 1}" + for idx, pos in enumerate(positions, start=1): + r = block_header + idx + cpd = pattern.get(pos) + set_cell(ws, r, 2, idx) + set_value_fill(ws, r, 3, initial.get(pos), initial.get(pos)) + set_value_fill(ws, r, 4, cpd, cpd) + for plate_idx in range(layout.n_plates): + value = "DMSO" if cpd is None else layout.plates[plate_idx]["compounds"][int(cpd.split()[1]) - 1] + set_value_fill(ws, r, 5 + plate_idx, value, "DMSO" if value == "DMSO" else None) + + # 3. QC sheet with before/after metrics + if QC_SHEET in wb.sheetnames: + del wb[QC_SHEET] + qc = wb.create_sheet(QC_SHEET) + new_layout = Layout(wells=pattern_treatments(layout, pattern), plates=layout.plates) + before = layout_metrics(layout, 0) + after = layout_metrics(new_layout, 0) + qc_rows = [ + ["Platemap QC - plaid redesign"], + [f"Generated {date.today().isoformat()} from {src.name}"], + [f"Model: pharmbio/plaid plate-design.mzn @ {solution.get('plaid_commit', PLAID_COMMIT)[:7]}, " + f"solver Gecode via MiniZinc, seed {solution['seed']}"], + [f"Compound wells relocated from out-of-zone plaid positions: {solution['displaced_wells_relocated']}"], + [f"Metrics identical for all {layout.n_plates} plates: one shared well pattern, relabeled per plate."], + [], + ["Metric", "Before", "After"], + ["Same-compound adjacent pairs (king distance 1), per plate", + before["same_compound_adjacent_pairs"], after["same_compound_adjacent_pairs"]], + ["Minimum same-compound distance (king), plate 1", + min(before["min_same_compound_distance"].values()), + min(after["min_same_compound_distance"].values())], + ["Rows used per compound, plate 1", + str(before["rows_used_per_compound"]), str(after["rows_used_per_compound"])], + ["Compound wells per row", + str(dict(sorted(Counter(parse_well(w)[0] for w in pattern).items()))), ""], + [], + ["Design notes"], + ["Compounds are spread across the full failing-heart region; DMSO fills the remaining " + "failing wells and all nonfailing wells (36 DMSO per plate, unchanged from V3)."], + ["Same-compound wells are at least 2 king-moves apart (one well between in every direction), " + "enforced by plaid's 2x2-block constraint plus post-solve relocation and verification."], + ["Corrected stale 'Visual Platemaps' content to match 'Platemap grid': plate headers 4/5 " + "and the hearts/compounds sketches (the V3 sketches were stylized and disagreed with the grid)."], + ["Per-well block formulas (E41:I64) were flattened to values matching the new layout."], + ] + for r, row in enumerate(qc_rows, start=1): + for c, v in enumerate(row, start=1): + cell = qc.cell(row=r, column=c, value=v) + if r in (1, 7, 11): + cell.font = Font(bold=True) + qc.column_dimensions["A"].width = 72 + wb.save(dst) + print(f"wrote {dst}") + print(f"QC: adjacent pairs {before['same_compound_adjacent_pairs']} -> " + f"{after['same_compound_adjacent_pairs']}; min distance " + f"{min(before['min_same_compound_distance'].values())} -> " + f"{min(after['min_same_compound_distance'].values())}") + + +# ---------------------------------------------------------------- verify + +def check_colors(ws, check: Layout, layout: Layout) -> list[str]: + """Verify the Visual Platemaps sketch fills match the learned legend. + + Checks the hearts sketch and the compounds sketch against the zones and + the plate-1 pattern. Returns a list of mismatch descriptions. + """ + problems = [] + legend = layout.color_legend or {} + by_pos = check.by_pos() + rows_in_region = sorted({w["row"] for w in check.wells}) + cols_in_region = sorted({w["col"] for w in check.wells}) + + hearts_anchor = find_anchor(ws, "hearts") + if hearts_anchor is not None and legend: + for row_letter in rows_in_region: + r = hearts_anchor + 1 + (ord(row_letter) - ord("A")) + for c in cols_in_region: + pos = well_label(row_letter, c) + expected_zone = "Failing" if by_pos[pos]["cell_type"] == "failing" else "Nonfailing" + cell = ws.cell(row=r, column=c + 1) + if fill_spec(cell.fill) != legend.get(expected_zone): + problems.append(f"hearts sketch {pos} ({cell.coordinate}): fill does not " + f"match {expected_zone} legend") + compounds_anchor = find_anchor(ws, "compounds") + if compounds_anchor is not None and legend: + # map well -> Cpd N using the plate-1 treatments of the checked layout + cpd_of_pos = {} + for w in check.wells: + t = w["treatments"][0] + if t != "DMSO": + try: + cpd_of_pos[w["pos"]] = f"Cpd {layout.plates[0]['compounds'].index(t) + 1}" + except ValueError: + pass + for row_letter in rows_in_region: + r = compounds_anchor + 1 + (ord(row_letter) - ord("A")) + for c in cols_in_region: + pos = well_label(row_letter, c) + expected_key = cpd_of_pos.get(pos, "DMSO") + cell = ws.cell(row=r, column=c + 1) + if fill_spec(cell.fill) != legend.get(expected_key): + problems.append(f"compounds sketch {pos} ({cell.coordinate}): fill does not " + f"match {expected_key} legend") + return problems + + +def cmd_verify(args: argparse.Namespace) -> None: + layout = load_layout(args.layout) + n_cpd = len(layout.plates[0]["compounds"]) + reps = len(layout.compound_wells(0)) // n_cpd + wb = openpyxl.load_workbook(args.input, data_only=True) + if GRID_SHEET not in wb.sheetnames: + raise SystemExit(f"no {GRID_SHEET} sheet in {args.input}") + wells, problems, n_plates = read_grid_sheet(wb[GRID_SHEET]) + check = Layout(wells=wells, plates=layout.plates, problems=problems) + failures = [] + if {w["pos"] for w in check.wells} != {w["pos"] for w in layout.wells}: + failures.append("well set differs from extracted layout") + if check.failing() != layout.failing(): + failures.append("failing zone changed") + patterns = [] + for i in range(n_plates): + m = layout_metrics(check, i) + print(f"plate {m['plate']}: adjacent_pairs={m['same_compound_adjacent_pairs']} " + f"min_dist={m['min_same_compound_distance']} counts={m['compound_counts']}") + counts = m["compound_counts"] + if sorted(counts) != sorted(layout.plates[i]["compounds"]): + failures.append(f"plate {i + 1}: compounds {sorted(counts)} != roster {sorted(layout.plates[i]['compounds'])}") + if counts and any(v != reps for v in counts.values()): + failures.append(f"plate {i + 1}: replicate counts {counts} != {reps}") + off = sorted(check.compound_wells(i) - check.failing()) + if off: + failures.append(f"plate {i + 1}: compounds on nonfailing wells {off}") + if m["same_compound_adjacent_pairs"]: + failures.append(f"plate {i + 1}: {m['same_compound_adjacent_pairs']} same-compound adjacent pairs") + if m["min_same_compound_distance"] and min(m["min_same_compound_distance"].values()) < args.min_dist: + failures.append(f"plate {i + 1}: min same-compound distance below {args.min_dist}") + patterns.append(check.compound_wells(i)) + if not all(p == patterns[0] for p in patterns): + failures.append("plates do not share one compound well pattern") + # color coding: every sketch cell's fill must match the learned legend + if layout.color_legend: + wsv = openpyxl.load_workbook(args.input)[VISUAL_SHEET] + color_fail = check_colors(wsv, check, layout) + failures += color_fail + print(f"color check: {len(color_fail)} mismatches") + if args.reference: + ref = load_layout(args.reference) + for i in range(n_plates): + b = layout_metrics(ref, i) + print(f"plate {i + 1} before: adjacent_pairs={b['same_compound_adjacent_pairs']} " + f"min_dist={b['min_same_compound_distance']}") + if args.output: + Path(args.output).write_text(json.dumps({"failures": failures}, indent=2)) + if failures: + raise SystemExit("VERIFY FAILED:\n - " + "\n - ".join(failures)) + print(f"VERIFY OK: {n_plates} plates, all invariants hold " + f"(min same-compound distance >= {args.min_dist})") + + +# ---------------------------------------------------------------- main + +def main(argv: list[str] | None = None) -> None: + parser = argparse.ArgumentParser(description=DESCRIPTION) + sub = parser.add_subparsers(dest="command", required=True) + + p = sub.add_parser("extract", help="read workbook -> layout JSON") + p.add_argument("--input", required=True) + p.add_argument("--output", required=True) + p.set_defaults(func=cmd_extract) + + p = sub.add_parser("dzn", help="layout JSON -> plaid .dzn (compounds-only model)") + p.add_argument("--layout", required=True) + p.add_argument("--output", required=True) + p.add_argument("--size-empty-edge", type=int, default=1) + p.set_defaults(func=cmd_dzn) + + default_mzn = Path(__file__).resolve().parent.parent / "assets" / "plate-design.mzn" + p = sub.add_parser("solve", help="run MiniZinc per seed until a verified layout passes") + p.add_argument("--mzn", default=str(default_mzn)) + p.add_argument("--dzn", required=True) + p.add_argument("--layout", required=True) + p.add_argument("--output", required=True) + p.add_argument("--min-dist", type=int, default=2, + help="minimum king distance between same-compound wells (default 2)") + p.add_argument("--first-seed", type=int, default=1) + p.add_argument("--max-seeds", type=int, default=25) + p.add_argument("--timeout", type=int, default=60, help="per-seed minizinc timeout in seconds") + p.set_defaults(func=cmd_solve) + + p = sub.add_parser("apply", help="write redesigned workbook from a solution") + p.add_argument("--input", required=True, help="original workbook (kept untouched)") + p.add_argument("--output", required=True) + p.add_argument("--layout", required=True) + p.add_argument("--solution", required=True) + p.set_defaults(func=cmd_apply) + + p = sub.add_parser("verify", help="independently re-check an output workbook") + p.add_argument("--input", required=True) + p.add_argument("--layout", required=True) + p.add_argument("--reference", help="layout JSON of the original, for a before/after report") + p.add_argument("--output", help="optional JSON report path") + p.add_argument("--min-dist", type=int, default=2) + p.set_defaults(func=cmd_verify) + + args = parser.parse_args(argv) + args.func(args) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/.agents/skills/plaid-platemap/tests/test_plaid_platemap.py b/.agents/skills/plaid-platemap/tests/test_plaid_platemap.py new file mode 100644 index 0000000..129ce18 --- /dev/null +++ b/.agents/skills/plaid-platemap/tests/test_plaid_platemap.py @@ -0,0 +1,302 @@ +# /// script +# requires-python = ">=3.11" +# dependencies = ["openpyxl>=3.1,<4", "pytest>=8,<9"] +# /// +"""Round-trip tests for plaid_platemap.py on a synthetic mini workbook. + +Builds an in-memory workbook with the V3 structure (grid + visual sheets, two +zones, two plates) and exercises extract/dzn/apply/verify and the pure-Python +helpers (relocate/check_pattern). The minizinc solve step is covered by an +end-to-end test that skips when minizinc or the vendored assets are missing. +""" +from __future__ import annotations + +import importlib.util +import json +import shutil +import subprocess +import sys +from pathlib import Path + +import openpyxl +import pytest + +SKILL = Path(__file__).resolve().parents[1] +SCRIPT = SKILL / "scripts" / "plaid_platemap.py" +ASSETS = SKILL / "assets" + +spec = importlib.util.spec_from_file_location("plaid_platemap", SCRIPT) +assert spec is not None and spec.loader is not None +pp = importlib.util.module_from_spec(spec) +sys.modules[spec.name] = pp # dataclass introspection requires registration +spec.loader.exec_module(pp) + +# Synthetic experiment: rows B-G x cols 2-5, 2 compounds x 2 reps per plate, +# 2 plates. Failing zone: rows D,E fully plus B02,B03 (mirrors the two-zone +# structure of the real workbook at a small scale). +ROWS = "BCDEFG" +COLS = [2, 3, 4, 5] +FAILING = {f"D{c:02d}" for c in COLS} | {f"E{c:02d}" for c in COLS} | {"B02", "B03"} +ROSTER = {1: ["C-100", "C-200"], 2: ["C-300", "C-400"]} + + +def build_workbook(path: Path) -> None: + wb = openpyxl.Workbook() + grid = wb.active + grid.title = pp.GRID_SHEET + grid.append(["WellRow", "WellCol", "well_position", "heart_number", "cell_type", + "heart_failure_type", "treatment_plate 1", "treatment_plate 2"]) + # V3-like bad pattern: same-compound wells adjacent, packed into rows D-E + bad = {"D02": 0, "D03": 0, "E04": 1, "E05": 1} + for row_letter in ROWS: + for col in COLS: + pos = f"{row_letter}{col:02d}" + failing = pos in FAILING + t = bad.get(pos) if failing else None + grid.append([row_letter, col, pos, 25 if failing else 7, + "failing" if failing else "nonfailing", + "dilated_cardiomyopathy" if failing else "NA", + "DMSO" if t is None else ROSTER[1][t], + "DMSO" if t is None else ROSTER[2][t]]) + vis = wb.create_sheet(pp.VISUAL_SHEET) + # roster table at row 22, cols N-S (14-19), same shape as the real workbook + vis.cell(row=22, column=14, value="Randomized order") + vis.cell(row=22, column=15, value="UCD") + vis.cell(row=22, column=16, value="New identifier") + vis.cell(row=22, column=17, value="Validation plate map #") + vis.cell(row=22, column=19, value="Dose (uM)") + order = 0 + for plate, cpds in sorted(ROSTER.items()): + for c in cpds: + order += 1 + vis.cell(row=22 + order, column=14, value=order) + vis.cell(row=22 + order, column=16, value=c) + vis.cell(row=22 + order, column=17, value=f"Plate {plate}") + vis.cell(row=22 + order, column=19, value=10) + # hearts + compounds sketches with a V3-like color legend + from openpyxl.styles import PatternFill + fills = { + "Failing": PatternFill("solid", fgColor="FFDD8368"), + "Nonfailing": PatternFill("solid", fgColor="FF92D050"), + "Cpd 1": PatternFill("solid", fgColor="FFFFC000"), + "Cpd 2": PatternFill("solid", fgColor="FF00B0F0"), + "DMSO": PatternFill("solid", fgColor="FFBFBFBF"), + } + vis.cell(row=1, column=1, value="Hearts") + vis.cell(row=11, column=1, value="Compounds") + vis.cell(row=22, column=3, value="Initial list") + vis.cell(row=22, column=4, value="Randomized") + for col in COLS: + vis.cell(row=1, column=col + 1, value=col) + vis.cell(row=11, column=col + 1, value=col) + for row_letter in ROWS: + r = 1 + ord(row_letter) - ord("A") + 1 + for col in COLS: + pos = f"{row_letter}{col:02d}" + failing = pos in FAILING + cell = vis.cell(row=r, column=col + 1, + value=f"Failing 25" if failing else "Nonfailing 7") + cell.fill = fills["Failing" if failing else "Nonfailing"] + r2 = 11 + ord(row_letter) - ord("A") + 1 + for col in COLS: + pos = f"{row_letter}{col:02d}" + t = bad.get(pos) if pos in FAILING else None + label = f"Cpd {t + 1}" if t is not None else "DMSO" + cell = vis.cell(row=r2, column=col + 1, value=label) + cell.fill = fills[label] + # per-well block rows 23+: raster well index 1..N with C/D labels and E.. plate values + positions = sorted([f"{r}{c:02d}" for r in ROWS for c in COLS], + key=lambda w: (w[0], int(w[1:]))) + for idx, pos in enumerate(positions, start=1): + r = 22 + idx + vis.cell(row=r, column=2, value=idx) + t = bad.get(pos) if pos in FAILING else None + if t is not None: + vis.cell(row=r, column=3, value=f"Cpd {t + 1}") + vis.cell(row=r, column=4, value=f"Cpd {t + 1}") + for plate, cpds in sorted(ROSTER.items()): + v = "DMSO" if t is None else cpds[t] + c = vis.cell(row=r, column=4 + plate, value=v) + if v == "DMSO": + c.fill = fills["DMSO"] + wb.save(path) + + +@pytest.fixture() +def workbook(tmp_path: Path) -> Path: + path = tmp_path / "mini.xlsx" + build_workbook(path) + return path + + +def run_cli(*argv: str, expect_ok: bool = True) -> subprocess.CompletedProcess: + proc = subprocess.run([sys.executable, str(SCRIPT), *argv], + capture_output=True, text=True, timeout=300) + if expect_ok and proc.returncode != 0: + raise AssertionError(f"command failed: {argv}\nstdout:\n{proc.stdout}\nstderr:\n{proc.stderr}") + return proc + + +def test_well_helpers() -> None: + assert pp.parse_well("C10") == ("C", 10) + assert pp.well_label("C", 10) == "C10" + assert pp.king_distance("D02", "E03") == 1 + assert pp.king_distance("D02", "E04") == 2 + assert pp.well_sort_key("B11") < pp.well_sort_key("C02") + with pytest.raises(ValueError): + pp.parse_well("C1") + with pytest.raises(ValueError): + pp.parse_well("I01") + + +def test_extract_reports_v3_style_problem(workbook: Path, tmp_path: Path) -> None: + layout_path = tmp_path / "layout.json" + run_cli("extract", "--input", str(workbook), "--output", str(layout_path)) + payload = json.loads(layout_path.read_text()) + assert payload["n_plates"] == 2 + assert len(payload["wells"]) == len(ROWS) * len(COLS) + before = payload["before_metrics"][0] + assert before["same_compound_adjacent_pairs"] > 0 + assert min(before["min_same_compound_distance"].values()) == 1 + assert [p["compounds"] for p in payload["plates"]] == [ROSTER[1], ROSTER[2]] + + +def test_dzn_matches_roster(workbook: Path, tmp_path: Path) -> None: + layout_path = tmp_path / "layout.json" + dzn_path = tmp_path / "plate.dzn" + run_cli("extract", "--input", str(workbook), "--output", str(layout_path)) + run_cli("dzn", "--layout", str(layout_path), "--output", str(dzn_path)) + text = dzn_path.read_text() + assert "compounds = 2;" in text + assert "compound_replicates = [2, 2];" in text + assert '"10 uM"' in text + # mini fixture region B-G x 2-5 with edge 1 -> 8 x 6 plate + assert "num_rows = 8;" in text + assert "num_cols = 6;" in text + assert "size_empty_edge = 1;" in text + + +def make_layout(): + return pp.Layout( + wells=[{"pos": f"{r}{c:02d}", "row": r, "col": c, + "heart": 25 if f"{r}{c:02d}" in FAILING else 7, + "cell_type": "failing" if f"{r}{c:02d}" in FAILING else "nonfailing", + "hf_type": "NA", "treatments": ["DMSO", "DMSO"]} + for r in ROWS for c in COLS], + plates=[{"index": 1, "compounds": ROSTER[1], "dose": 10}, + {"index": 2, "compounds": ROSTER[2], "dose": 10}], + ) + + +def test_relocate_and_check_pattern() -> None: + layout = make_layout() + # plaid-style assignment: two wells land on nonfailing wells and must be relocated + assign = {"D02": "Cpd 1", "C04": "Cpd 2", "B11": "Cpd 1", "G04": "Cpd 2"} + placed = pp.relocate(assign, layout, min_dist=2) + assert "B11" not in placed and "G04" not in placed + assert set(placed) <= layout.failing() + problems, metrics = pp.check_pattern(placed, layout, n_cpd=2, reps=2, min_dist=2) + assert problems == [] + assert min(metrics["min_same_compound_distance"].values()) >= 2 + + +def test_check_pattern_flags_violations() -> None: + layout = make_layout() + bad = {"D02": "Cpd 1", "D03": "Cpd 1", "E02": "Cpd 2", "E03": "Cpd 2"} + problems, _ = pp.check_pattern(bad, layout, n_cpd=2, reps=2, min_dist=2) + assert any("closer than 2" in p for p in problems) + + +def test_apply_and_verify_roundtrip(workbook: Path, tmp_path: Path) -> None: + layout_path = tmp_path / "layout.json" + solution_path = tmp_path / "solution.json" + out_path = tmp_path / "mini-plaid.xlsx" + run_cli("extract", "--input", str(workbook), "--output", str(layout_path)) + solution = { + "seed": 1, "min_dist": 2, "displaced_wells_relocated": 0, + "plaid_commit": pp.PLAID_COMMIT, "n_compounds": 2, "replicates": 2, + "pattern": {"B02": "Cpd 1", "D04": "Cpd 1", "E02": "Cpd 2", "E05": "Cpd 2"}, + "metrics": {}, + } + solution_path.write_text(json.dumps(solution)) + run_cli("apply", "--input", str(workbook), "--output", str(out_path), + "--layout", str(layout_path), "--solution", str(solution_path)) + proc = run_cli("verify", "--input", str(out_path), "--layout", str(layout_path), + "--reference", str(layout_path)) + assert "VERIFY OK" in proc.stdout + # grid treatments relabeled per plate from the same pattern + wb = openpyxl.load_workbook(out_path, data_only=True) + ws = wb[pp.GRID_SHEET] + rows = {r[2]: r for r in ws.iter_rows(values_only=True) if r[2]} + assert rows["B02"][6] == "C-100" and rows["B02"][7] == "C-300" + assert rows["E02"][6] == "C-200" and rows["E02"][7] == "C-400" + assert rows["D04"][6] == "C-100" and rows["E05"][7] == "C-400" + assert rows["B04"][6] == "DMSO" # nonfailing stays DMSO + # QC sheet exists with before/after + assert pp.QC_SHEET in wb.sheetnames + # original untouched + proc = run_cli("verify", "--input", str(workbook), "--layout", str(layout_path), + expect_ok=False) + assert proc.returncode != 0 # original still has adjacent pairs + + +def test_verify_catches_wrong_colors(workbook: Path, tmp_path: Path) -> None: + layout_path = tmp_path / "layout.json" + solution_path = tmp_path / "solution.json" + out_path = tmp_path / "mini-plaid.xlsx" + run_cli("extract", "--input", str(workbook), "--output", str(layout_path)) + solution = { + "seed": 1, "min_dist": 2, "displaced_wells_relocated": 0, + "plaid_commit": pp.PLAID_COMMIT, "n_compounds": 2, "replicates": 2, + "pattern": {"B02": "Cpd 1", "D04": "Cpd 1", "E02": "Cpd 2", "E05": "Cpd 2"}, + "metrics": {}, + } + solution_path.write_text(json.dumps(solution)) + run_cli("apply", "--input", str(workbook), "--output", str(out_path), + "--layout", str(layout_path), "--solution", str(solution_path)) + # find a flipped cell: a sketch cell holding Cpd 2 whose old fill was DMSO + wb = openpyxl.load_workbook(out_path) + vis = wb[pp.VISUAL_SHEET] + anchor = pp.find_anchor(vis, "compounds") + assert anchor is not None + flipped = None + for row_letter in "BCDEFG": + for col in COLS: + cell = vis.cell(row=anchor + 1 + ord(row_letter) - ord("A"), column=col + 1) + if cell.value == "Cpd 2": + flipped = cell + break + if flipped is not None: + break + assert flipped is not None + # give it the DMSO fill (wrong color) and expect verify to fail + dmso_cell = next(c for r in vis.iter_rows() for c in r if c.value == "DMSO") + flipped.fill = pp.fill_from_spec(pp.fill_spec(dmso_cell.fill)) + wb.save(out_path) + proc = run_cli("verify", "--input", str(out_path), "--layout", str(layout_path), + expect_ok=False) + assert proc.returncode != 0 + assert "color check: 1 mismatches" in proc.stdout + assert "does not match Cpd 2 legend" in proc.stderr + + +@pytest.mark.skipif(shutil.which("minizinc") is None or not (ASSETS / "plate-design.mzn").exists(), + reason="minizinc or vendored plaid assets not available") +def test_solve_end_to_end(workbook: Path, tmp_path: Path) -> None: + layout_path = tmp_path / "layout.json" + dzn_path = tmp_path / "plate.dzn" + solution_path = tmp_path / "solution.json" + out_path = tmp_path / "mini-plaid.xlsx" + run_cli("extract", "--input", str(workbook), "--output", str(layout_path)) + run_cli("dzn", "--layout", str(layout_path), "--output", str(dzn_path)) + run_cli("solve", "--mzn", str(ASSETS / "plate-design.mzn"), "--dzn", str(dzn_path), + "--layout", str(layout_path), "--output", str(solution_path), "--timeout", "30") + run_cli("apply", "--input", str(workbook), "--output", str(out_path), + "--layout", str(layout_path), "--solution", str(solution_path)) + proc = run_cli("verify", "--input", str(out_path), "--layout", str(layout_path)) + assert "VERIFY OK" in proc.stdout + # original workbook must be untouched by apply + before = workbook.read_bytes() + run_cli("apply", "--input", str(workbook), "--output", str(out_path), + "--layout", str(layout_path), "--solution", str(solution_path)) + assert workbook.read_bytes() == before \ No newline at end of file