import faps as fp
import numpy as np

np.random.seed(27) # this ensures you get exactly the same answers as I do.
allele_freqs = np.random.uniform(0.3,0.5, 50)
mypop        = fp.make_parents(4, allele_freqs, family_name='my_population')
progeny      = fp.make_sibships(mypop, 0, [1,2], 3, 'myprogeny')

"""
Explanation: Paternity arrays
Tom Ellis, March 2017
Paternity arrays are the what sibship clustering is built on in FAPS. They contain information about the probability that each candidate male is the father of each individual offspring. This information is stored in a paternityArray object, along with other related information. A paternityArray can either be imported directly, or created from genotype data.
This notebook will examine how to:

Create a paternityArray from marker data.
Examine what information it contains.
Read and write a paternityArray to disk, or import a custom paternityArray.

Once you have made your paternityArray, the [next step]((https://github.com/ellisztamas/faps/blob/master/docs/04%20Sibship%20clustering.ipynb) is to cluster the individuals in your array into full sibship groups.
Creating a paternityArray from genotype data
To create a paternityArray from genotype data we need to specficy genotypeArrays for the offspring, mothers and candidate males. Currently only biallelic SNP data are supported.
We will illustrate this with a small simulated example again with four adults and six offspring typed at 50 loci.
End of explanation
"""


mum_index = progeny.parent_index('mother', mypop.names) # positions in the mothers in the array of adults
mothers   = mypop.subset(mum_index) # genotypeArray of the mothers

"""
Explanation: We need to supply a genotypeArray for the mothers. This needs to have an entry for for every offspring, i.e. six replicates of the mother.
End of explanation
"""


error_rate = 0.0015

sample_af = mypop.allele_freqs()
patlik = fp.paternity_array(
    offspring = progeny,
    mothers = mothers,
    males= mypop,
    mu=error_rate)

"""
Explanation: To create the paternityArray we also need to supply information on the genotyping error rate (mu), and population allele frequencies. In toy example we know the error rate to be zero. However, in reality this will almost never be true, and moreover, sibship clustering becomes unstable when errors are zero. With that in mind, fp.paternity_array will throw an error when this is included.
For allele frequencies, we can either take the population allele frequencies defined above, or estimate them from the data, which will give slightly different answers. The function paternity_array creates an object of class paternityArray.
End of explanation
"""


print(patlik.candidates)
print(patlik.mothers)
print(patlik.offspring)

"""
Explanation: Note that if you try running the above cell but with mu=0, it will run, but throw a warning. The reason is that setting mu to zero tends to make sibship clustering unstable, so it automatically sets mu to a very small number.
paternityArray structure
A paternityArray inherits information about individuals from found in a genotypeArray. For example, labels of the candidates, mothers and offspring.
End of explanation
"""


patlik.lik_array

"""
Explanation: The most important part of the paternityArray is the likelihood array, which represent the log likelihood that each candidate male is the true father of each offspring individual. In this case it will be a 6x4 dimensional array with a row for each offspring and a column for each candidate.
End of explanation
"""


patlik.lik_absent

"""
Explanation: You can see that the log likelihoods of paternity for the first individual are much lower than the other candidates. This individual is the mother, so this makes sense. You can also see that the highest log likelihoods are in the columns for the real fathers (the 2nd column in rows one to three, and the third column in rows four to six).
The paternityArray also includes information that the true sire is not in the sample of candidate males. In this case this is not helpful, because we know sampling is complete, but in real examples is seldom the case. By default this is defined as the liklihood of generating the offspring genotypes given the known mothers genotype and alleles drawn from population allele frequencies. Here, values for the six offspring are higher than the likelihoods for the non-sires, indicating that they are no more likely to be the true sire than a random unrelated individual.
End of explanation
"""


patlik.prob_array.shape

"""
Explanation: This becomes more informative when we combine likelihoods about sampled and unsampled fathers. For fractional analyses we really want to know the probability that the father was unsampled vs sampled, and how probable it is that a single candidate is the true sire. To do this, patlik.lik_array is concatenate with values in patlik.lik_absent, then normalises the array so that values in each row sum to one. Printing the shape of the array demonstrates that we have gained a column.
End of explanation
"""


np.exp(patlik.prob_array).sum(axis=1)

"""
Explanation: If we sum the rows we see that they do indeed add up to one now. Probabilities are stored as log probabilities, so we have to exponentiate first.
End of explanation
"""


patlik.adjust_prob_array(selfing_rate=0)[:, :3]

"""
Explanation: Modifying a paternityArray
We can alter the information in patlik.prob_array to reflect different prior beliefs about the dataset. (In contrast, it's seldom a good idea to manipulate the likelihoods from genetic data contained in patlik.lik_array).
For example, often the mother is included in the sample of candidate males, either because you are using the same array for multiple families, or self-fertilisation is a biological possibility. In a lot of cases though the mother cannot simultaneously be the sperm/pollen donor, and it is necessary to set the rate of self-fertilisation to zero (the natural logarithm of zero is negative infinity). Here, only the first three columns are shown.
End of explanation
"""


patlik.prob_array = patlik.adjust_prob_array(selfing_rate=0)

"""
Explanation: The likelihoods for the mother have changed to 0 (negative infinity on the log scale). You can set any selfing rate between zero and one if you have a good idea of what the value should be and how much it varies. Otherwise it may be better to estimate the selfing rate from the data, or else estimate it some other way.
adjust_prob_array always refers back to the original patlik.lik_array and patlik.lik_absent, which remain unchanged. Calling adjust_prob_array will not alter the data stored for patlik.prob_array unless you assign it yourself.
End of explanation
"""


patlik.adjust_prob_array(purge = 'my_population_0')

"""
Explanation: You can also set likelihoods for particular individuals to zero manually. You might want to do this if you wanted to test the effects of incomplete sampling on your results, or if you had a good reason to suspect that some candidates could not possibly be the sire (for example, if the data are multigenerational, and the candidate was born after the offspring).
End of explanation
"""


patprob2 = patlik.adjust_prob_array(purge=0.4)
np.isinf(patprob2).mean(1) # proportion missing along each row.

"""
Explanation: This has removed the first individual (notice that this is identical to the previous example, because in this case the first individual is the mother). Alternatively you can supply a float between zero and one, which will be interpreted as a proportion of the candidates to be removed at random, which can be useful for simulations.
End of explanation
"""


np.exp(patprob2).sum(1)

patlik.adjust_prob_array(missing_parents=0.1)

"""
Explanation: You can specify the proportion $\theta$ of the population of candidate males which are missing with the option missing_parents. The likelihoods for non-sampled parents will be weighted by $\theta$, and likelihoods for sampled candidates by $1-\theta$.
Of course, rows still need to sum to one. Luckily adjust_prob_array does that automatically.
End of explanation
"""


patlik.adjust_prob_array(max_clashes=3)

"""
Explanation: You might want to remove candidates who have an a priori very low probability of paternity, for example to reduce the memory requirements of the paternityArray. One simple rule is to exclude any candidates with more than some arbritray number of loci with opposing homozygous genotypes relative to the offspring (you want to allow for a small number, in case there are genotyping errors). This is done with max_clashes.
End of explanation
"""


patlik.clashes

"""
Explanation: The option max_clashes refers back to a matrix that counts the number of such incompatibilities for each offspring-candidate pair. When you create a paternityArray from genotypeArray objects, this matrix is created automatically ad can be called with:
End of explanation
"""


fp.incompatibilities(mypop, progeny)

"""
Explanation: You can recreate this manually with:
End of explanation
"""


patlik.write('../data/mypatlik.csv')

"""
Explanation: Notice that this array has a row for each offspring, and a column for each candidate father. The first column is for the mother, which is why everything is zero.
Importing a paternityArray
Frequently you may wish to save an array and reload it. Otherwise, you may be working with a more exotic system than FAPS currently supports, such as microsatellite markers or a funky ploidy system. In this case you can create your own matrix of paternity likelihoods and import this directly as a paternityArray. Firstly, we can save the array we made before to disk by supplying a path to save to:
End of explanation
"""


patlik = fp.read_paternity_array(
    path = '../data/mypatlik.csv',
    mothers_col=1,
    likelihood_col=2)

"""
Explanation: We can reimport it again using read_paternity_array. This function is similar to the function for importing a genotypeArray, and the data need to have a specific structure:

Offspring names should be given in the first column
Names of the mothers are usually given in the second column.
If known for some reason, names of fathers can be given as well.
Likelihood information should be given to the right of columns indicating individual or parental names, with candidates' names in the column headers.
The final column should specify a likelihood that the true sire of an individual has not been sampled. Usually this is given as the likelihood of drawing the paternal alleles from population allele frequencies.
End of explanation
"""


fp.incompatibilities(mypop, progeny)

"""
Explanation: Of course, you can of course generate your own paternityArray and import it in the same way. This is especially useful if your study system has some specific marker type or genetic system not supported by FAPS.
One caveat with importing data is that the array of opposing homozygous loci is not imported automatically. You can either import this as a separate text file, or you can recreate this as above:
End of explanation
"""


