%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import numpy as np
import pylab as pb
import GPy
import pandas
from pandas import read_csv

"""
Explanation: Partial Differential Equations and Gaussian Process with GPy
presented at the EBI BioPreDyn Course 'The Systems Biology Modelling Cycle'
Mu Niu, Neil Lawrence, 12th May 20014, University of Sheffield
The Spatio-Temporal Model
In this notebook we consider the latent force model paradigm (<a href="http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=6514873">Alvarez et al, 2013</a>, <a href="http://arxiv.org/abs/1107.2699">Alvarez et al, 2011</a>, <a href="http://jmlr.org/proceedings/papers/v5/alvarez09a/alvarez09a.pdf">Alvarez et al, 2009</a>) to drive spatio-temporal differential equation with a Gaussian process. Latent force models are differential equations whose initial conditions or driving forces are given by a stochastic process. Linear latent force models are linear (partial or ordinary) differential equations. If a linear latent force model is driven by a Gaussian process latent force, this describes a joint Gaussian process distribution across all the variables of interest. The Gaussian process covariance function encodes the relationships between the variables, as proscribed by the differential equation, through the covariance function. 
This direction of research is part of our ideas as to how to merge mechanistic and statistical models of data. It also maps onto our ideas about 'Gaussian processes over everything'. The covariance function interelates the protein and mRNA concentrations in the same model.
Differential Equation Model
A model of post-transcriptional processing is formulated to describe the spatio-temporal Drosophila
protein expression data (<a href="http://www.ploscompbiol.org/article/info%3Adoi%2F10.1371%2Fjournal.pcbi.1003281">Becker et al, 2013</a>, <a href="http://arxiv.org/abs/1107.2699">Alvarez et al, 2011</a>). Protein production is considered to be linearly dependent on the concentration of mRNA at an earlier time point. The model also allows for diffusion of protein between nuclei and linear protein decay. These processes are dependent on the diffusion parameter and the degradation rate of protein respectively.
\begin{equation}
 a \frac{\partial ^2 y_{x,t}}{\partial x^2} + b \frac{\partial y_{x,t}}{\partial t} + c y_{x,t}=  f_{x,t} 
\end{equation}
The coefficients $a$, $b$ and $c$ are unknown. In this study, we use Gaussian process with an exponentiated quadratic kernel as a prior over $y_{x,t}$ (protein). The kernel of $f_{x,t}$ (mRNA) is derived by applying the partial differential operator on the spatial-temporal kernel of protein. The multi-output Gaussian process are developed by combining the covariance matrix of mRNA and protein and their cross covariance.
End of explanation
"""


data = GPy.util.datasets.drosophila_knirps()

"""
Explanation: The spatio-temporal multi-output partial differential equation covariance functions have been developed with the kernel name GPy.kern.ODE_st in GPy. The inputs are one dimension spatial data, one dimension temporal data and one dimension index which is used to indicate $f$ and $y$.
End of explanation
"""


kern = GPy.kern.ODE_st(input_dim=3,
                       a=1., b=1., c=1.,
                       variance_Yx=1., variance_Yt=1.,
                       lengthscale_Yx=15.,
                       lengthscale_Yt=15.) 

"""
Explanation: The next thing to do is to compose the data ready for presentation to GPy. Here we need to use the time field as the input and a concatanation of the expression levels as the output. We need to provide a corresponding index for to each input to describe what output is being represented.
Now we set up the covariance function with the relevant parameters.
End of explanation
"""


data['X']
data['Y']

m = GPy.models.GPRegression(data['X'],data['Y'],kern) 

"""
Explanation: With the data correctly presented and the covariance function defined, we are ready to proceed with the Gaussian process regression.
End of explanation
"""


leng = data['X'].shape[0]
m.plot(fixed_inputs=[(2,0)], which_data_rows = slice(0,leng*2/2))
m.plot(fixed_inputs=[(2,1)], which_data_rows = slice(leng*2/2,leng*2))

"""
Explanation: Initial Fit
The initial value of $a$, $b$ and $c$ are 1. For these choices of covariance function parameters, we can plot the random field of $f$ and $y$ separately.
End of explanation
"""


m.optimize(messages=True)

"""
Explanation: Now we optimize the model, this will take a few minutes.
End of explanation
"""


print m

"""
Explanation: After optimization, the estimated value of $a$, $b$ and $c$ can be printed.
End of explanation
"""


m.plot(fixed_inputs=[(2,0)], which_data_rows = slice(0,leng*2/2))
pb.savefig("gene.pdf")
GPy.plotting.show(m.plot(fixed_inputs=[(2,1)], which_data_rows = slice(leng*2/2,leng*2)), filename='test')
pb.savefig("protein.pdf")

"""
Explanation: The plot of the random fields are plotted below.
End of explanation
"""


