CODISINCEST(handles,KC)

function CODISINCEST(handles,KC) 

%each variable is inheredited by CODIS
global I M C HM HC ED WE matM Amb NL CODISpop
%I is the for...end variable that can be passed to Edwards'es and Wenk's functions
%M, C are the matrices of mother's, child's and alleged father's STRs data;
%HM, HC are the flag matrices of homozygosity (HM - mother; HC - child;) 
%HX(I) is 1 if the I-th locus is homozygous, 0 otherwise
%matM is the match-flag matrices: matM(I)=1 if there is a match between child's
%and mother's locus, 0 otherwise
%Amb is the ambigous match-flag matrix: if M=PQ C=PQ there is a match, but
%it is impossible to extablish which allele is inheredited allele from the 
%mother and which from the father, so (Amb(I)=1);
%for the I-th STR and selected population dataset (CODISpop matrix):
%p is the alleles frequencies array. Remember that in the CODIS GUI you can
%select the 0 allele (when the STR was not tested); so if you selected the J-th
%allele from the popupmenu, the correct frequency will be J-1.
%AIC relatedness coefficient

%Set the match-flag vectors:
%starting from the latin assumption "Mater semper certa est, pater numquam"
%the matM array is all 
matM=ones(18,1); Amb=zeros(18,1);
%set the Incest Indexes
WE=matM;
%sort M C vectors
M=sort(M,2); C=sort(C,2);
%highlight the homozygous loci
HM=M(:,1)==M(:,2); HC=C(:,1)==C(:,2);

NL=length(CODISpop);
for I=1:NL %there are 18 loci in the dataset
    % Check for errors: the child's and the mother's loci are required
    if (all(C(I,:)) && all(M(I,:)))
        if sum(ismember(M(I,:),C(I,:)))==0
            matM(I)=0; %switch off the flag and switch on CODIS2
        else %if there is the match between mother and child...
            %tidy up the loci
            %Check if M=PQ & C=PQ
            Amb(I)=(~HM(I) && ~HC(I) && isequal(M(I,:),C(I,:)));
            if HC(I) %the child's locus is homozygous
                %the mother's locus is heterozygous and it need to be swapped
                if ~HM(I) && M(I,1)==C(I,1) %M=PQ F=PP
                    CODISswap(1,0,0) %M=QP F=PP
                end
            else %the child's locus is heterozygous
                if HM(I) %the mother's locus is homozygous
                    %and the child's locus need to be swapped
                    if C(I,2)==M(I,2) %M=PP C=QP
                        CODISswap(0,1,0) %M=PP C=PQ
                    end
                else %the mother's locus is heterozygous
                    if ~Amb(I) %if we are not in the situation M=PQ C=PQ
                        if M(I,1)==C(I,1)
                            CODISswap(1,0,0) %swap only the mother's locus
                        elseif M(I,1)==C(I,2)
                            CODISswap(1,1,0) %swap both loci
                        elseif M(I,2)==C(I,2)
                            CODISswap(0,1,0) %swap only the child's locus
                        end
                    end
                end
            end
            WE(I)=Wenk(handles);
        end
    end
end
ED=Edwards(handles,KC);
return

function ED=Edwards(handles,KC) %#ok<INUSL>
global CODISpop HC C NL
ED=ones(NL,1);
Cidx=(C(:,1)~=0 & C(:,2)~=0);
HET=sum(HC(Cidx)==0);
hom=HC(Cidx)==1;
K=1;
for J=1:18
    if Cidx(J)
        if hom(K)
            eval(sprintf('locus=get(handles.text%i,''String'');',J))
            TF=sum(strcmp(locus,CODISpop(:,7)).*(1:1:NL)');
            x=CODISpop{TF,1}; y=CODISpop{TF,2};
            z=round(x.*y); y=y+2; z(C(J,:)-1)=z(C(J,:)-1)+1+HC(J); x=z./y;
            P=x(C(J,1)-1);
            ED(J)=(KC+((1-KC)*P))/P;
        end
        K=K+1;
    end
end
ED(19)=10^(sum(log10(ED(1:NL)))+log10(1-KC)*HET);
return

function Z=Wenk(handles) %#ok<INUSD>
global HC HM M I CODISpop Amb NL
eval(sprintf('locus=get(handles.text%i,''String'');',I))
TF=sum(strcmp(locus,CODISpop(:,7)).*(1:1:NL)');
x=CODISpop{TF,1}; y=CODISpop{TF,2};
z=round(x.*y); y=y+2; z(M(I,:)-1)=z(M(I,:)-1)+1+HM(I); x=z./y;
%Frequency of inherited maternal allele
Q=x(M(I,2)-1);
if HM(I) %mother's locus is homozygous
    if HC(I) %child's locus is homozygous
        Z=(0.5+(0.5*Q))/Q;
    else %child's locus is heterozygous
        Z=(0.5-(0.5*Q))/Q;
    end
else %mother's locus is heterozygous
    P=x(M(I,1)-1); %Frequency of not inherited maternal allele
    if HC(I) %child's locus is homozygous
        Z=(0.25+(0.25*(P+Q)))/(0.5*(P+Q));
    else %child's locus is heterozygous
        if Amb(I) %M=PQ C=PQ
            Z=(0.25+(0.25*(P+Q)))/(0.5*(P+Q));
        else %M=XP C=PQ
            Z=(0.5-(0.5*(P+Q)))/(1-(P+Q));
        end
    end
end
return