Computes the completion and concordance of each samples.

    :param tped: the ``tped`` containing duplicated samples.
    :param tfam: the ``tfam`` containing duplicated samples.
    :param samples: the updated position of the samples in the tped containing
                    only duplicated samples.
    :param oldSamples: the original duplicated sample positions.
    :param prefix: the prefix of all the files.

    :type tped: :py:class:`numpy.array`
    :type tfam: :py:class:`numpy.array`
    :type samples: dict
    :type oldSamples: dict
    :type prefix: str

    :returns: a tuple containing the completion (:py:class:`numpy.array`) as
              first element, and the concordance (:py:class:`dict`) as last
              element.

    Reads the ``tped`` file and compute the completion for each duplicated
    samples and the pairwise concordance between duplicated samples.

    .. note::
        The completion and concordance computation excludes a markers if it's
        on chromosome 24 and if the sample is a female.

    .. note::
        A missing genotype is encoded by ``0``.

    .. note::
        No percentage is computed here, only the numbers. Percentages are
        computing in other functions: :py:func:`printStatistics`, for
        completion, and :py:func:`printConcordance`, for concordance.


    **Completion**

    Computes the completion of none zero values (where all genotypes of at
    least one duplicated sample are no call [*i.e.* ``0``]). The completion of
    sample :math:`i` (*i.e.* :math:`Comp_i`) is the number of genotypes
    that have a call divided by the total number of genotypes (the set
    :math:`G_i`):

    .. math::
        Comp_i = \\frac{||g \\in G_i\\textrm{ where }g \\neq 0||}{||G_i||}

    .. note::
        We consider a genotype as being missing if the sample is a male and if
        a marker on chromosome 23 or 24 is heterozygous.


    **Concordance**

    Computes the pairwise concordance between duplicated samples. For each
    marker, if both genotypes are not missing, we add one to the total number
    of compared markers. If both genotypes are the same, we add one to the
    number of concordant calls. We write the observed genotype difference in
    the file ``prefix.diff``. The concordance between sample :math:`i` and
    :math:`j` (*i.e.* :math:`Concordance_{i,j}`) is the number of genotypes
    that are equal divided by the total number of genotypes (excluding the no
    calls):

    .. math::
        Concordance_{i,j} = \\frac{
            ||g \\in G_i \\cup G_j \\textrm{ where } g_i = g_j \\neq 0||
        }{
            ||g \\in G_i \\cup G_j \\textrm{ where } g \\neq 0||
        }

def computeStatistics(tped, tfam, samples, oldSamples, prefix):
    """Computes the completion and concordance of each samples.

    :param tped: the ``tped`` containing duplicated samples.
    :param tfam: the ``tfam`` containing duplicated samples.
    :param samples: the updated position of the samples in the tped containing
                    only duplicated samples.
    :param oldSamples: the original duplicated sample positions.
    :param prefix: the prefix of all the files.

    :type tped: :py:class:`numpy.array`
    :type tfam: :py:class:`numpy.array`
    :type samples: dict
    :type oldSamples: dict
    :type prefix: str

    :returns: a tuple containing the completion (:py:class:`numpy.array`) as
              first element, and the concordance (:py:class:`dict`) as last
              element.

    Reads the ``tped`` file and compute the completion for each duplicated
    samples and the pairwise concordance between duplicated samples.

    .. note::
        The completion and concordance computation excludes a markers if it's
        on chromosome 24 and if the sample is a female.

    .. note::
        A missing genotype is encoded by ``0``.

    .. note::
        No percentage is computed here, only the numbers. Percentages are
        computing in other functions: :py:func:`printStatistics`, for
        completion, and :py:func:`printConcordance`, for concordance.


    **Completion**

    Computes the completion of none zero values (where all genotypes of at
    least one duplicated sample are no call [*i.e.* ``0``]). The completion of
    sample :math:`i` (*i.e.* :math:`Comp_i`) is the number of genotypes
    that have a call divided by the total number of genotypes (the set
    :math:`G_i`):

    .. math::
        Comp_i = \\frac{||g \\in G_i\\textrm{ where }g \\neq 0||}{||G_i||}

    .. note::
        We consider a genotype as being missing if the sample is a male and if
        a marker on chromosome 23 or 24 is heterozygous.


    **Concordance**

    Computes the pairwise concordance between duplicated samples. For each
    marker, if both genotypes are not missing, we add one to the total number
    of compared markers. If both genotypes are the same, we add one to the
    number of concordant calls. We write the observed genotype difference in
    the file ``prefix.diff``. The concordance between sample :math:`i` and
    :math:`j` (*i.e.* :math:`Concordance_{i,j}`) is the number of genotypes
    that are equal divided by the total number of genotypes (excluding the no
    calls):

    .. math::
        Concordance_{i,j} = \\frac{
            ||g \\in G_i \\cup G_j \\textrm{ where } g_i = g_j \\neq 0||
        }{
            ||g \\in G_i \\cup G_j \\textrm{ where } g \\neq 0||
        }

    """
    # The diff file containing the genotype differences between a pair of
    # duplicated samples
    diffOutput = None
    try:
        diffOutput = open(prefix + ".diff", "w")
    except IOError:
        msg = "%(prefix)s.diff: can't write file" % locals()
        raise ProgramError(msg)
    print >>diffOutput, "\t".join(["name", "famID", "indID", "dupIndex_1",
                                   "dupIndex_2", "genotype_1", "genotype_2"])

    # The completion data type
    completion = np.array([[0 for i in xrange(len(tped[0][4:]))],
                           [0 for i in xrange(len(tped[0][4:]))]])

    # The concordance data type
    concordance = {}
    for sampleID in samples.keys():
        nbDup = len(samples[sampleID])
        concordance[sampleID] = [
            np.asmatrix(np.zeros((nbDup, nbDup), dtype=int)),
            np.asmatrix(np.zeros((nbDup, nbDup), dtype=int)),
        ]

    #####################################################################
    # Add options for concordance only for hetero SNPs... (all samples) #
    #####################################################################

    for row in tped:
        genotype = row[4:]
        chromosome = row[0]
        snpName = row[1]
        for key, indexes in samples.iteritems():
            for i, index in enumerate(indexes):
                sex = tfam[oldSamples[key][i]][4]
                genotype1 = set(genotype[index].split(" "))

                # Updating the completion
                if not ((chromosome == "24") and (sex == "2")):
                    if (sex == "1") and (chromosome in ["23", "24"]):
                        if ("0" not in genotype1) and (len(genotype1) != 2):
                            completion[0][index] += 1
                    elif "0" not in genotype1:
                        completion[0][index] += 1
                    completion[1][index] += 1

                    # Updating the concordance
                    for j in xrange(i + 1, len(indexes)):
                        genotype2 = set(genotype[indexes[j]].split(" "))

# This is for man on chr 23 and 24, if heterozygous!!! #
# if (sex == "1") and (chromosome in ["23", "24"]):
#     if ("0" not in genotype1) and ("0" not in genotype2):
#         if (len(genotype1) != 2) and (len(genotype2) != 2):
#             concordance[key][1][i,j] += 1
#             concordance[key][1][j,i] += 1
#             if genotype1 == genotype2:
#                 concordance[key][0][i,j] += 1
#                 concordance[key][0][j,i] += 1
#             else:
#                 # Print to diff file
#                 toPrint = [snpName, key[0], key[1],
#                            str(index+1),
#                            str(indexes[j]+1)]
#                 toPrint.append(" ".join(genotype1))
#                 if len(genotype1) == 1:
#                     toPrint[-1] += " %s" % toPrint[-1]
#                 toPrint.append(" ".join(genotype2))
#                 if len(genotype2) == 1:
#                     toPrint[-1] += " %s" % toPrint[-1]
#                 print >>diffOutput, "\t".join(toPrint)
#
# elif ("0" not in genotype1) and ("0" not in genotype2):

                        if ("0" not in genotype1) and ("0" not in genotype2):
                            # Both have calls
                            concordance[key][1][i, j] += 1
                            concordance[key][1][j, i] += 1
                            if genotype1 == genotype2:
                                concordance[key][0][i, j] += 1
                                concordance[key][0][j, i] += 1
                            else:
                                # Print to diff file
                                toPrint = [snpName, key[0], key[1],
                                           str(index+1), str(indexes[j]+1)]
                                toPrint.append(" ".join(genotype1))
                                if len(genotype1) == 1:
                                    toPrint[-1] += " %s" % toPrint[-1]
                                toPrint.append(" ".join(genotype2))
                                if len(genotype2) == 1:
                                    toPrint[-1] += " %s" % toPrint[-1]
                                print >>diffOutput, "\t".join(toPrint)

    diffOutput.close()

    for key in concordance.iterkeys():
        for i in range(len(concordance[key][0])):
            concordance[key][0][i, i] = 1
            concordance[key][1][i, i] = 1

    return completion, concordance