% 
% Mendelian error detection in complex pedigrees in MiniZinc.
%
% M. Sanchez, S. de Givry, and T. Schiex
% Mendelian error detection in complex pedigrees using weighted constraint satisfaction techniques
% Constraints, 13(1):130-154, 2008

% The  problem is  to  identify genotype  incompatibilities in  pedigree
% data.  A  pedigree data consists in  a list of  individuals with their
% parental relationship.  For each  individual, there are two alleles at
% a given locus  of interest, one coming from the  father, the other one
% from the mother.  Zero or two unordered alleles are known depending on
% genotyping  data.    For  each  nuclear  family   (parents  and  their
% children), Mendelien rules define  the possible combination of alleles
% between parents and children: for  every child, its father allele must
% be identical to one allele of its father and its mother allele must be
% identical  to one  allele of  its mother.
% If the problem is inconsistent, then the goal is to find the
% minimum number  of genotyping data to  be removed so  that the problem
% becomes satisfiable.
% 
% This MiniZinc model was created by Simon de Givry, degivry@toulouse.inra.fr
%

include "globals.mzn"; 

int: num_individuals;
int: num_alleles;

array[1..num_individuals] of int: father;
array[1..num_individuals] of int: mother;
array[1..num_individuals] of int: sex;
array[1..num_individuals] of int: genotyping1;
array[1..num_individuals] of int: genotyping2;

% the assignments
array[1..num_individuals] of var 1..num_alleles: allele1;
array[1..num_individuals] of var 1..num_alleles: allele2;
var int: objective;

% solve minimize objective;
solve :: int_search(allele1 ++ allele2, first_fail, indomain_min, complete) minimize objective;

constraint
   % parsimony criterion: minimizes number of incorrect genotyping data
   objective = sum(j in 1..num_individuals where genotyping1[j] * genotyping2[j] > 0) (bool2int(allele1[j] != min(genotyping1[j],genotyping2[j]) \/ allele2[j] != max(genotyping1[j],genotyping2[j]))) 
               + sum(j in 1..num_individuals where genotyping1[j] > 0 /\ genotyping2[j] = 0) (bool2int(allele1[j] != genotyping1[j] /\ allele2[j] != genotyping1[j]))
               + sum(j in 1..num_individuals where genotyping1[j] = 0 /\ genotyping2[j] > 0) (bool2int(allele1[j] != genotyping2[j] /\ allele2[j] != genotyping2[j]))
   %
   /\ % hard binary constraints for breaking symmetries
   forall(j in 1..num_individuals) (
       allele1[j] <= allele2[j]
   )
   /\ % hard mendelian rule constraints
   forall(j in 1..num_individuals where father[j] * mother[j] > 0) (
       (allele1[j] = allele1[father[j]] /\ allele2[j] = allele1[mother[j]]) \/ 
       (allele2[j] = allele1[father[j]] /\ allele1[j] = allele1[mother[j]]) \/ 
       (allele1[j] = allele1[father[j]] /\ allele2[j] = allele2[mother[j]]) \/ 
       (allele2[j] = allele1[father[j]] /\ allele1[j] = allele2[mother[j]]) \/ 
       (allele1[j] = allele2[father[j]] /\ allele2[j] = allele1[mother[j]]) \/ 
       (allele2[j] = allele2[father[j]] /\ allele1[j] = allele1[mother[j]]) \/ 
       (allele1[j] = allele2[father[j]] /\ allele2[j] = allele2[mother[j]]) \/ 
       (allele2[j] = allele2[father[j]] /\ allele1[j] = allele2[mother[j]])
   )
   /\
   forall(j in 1..num_individuals where father[j] > 0 /\ mother[j] = 0) (
       (allele1[j] = allele1[father[j]]) \/ 
       (allele2[j] = allele1[father[j]]) \/ 
       (allele1[j] = allele2[father[j]]) \/ 
       (allele2[j] = allele2[father[j]])
   )
   /\
   forall(j in 1..num_individuals where father[j] = 0 /\ mother[j] > 0) (
       (allele1[j] = allele1[mother[j]]) \/ 
       (allele2[j] = allele1[mother[j]]) \/ 
       (allele1[j] = allele2[mother[j]]) \/ 
       (allele2[j] = allele2[mother[j]])
   )
%
% solution checker for connell.dzn: objective=1
%   /\
%   forall(j in 1..12) (
%   	   allele1[j] = [1, 2, 2, 1, 1, 2, 2, 1, 2, 2, 1, 2][j]
%   	   allele2[j] = [2, 3, 2, 2, 3, 2, 2, 2, 3, 2, 2, 3][j]
%   )
;

output
[
  "alleles:\n" ++ show(allele1) ++ "\n" ++ show(allele2) ++ "\n" ++
  "objective: " ++ show(objective)
];