def computeStatistics(tped, tfam, snps):
    """Computes the completion and concordance of each SNPs.

    :param tped: a representation of the ``tped``.
    :param tfam: a representation of the ``tfam``
    :param snps: the position of the duplicated markers in the ``tped``.

    :type tped: numpy.array
    :type tfam: list
    :type snps: dict

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

    A marker's completion is compute using this formula (where :math:`G_i` is
    the set of genotypes for the marker :math:`i`):

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

    The pairwise concordance between duplicated markers is compute as follow
    (where :math:`G_i` and :math:`G_j` are the sets of genotypes for markers
    :math:`i` and :math:`j`, respectively):

    .. 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||
        }

    Hence, we only computes the numerators and denominators of the completion
    and concordance, for future reference.

    .. note::
        When the genotypes are not comparable, the function tries to flip one
        of the genotype to see if it becomes comparable.

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

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

    # The women and the no sex
    menIndex = np.where(tfam[:, 4] == "1")
    womenIndex = np.where(tfam[:, 4] == "2")
    noSexIndex = np.where(tfam[:, 4] == "0")

    for snpID, indexes in snps.iteritems():
        nbDup = len(indexes)
        currGenotypes = tped[indexes, 4:]
        chromosome, position = snpID

#         if chromosome == "24":
#             # Remove the heterozygous men
#             menToRemove = getIndexOfHeteroMen(currGenotypes, menIndex)
#             # Remove the women and the no sex
#             currGenotypes = np.delete(currGenotypes,
#                                        np.hstack((womenIndex, noSexIndex,
#                                                    menToRemove)), 1)
#         elif chromosome == "23":
#             # Remove the heterozygous men
#             menToRemove = getIndexOfHeteroMen(currGenotypes, menIndex)
#             # Remove the no sex
#             currGenotypes = np.delete(currGenotypes,
#                                        np.hstack((noSexIndex, menToRemove)),
#                                        1)

        for i in xrange(nbDup):
            # Compute completion here
            completion[0][indexes[i]] = len(
                np.where(currGenotypes[i] != "0 0")[0]
            )
            completion[1][indexes[i]] = len(currGenotypes[i])
            for j in xrange(i+1, nbDup):
                # Compute concordance here
                # Removing samples with at least one null genotype
                nullGenotypeIndexes = np.where(
                    np.any(currGenotypes[[i, j]] == "0 0", 0)
                )
                subGenotypes = np.delete(
                    currGenotypes,
                    nullGenotypeIndexes,
                    1,
                )

                # Finding the errors in the subseted genotypes
                errorIndexes = np.where(subGenotypes[i] != subGenotypes[j])[0]
                nbDiff = len(errorIndexes)

                for k in errorIndexes:
                    # Getting the genotypes
                    genotype1 = set(subGenotypes[i, k].split(" "))
                    genotype2 = set(subGenotypes[j, k].split(" "))

                    # Checking for flips
                    if len(genotype1) == len(genotype2):
                        # Both have the same number of different alleles,
                        # so they might be flipped
                        genotype2 = flipGenotype(genotype2)
                        if genotype1 == genotype2:
                            # The genotypes are equivalent after the flip
                            nbDiff -= 1

                # Updating the concordance
                nbTot = len(subGenotypes[i])
                concordance[snpID][0][i, j] = nbTot - nbDiff
                concordance[snpID][0][j, i] = nbTot - nbDiff
                if nbTot == 0:
                    # We will have a division by 0...
                    nbTot = 1
                concordance[snpID][1][i, j] = nbTot
                concordance[snpID][1][j, i] = nbTot

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

    return completion, concordance