from pomegranate import *
%pylab inline
import numpy as np

data = np.concatenate( (np.random.randn(250, 1) * 2.75 + 1.25, np.random.randn(500, 1) * 1.2 + 7.85) )
np.random.shuffle(data)
plt.hist( data, edgecolor='c', color='c', bins=20 )

"""
Explanation: General Mixture Models
author: Jacob Schreiber <br>
contact: jmschreiber91@gmail.com
It is frequently the case that the data you have is not explained by a single underlying distribution. If we want to try to recover the underlying distributions, we need to have a model which has multiple components. An example is the following data.
End of explanation
"""


d = GeneralMixtureModel( [NormalDistribution(2.5, 1), NormalDistribution(8, 1)] )

"""
Explanation: We can create our initial estimate of what this distribution is a General Mixture Model. This is a model which is comprised of multiple distributions, and weights on those distributions representing the prior probability of a point falling under that distribution given no knowledge of the point itself (defaults to equal). We can have univariate mixture models by using univariate distributions, or multivariate distributions by using multivariate distributions.
Gaussian Mixture Models
End of explanation
"""


labels = d.predict( data )
print labels[:5]
print "{} 1 labels, {} 0 labels".format( labels.sum(), labels.shape[0] - labels.sum() )

"""
Explanation: We can now predict the class labels of each point under this mixture.
End of explanation
"""


plt.hist( data[ labels == 0 ], edgecolor='r', color='r', bins=20 )
plt.hist( data[ labels == 1 ], edgecolor='c', color='c', bins=20 )

"""
Explanation: This is fairly close to the number of underlying points from each distribution, off by 17 in each label. We still don't know if the labels are accurate, just the number of labels.
End of explanation
"""


labels = d.predict_proba( data )
print labels[:5]
print labels.sum(axis=0)

"""
Explanation: It is slightly more difficult to update the underlying components of the model because we don't have labels indicating which point came from which distribution. We could try to use the labels inferred from the model. It seems to cleanly split it, but what if our initial estimate was not very good? It could be difficult to get a good update if we had a bad prior. 
Another possibility is to predict the probability of each point under each component, to get a softer estimate of the labels. Lets take a look.
End of explanation
"""


d.fit( data, verbose=True )

"""
Explanation: This is slightly closer to the truth, with 15.2 off instead of 17, around 10% closer.
This is the beginning of a common unsupervised training algorithm called <b>expectation maximization</b>. It has two steps, <b>expectation</b> and <b>maximization</b>. The <b>expectation</b> step involves what we just did--assigning weights based on the probability of each point being generated by each component. The next step, <b>maximization</b>, is maximizing the probability that the distribution generated these points but performing weighted MLE.
This process must be iterated until convergence. Sometimes this requires only a single update, but for overlapping distributions (such as this one) it can sometimes take many iterations.
End of explanation
"""


mu = np.arange(5)
cov = np.eye(5)

mgs = [ MultivariateGaussianDistribution( mu*i, cov ) for i in range(5) ]
gmm = GeneralMixtureModel( mgs )

data = numpy.random.randn(1000, 5) * 5
for i in range(5):
    data[i::5] += np.arange(5)*i

"""
Explanation: Multivariate Gaussian Mixture Models
We can do the same with multivariate distributions just as easily.
End of explanation
"""


for i in range(10):
    print "Point {}: logp {}".format( i, gmm.log_probability(data[i]) )

gmm.fit(data, verbose=True, stop_threshold=1)

"""
Explanation: Lets see how well some points fit under the mixture model.
End of explanation
"""


for i in range(10):
    print "Point {}: logp {}".format( i, gmm.log_probability(data[i]) )

"""
Explanation: Now lets see how well the previous points fit.
End of explanation
"""


d1 = DiscreteDistribution( {'A' : 0.25, 'C': 0.25, 'G' : 0.25, 'T': 0.25 } ) # Background
d2 = DiscreteDistribution( {'A' : 0.05, 'C': 0.45, 'G' : 0.45, 'T': 0.05 } ) # CG rich regions
gmm = GeneralMixtureModel( [d1, d2] )

seq = numpy.array(list('CGACATCTGACTACGGCGCGCCTACTACTTGATCGATACGGCGTCAGCGACGACGATGATCGGCATCAGTCACTAC'))
gmm.fit(seq)

print gmm.distributions
print
print numpy.exp(gmm.weights)

"""
Explanation: Looks like they're being fit significantly better than before! Training works.
Discrete Mixture Models
In addition to having general mixture models over continuous distributions, we can also have mixture models over discrete distributions. This is useful in many bioinformatics contexts, specifically sequence analysis. Lets use the toy analysis of trying to analyze CG island distribution.
The problem is the following; DNA is made up of long sequences the four canonical nucleotides, abbreviated 'A', 'C', 'G', and 'T'. These nucleotides are not distributed randomly, and there is significant amounts of structure in the sequence. A major field in bioinformatics is trying to interpret this structure. One structured element is the CG content, where the Cs and the Gs appear more commonly than in the background. If we want to try to determine CG percentages in these islands, we can use a mixture model.
End of explanation
"""


print gmm.to_json()

"""
Explanation: Looks like in this case the concept was sound, that there many CG rich regions, but our initial estimates of the  percentages were off. We can use a GMM like the one above to both identify and study the composition of these regions at the same time, updating the parameters of the distributions using expectation-maximization. We will go into a more complex way of dong this using HMMs in the next tutorial.
Serialization
General Mixture Models support serialization to JSONs using to_json() and from_json( json ). This is useful is you want to train a GMM on large amounts of data, taking a significant amount of time, and then use this model in the future without having to repeat this computationally intensive step.
End of explanation
"""


gmm_2 = GeneralMixtureModel.from_json( gmm.to_json() )
gmm_2.distributions

"""
Explanation: Not the prettiest thing to look at right now. However, we can easily load it and use it immediately.
End of explanation
"""


