Merge pull request #359 from SheffieldML/devel

Minor patch
This commit is contained in:
Max Zwiessele 2016-04-05 13:46:21 +01:00
commit 8d08da3348
15 changed files with 201 additions and 182 deletions

View file

@ -1 +1 @@
__version__ = "1.0.1" __version__ = "1.0.2"

View file

@ -437,15 +437,22 @@ class GP(Model):
warnings.warn("Wrong naming, use predict_wishart_embedding instead. Will be removed in future versions!", DeprecationWarning) warnings.warn("Wrong naming, use predict_wishart_embedding instead. Will be removed in future versions!", DeprecationWarning)
return self.predict_wishart_embedding(Xnew, kern, mean, covariance) return self.predict_wishart_embedding(Xnew, kern, mean, covariance)
def predict_magnification(self, Xnew, kern=None, mean=True, covariance=True): def predict_magnification(self, Xnew, kern=None, mean=True, covariance=True, dimensions=None):
""" """
Predict the magnification factor as Predict the magnification factor as
sqrt(det(G)) sqrt(det(G))
for each point N in Xnew for each point N in Xnew.
:param bool mean: whether to include the mean of the wishart embedding.
:param bool covariance: whether to include the covariance of the wishart embedding.
:param array-like dimensions: which dimensions of the input space to use [defaults to self.get_most_significant_input_dimensions()[:2]]
""" """
G = self.predict_wishard_embedding(Xnew, kern, mean, covariance) G = self.predict_wishard_embedding(Xnew, kern, mean, covariance)
if dimensions is None:
dimensions = self.get_most_significant_input_dimensions()[:2]
G = G[:, dimensions][:,:,dimensions]
from ..util.linalg import jitchol from ..util.linalg import jitchol
mag = np.empty(Xnew.shape[0]) mag = np.empty(Xnew.shape[0])
for n in range(Xnew.shape[0]): for n in range(Xnew.shape[0]):

View file

@ -14,10 +14,10 @@ import scipy as sp
class sde_RBF(RBF): class sde_RBF(RBF):
""" """
Class provide extra functionality to transfer this covariance function into Class provide extra functionality to transfer this covariance function into
SDE form. SDE form.
Radial Basis Function kernel: Radial Basis Function kernel:
.. math:: .. math::
@ -30,90 +30,90 @@ class sde_RBF(RBF):
Update gradient in the order in which parameters are represented in the Update gradient in the order in which parameters are represented in the
kernel kernel
""" """
self.variance.gradient = gradients[0] self.variance.gradient = gradients[0]
self.lengthscale.gradient = gradients[1] self.lengthscale.gradient = gradients[1]
def sde(self): def sde(self):
""" """
Return the state space representation of the covariance. Return the state space representation of the covariance.
""" """
N = 10# approximation order ( number of terms in exponent series expansion) N = 10# approximation order ( number of terms in exponent series expansion)
roots_rounding_decimals = 6 roots_rounding_decimals = 6
fn = np.math.factorial(N) fn = np.math.factorial(N)
kappa = 1.0/2.0/self.lengthscale**2 kappa = 1.0/2.0/self.lengthscale**2
Qc = np.array((self.variance*np.sqrt(np.pi/kappa)*fn*(4*kappa)**N,),) Qc = np.array((self.variance*np.sqrt(np.pi/kappa)*fn*(4*kappa)**N,),)
pp = np.zeros((2*N+1,)) # array of polynomial coefficients from higher power to lower pp = np.zeros((2*N+1,)) # array of polynomial coefficients from higher power to lower
for n in range(0, N+1): # (2N+1) - number of polynomial coefficients for n in range(0, N+1): # (2N+1) - number of polynomial coefficients
pp[2*(N-n)] = fn*(4.0*kappa)**(N-n)/np.math.factorial(n)*(-1)**n pp[2*(N-n)] = fn*(4.0*kappa)**(N-n)/np.math.factorial(n)*(-1)**n
pp = sp.poly1d(pp) pp = sp.poly1d(pp)
roots = sp.roots(pp) roots = sp.roots(pp)
neg_real_part_roots = roots[np.round(np.real(roots) ,roots_rounding_decimals) < 0] neg_real_part_roots = roots[np.round(np.real(roots) ,roots_rounding_decimals) < 0]
aa = sp.poly1d(neg_real_part_roots, r=True).coeffs aa = sp.poly1d(neg_real_part_roots, r=True).coeffs
F = np.diag(np.ones((N-1,)),1) F = np.diag(np.ones((N-1,)),1)
F[-1,:] = -aa[-1:0:-1] F[-1,:] = -aa[-1:0:-1]
L= np.zeros((N,1)) L= np.zeros((N,1))
L[N-1,0] = 1 L[N-1,0] = 1
H = np.zeros((1,N)) H = np.zeros((1,N))
H[0,0] = 1 H[0,0] = 1
# Infinite covariance: # Infinite covariance:
Pinf = sp.linalg.solve_lyapunov(F, -np.dot(L,np.dot( Qc[0,0],L.T))) Pinf = sp.linalg.solve_lyapunov(F, -np.dot(L,np.dot( Qc[0,0],L.T)))
Pinf = 0.5*(Pinf + Pinf.T) Pinf = 0.5*(Pinf + Pinf.T)
# Allocating space for derivatives # Allocating space for derivatives
dF = np.empty([F.shape[0],F.shape[1],2]) dF = np.empty([F.shape[0],F.shape[1],2])
dQc = np.empty([Qc.shape[0],Qc.shape[1],2]) dQc = np.empty([Qc.shape[0],Qc.shape[1],2])
dPinf = np.empty([Pinf.shape[0],Pinf.shape[1],2]) dPinf = np.empty([Pinf.shape[0],Pinf.shape[1],2])
# Derivatives: # Derivatives:
dFvariance = np.zeros(F.shape) dFvariance = np.zeros(F.shape)
dFlengthscale = np.zeros(F.shape) dFlengthscale = np.zeros(F.shape)
dFlengthscale[-1,:] = -aa[-1:0:-1]/self.lengthscale * np.arange(-N,0,1) dFlengthscale[-1,:] = -aa[-1:0:-1]/self.lengthscale * np.arange(-N,0,1)
dQcvariance = Qc/self.variance dQcvariance = Qc/self.variance
dQclengthscale = np.array(((self.variance*np.sqrt(2*np.pi)*fn*2**N*self.lengthscale**(-2*N)*(1-2*N,),))) dQclengthscale = np.array(((self.variance*np.sqrt(2*np.pi)*fn*2**N*self.lengthscale**(-2*N)*(1-2*N,),)))
dPinf_variance = Pinf/self.variance dPinf_variance = Pinf/self.variance
lp = Pinf.shape[0] lp = Pinf.shape[0]
coeff = np.arange(1,lp+1).reshape(lp,1) + np.arange(1,lp+1).reshape(1,lp) - 2 coeff = np.arange(1,lp+1).reshape(lp,1) + np.arange(1,lp+1).reshape(1,lp) - 2
coeff[np.mod(coeff,2) != 0] = 0 coeff[np.mod(coeff,2) != 0] = 0
dPinf_lengthscale = -1/self.lengthscale*Pinf*coeff dPinf_lengthscale = -1/self.lengthscale*Pinf*coeff
dF[:,:,0] = dFvariance dF[:,:,0] = dFvariance
dF[:,:,1] = dFlengthscale dF[:,:,1] = dFlengthscale
dQc[:,:,0] = dQcvariance dQc[:,:,0] = dQcvariance
dQc[:,:,1] = dQclengthscale dQc[:,:,1] = dQclengthscale
dPinf[:,:,0] = dPinf_variance dPinf[:,:,0] = dPinf_variance
dPinf[:,:,1] = dPinf_lengthscale dPinf[:,:,1] = dPinf_lengthscale
P0 = Pinf.copy() P0 = Pinf.copy()
dP0 = dPinf.copy() dP0 = dPinf.copy()
# Benefits of this are not very sound. Helps only in one case: # Benefits of this are not very sound. Helps only in one case:
# SVD Kalman + RBF kernel # SVD Kalman + RBF kernel
import GPy.models.state_space_main as ssm import GPy.models.state_space_main as ssm
(F, L, Qc, H, Pinf, P0, dF, dQc, dPinf,dP0, T) = ssm.balance_ss_model(F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0 ) (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf,dP0, T) = ssm.balance_ss_model(F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0 )
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0) return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
class sde_Exponential(Exponential): class sde_Exponential(Exponential):
""" """
Class provide extra functionality to transfer this covariance function into Class provide extra functionality to transfer this covariance function into
SDE form. SDE form.
Exponential kernel: Exponential kernel:
.. math:: .. math::
@ -121,53 +121,53 @@ class sde_Exponential(Exponential):
k(r) = \sigma^2 \exp \\bigg(- \\frac{1}{2} r \\bigg) \\ \\ \\ \\ \text{ where } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} } k(r) = \sigma^2 \exp \\bigg(- \\frac{1}{2} r \\bigg) \\ \\ \\ \\ \text{ where } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} }
""" """
def sde_update_gradient_full(self, gradients): def sde_update_gradient_full(self, gradients):
""" """
Update gradient in the order in which parameters are represented in the Update gradient in the order in which parameters are represented in the
kernel kernel
""" """
self.variance.gradient = gradients[0] self.variance.gradient = gradients[0]
self.lengthscale.gradient = gradients[1] self.lengthscale.gradient = gradients[1]
def sde(self): def sde(self):
""" """
Return the state space representation of the covariance. Return the state space representation of the covariance.
""" """
variance = float(self.variance.values) variance = float(self.variance.values)
lengthscale = float(self.lengthscale) lengthscale = float(self.lengthscale)
F = np.array(((-1.0/lengthscale,),)) F = np.array(((-1.0/lengthscale,),))
L = np.array(((1.0,),)) L = np.array(((1.0,),))
Qc = np.array( ((2.0*variance/lengthscale,),) ) Qc = np.array( ((2.0*variance/lengthscale,),) )
H = np.array(((1.0,),)) H = np.array(((1.0,),))
Pinf = np.array(((variance,),)) Pinf = np.array(((variance,),))
P0 = Pinf.copy() P0 = Pinf.copy()
dF = np.zeros((1,1,2)); dF = np.zeros((1,1,2));
dQc = np.zeros((1,1,2)); dQc = np.zeros((1,1,2));
dPinf = np.zeros((1,1,2)); dPinf = np.zeros((1,1,2));
dF[:,:,0] = 0.0 dF[:,:,0] = 0.0
dF[:,:,1] = 1.0/lengthscale**2 dF[:,:,1] = 1.0/lengthscale**2
dQc[:,:,0] = 2.0/lengthscale dQc[:,:,0] = 2.0/lengthscale
dQc[:,:,1] = -2.0*variance/lengthscale**2 dQc[:,:,1] = -2.0*variance/lengthscale**2
dPinf[:,:,0] = 1.0 dPinf[:,:,0] = 1.0
dPinf[:,:,1] = 0.0 dPinf[:,:,1] = 0.0
dP0 = dPinf.copy() dP0 = dPinf.copy()
return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0) return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
class sde_RatQuad(RatQuad): class sde_RatQuad(RatQuad):
""" """
Class provide extra functionality to transfer this covariance function into Class provide extra functionality to transfer this covariance function into
SDE form. SDE form.
Rational Quadratic kernel: Rational Quadratic kernel:
.. math:: .. math::
@ -177,16 +177,16 @@ class sde_RatQuad(RatQuad):
""" """
def sde(self): def sde(self):
""" """
Return the state space representation of the covariance. Return the state space representation of the covariance.
""" """
assert False, 'Not Implemented' assert False, 'Not Implemented'
# Params to use: # Params to use:
# self.lengthscale # self.lengthscale
# self.variance # self.variance
#self.power #self.power
#return (F, L, Qc, H, Pinf, dF, dQc, dPinf) #return (F, L, Qc, H, Pinf, dF, dQc, dPinf)

View file

@ -15,52 +15,43 @@
# #
import numpy as np import numpy as np
from scipy import linalg
from scipy import stats from scipy import stats
from ..core import Model
from .. import kern
#from GPy.plotting.matplot_dep.models_plots import gpplot
#from GPy.plotting.matplot_dep.base_plots import x_frame1D
#from GPy.plotting.matplot_dep import Tango
#import pylab as pb
from GPy.core.parameterization.param import Param
import GPy
from .. import likelihoods from .. import likelihoods
#from . import state_space_setup as ss_setup
from ..core import Model
from . import state_space_main as ssm from . import state_space_main as ssm
from . import state_space_setup as ss_setup from . import state_space_setup as ss_setup
class StateSpace(Model): class StateSpace(Model):
def __init__(self, X, Y, kernel=None, noise_var=1.0, kalman_filter_type = 'regular', use_cython = False, name='StateSpace'): def __init__(self, X, Y, kernel=None, noise_var=1.0, kalman_filter_type = 'regular', use_cython = False, name='StateSpace'):
super(StateSpace, self).__init__(name=name) super(StateSpace, self).__init__(name=name)
if len(X.shape) == 1: if len(X.shape) == 1:
X = np.atleast_2d(X).T X = np.atleast_2d(X).T
self.num_data, input_dim = X.shape self.num_data, self.input_dim = X.shape
if len(Y.shape) == 1: if len(Y.shape) == 1:
Y = np.atleast_2d(Y).T Y = np.atleast_2d(Y).T
assert input_dim==1, "State space methods are only for 1D data" assert self.input_dim==1, "State space methods are only for 1D data"
if len(Y.shape)==2: if len(Y.shape)==2:
num_data_Y, self.output_dim = Y.shape num_data_Y, self.output_dim = Y.shape
ts_number = None ts_number = None
elif len(Y.shape)==3: elif len(Y.shape)==3:
num_data_Y, self.output_dim, ts_number = Y.shape num_data_Y, self.output_dim, ts_number = Y.shape
self.ts_number = ts_number self.ts_number = ts_number
assert num_data_Y == self.num_data, "X and Y data don't match" assert num_data_Y == self.num_data, "X and Y data don't match"
assert self.output_dim == 1, "State space methods are for single outputs only" assert self.output_dim == 1, "State space methods are for single outputs only"
self.kalman_filter_type = kalman_filter_type self.kalman_filter_type = kalman_filter_type
#self.kalman_filter_type = 'svd' # temp test #self.kalman_filter_type = 'svd' # temp test
ss_setup.use_cython = use_cython ss_setup.use_cython = use_cython
#import pdb; pdb.set_trace() #import pdb; pdb.set_trace()
global ssm global ssm
#from . import state_space_main as ssm #from . import state_space_main as ssm
if (ssm.cython_code_available) and (ssm.use_cython != ss_setup.use_cython): if (ssm.cython_code_available) and (ssm.use_cython != ss_setup.use_cython):
@ -72,13 +63,13 @@ class StateSpace(Model):
# Noise variance # Noise variance
self.likelihood = likelihoods.Gaussian(variance=noise_var) self.likelihood = likelihoods.Gaussian(variance=noise_var)
# Default kernel # Default kernel
if kernel is None: if kernel is None:
raise ValueError("State-Space Model: the kernel must be provided.") raise ValueError("State-Space Model: the kernel must be provided.")
else: else:
self.kern = kernel self.kern = kernel
self.link_parameter(self.kern) self.link_parameter(self.kern)
self.link_parameter(self.likelihood) self.link_parameter(self.likelihood)
self.posterior = None self.posterior = None
@ -92,14 +83,14 @@ class StateSpace(Model):
""" """
Parameters have now changed Parameters have now changed
""" """
#np.set_printoptions(16) #np.set_printoptions(16)
#print(self.param_array) #print(self.param_array)
#import pdb; pdb.set_trace() #import pdb; pdb.set_trace()
# Get the model matrices from the kernel # Get the model matrices from the kernel
(F,L,Qc,H,P_inf, P0, dFt,dQct,dP_inft, dP0t) = self.kern.sde() (F,L,Qc,H,P_inf, P0, dFt,dQct,dP_inft, dP0t) = self.kern.sde()
# necessary parameters # necessary parameters
measurement_dim = self.output_dim measurement_dim = self.output_dim
grad_params_no = dFt.shape[2]+1 # we also add measurement noise as a parameter grad_params_no = dFt.shape[2]+1 # we also add measurement noise as a parameter
@ -109,30 +100,30 @@ class StateSpace(Model):
dQc = np.zeros([dQct.shape[0],dQct.shape[1],grad_params_no]) dQc = np.zeros([dQct.shape[0],dQct.shape[1],grad_params_no])
dP_inf = np.zeros([dP_inft.shape[0],dP_inft.shape[1],grad_params_no]) dP_inf = np.zeros([dP_inft.shape[0],dP_inft.shape[1],grad_params_no])
dP0 = np.zeros([dP0t.shape[0],dP0t.shape[1],grad_params_no]) dP0 = np.zeros([dP0t.shape[0],dP0t.shape[1],grad_params_no])
# Assign the values for the kernel function # Assign the values for the kernel function
dF[:,:,:-1] = dFt dF[:,:,:-1] = dFt
dQc[:,:,:-1] = dQct dQc[:,:,:-1] = dQct
dP_inf[:,:,:-1] = dP_inft dP_inf[:,:,:-1] = dP_inft
dP0[:,:,:-1] = dP0t dP0[:,:,:-1] = dP0t
# The sigma2 derivative # The sigma2 derivative
dR = np.zeros([measurement_dim,measurement_dim,grad_params_no]) dR = np.zeros([measurement_dim,measurement_dim,grad_params_no])
dR[:,:,-1] = np.eye(measurement_dim) dR[:,:,-1] = np.eye(measurement_dim)
# Balancing # Balancing
#(F,L,Qc,H,P_inf,P0, dF,dQc,dP_inf,dP0) = ssm.balance_ss_model(F,L,Qc,H,P_inf,P0, dF,dQc,dP_inf, dP0) #(F,L,Qc,H,P_inf,P0, dF,dQc,dP_inf,dP0) = ssm.balance_ss_model(F,L,Qc,H,P_inf,P0, dF,dQc,dP_inf, dP0)
# Use the Kalman filter to evaluate the likelihood # Use the Kalman filter to evaluate the likelihood
grad_calc_params = {} grad_calc_params = {}
grad_calc_params['dP_inf'] = dP_inf grad_calc_params['dP_inf'] = dP_inf
grad_calc_params['dF'] = dF grad_calc_params['dF'] = dF
grad_calc_params['dQc'] = dQc grad_calc_params['dQc'] = dQc
grad_calc_params['dR'] = dR grad_calc_params['dR'] = dR
grad_calc_params['dP_init'] = dP0 grad_calc_params['dP_init'] = dP0
kalman_filter_type = self.kalman_filter_type kalman_filter_type = self.kalman_filter_type
# The following code is required because sometimes the shapes of self.Y # The following code is required because sometimes the shapes of self.Y
# becomes 3D even though is must be 2D. The reason is undescovered. # becomes 3D even though is must be 2D. The reason is undescovered.
Y = self.Y Y = self.Y
@ -140,63 +131,63 @@ class StateSpace(Model):
Y.shape = (self.num_data,1) Y.shape = (self.num_data,1)
else: else:
Y.shape = (self.num_data,1,self.ts_number) Y.shape = (self.num_data,1,self.ts_number)
(filter_means, filter_covs, log_likelihood, (filter_means, filter_covs, log_likelihood,
grad_log_likelihood,SmootherMatrObject) = ssm.ContDescrStateSpace.cont_discr_kalman_filter(F,L,Qc,H, grad_log_likelihood,SmootherMatrObject) = ssm.ContDescrStateSpace.cont_discr_kalman_filter(F,L,Qc,H,
float(self.Gaussian_noise.variance),P_inf,self.X,Y,m_init=None, float(self.Gaussian_noise.variance),P_inf,self.X,Y,m_init=None,
P_init=P0, p_kalman_filter_type = kalman_filter_type, calc_log_likelihood=True, P_init=P0, p_kalman_filter_type = kalman_filter_type, calc_log_likelihood=True,
calc_grad_log_likelihood=True, calc_grad_log_likelihood=True,
grad_params_no=grad_params_no, grad_params_no=grad_params_no,
grad_calc_params=grad_calc_params) grad_calc_params=grad_calc_params)
if np.any( np.isfinite(log_likelihood) == False): if np.any( np.isfinite(log_likelihood) == False):
#import pdb; pdb.set_trace() #import pdb; pdb.set_trace()
print("State-Space: NaN valkues in the log_likelihood") print("State-Space: NaN valkues in the log_likelihood")
if np.any( np.isfinite(grad_log_likelihood) == False): if np.any( np.isfinite(grad_log_likelihood) == False):
#import pdb; pdb.set_trace() #import pdb; pdb.set_trace()
print("State-Space: NaN valkues in the grad_log_likelihood") print("State-Space: NaN valkues in the grad_log_likelihood")
#print(grad_log_likelihood) #print(grad_log_likelihood)
grad_log_likelihood_sum = np.sum(grad_log_likelihood,axis=1) grad_log_likelihood_sum = np.sum(grad_log_likelihood,axis=1)
grad_log_likelihood_sum.shape = (grad_log_likelihood_sum.shape[0],1) grad_log_likelihood_sum.shape = (grad_log_likelihood_sum.shape[0],1)
self._log_marginal_likelihood = np.sum( log_likelihood,axis=1 ) self._log_marginal_likelihood = np.sum( log_likelihood,axis=1 )
self.likelihood.update_gradients(grad_log_likelihood_sum[-1,0]) self.likelihood.update_gradients(grad_log_likelihood_sum[-1,0])
self.kern.sde_update_gradient_full(grad_log_likelihood_sum[:-1,0]) self.kern.sde_update_gradient_full(grad_log_likelihood_sum[:-1,0])
def log_likelihood(self): def log_likelihood(self):
return self._log_marginal_likelihood return self._log_marginal_likelihood
def _raw_predict(self, Xnew=None, Ynew=None, filteronly=False): def _raw_predict(self, Xnew=None, Ynew=None, filteronly=False, **kw):
""" """
Performs the actual prediction for new X points. Performs the actual prediction for new X points.
Inner function. It is called only from inside this class. Inner function. It is called only from inside this class.
Input: Input:
--------------------- ---------------------
Xnews: vector or (n_points,1) matrix Xnews: vector or (n_points,1) matrix
New time points where to evaluate predictions. New time points where to evaluate predictions.
Ynews: (n_train_points, ts_no) matrix Ynews: (n_train_points, ts_no) matrix
This matrix can substitude the original training points (in order This matrix can substitude the original training points (in order
to use only the parameters of the model). to use only the parameters of the model).
filteronly: bool filteronly: bool
Use only Kalman Filter for prediction. In this case the output does Use only Kalman Filter for prediction. In this case the output does
not coincide with corresponding Gaussian process. not coincide with corresponding Gaussian process.
Output: Output:
-------------------- --------------------
m: vector m: vector
Mean prediction Mean prediction
V: vector V: vector
Variance in every point Variance in every point
""" """
# Set defaults # Set defaults
if Ynew is None: if Ynew is None:
Ynew = self.Y Ynew = self.Y
@ -209,8 +200,8 @@ class StateSpace(Model):
else: else:
X = self.X X = self.X
Y = Ynew Y = Ynew
predict_only_training = True predict_only_training = True
# Sort the matrix (save the order) # Sort the matrix (save the order)
_, return_index, return_inverse = np.unique(X,True,True) _, return_index, return_inverse = np.unique(X,True,True)
X = X[return_index] # TODO they are not used X = X[return_index] # TODO they are not used
@ -218,37 +209,37 @@ class StateSpace(Model):
# Get the model matrices from the kernel # Get the model matrices from the kernel
(F,L,Qc,H,P_inf, P0, dF,dQc,dP_inf,dP0) = self.kern.sde() (F,L,Qc,H,P_inf, P0, dF,dQc,dP_inf,dP0) = self.kern.sde()
state_dim = F.shape[0] state_dim = F.shape[0]
#Y = self.Y[:, 0,0] #Y = self.Y[:, 0,0]
# Run the Kalman filter # Run the Kalman filter
#import pdb; pdb.set_trace() #import pdb; pdb.set_trace()
kalman_filter_type = self.kalman_filter_type kalman_filter_type = self.kalman_filter_type
(M, P, log_likelihood, (M, P, log_likelihood,
grad_log_likelihood,SmootherMatrObject) = ssm.ContDescrStateSpace.cont_discr_kalman_filter( grad_log_likelihood,SmootherMatrObject) = ssm.ContDescrStateSpace.cont_discr_kalman_filter(
F,L,Qc,H,float(self.Gaussian_noise.variance),P_inf,X,Y,m_init=None, F,L,Qc,H,float(self.Gaussian_noise.variance),P_inf,X,Y,m_init=None,
P_init=P0, p_kalman_filter_type = kalman_filter_type, P_init=P0, p_kalman_filter_type = kalman_filter_type,
calc_log_likelihood=False, calc_log_likelihood=False,
calc_grad_log_likelihood=False) calc_grad_log_likelihood=False)
# (filter_means, filter_covs, log_likelihood, # (filter_means, filter_covs, log_likelihood,
# grad_log_likelihood,SmootherMatrObject) = ssm.ContDescrStateSpace.cont_discr_kalman_filter(F,L,Qc,H, # grad_log_likelihood,SmootherMatrObject) = ssm.ContDescrStateSpace.cont_discr_kalman_filter(F,L,Qc,H,
# float(self.Gaussian_noise.variance),P_inf,self.X,self.Y,m_init=None, # float(self.Gaussian_noise.variance),P_inf,self.X,self.Y,m_init=None,
# P_init=P0, p_kalman_filter_type = kalman_filter_type, calc_log_likelihood=True, # P_init=P0, p_kalman_filter_type = kalman_filter_type, calc_log_likelihood=True,
# calc_grad_log_likelihood=True, # calc_grad_log_likelihood=True,
# grad_params_no=grad_params_no, # grad_params_no=grad_params_no,
# grad_calc_params=grad_calc_params) # grad_calc_params=grad_calc_params)
# Run the Rauch-Tung-Striebel smoother # Run the Rauch-Tung-Striebel smoother
if not filteronly: if not filteronly:
(M, P) = ssm.ContDescrStateSpace.cont_discr_rts_smoother(state_dim, M, P, (M, P) = ssm.ContDescrStateSpace.cont_discr_rts_smoother(state_dim, M, P,
p_dynamic_callables=SmootherMatrObject, X=X, F=F,L=L,Qc=Qc) p_dynamic_callables=SmootherMatrObject, X=X, F=F,L=L,Qc=Qc)
# remove initial values # remove initial values
M = M[1:,:,:] M = M[1:,:,:]
P = P[1:,:,:] P = P[1:,:,:]
# Put the data back in the original order # Put the data back in the original order
M = M[return_inverse,:,:] M = M[return_inverse,:,:]
P = P[return_inverse,:,:] P = P[return_inverse,:,:]
@ -257,40 +248,41 @@ class StateSpace(Model):
if not predict_only_training: if not predict_only_training:
M = M[self.num_data:,:,:] M = M[self.num_data:,:,:]
P = P[self.num_data:,:,:] P = P[self.num_data:,:,:]
# Calculate the mean and variance # Calculate the mean and variance
# after einsum m has dimension in 3D (sample_num, dim_no,time_series_no) # after einsum m has dimension in 3D (sample_num, dim_no,time_series_no)
m = np.einsum('ijl,kj', M, H)# np.dot(M,H.T) m = np.einsum('ijl,kj', M, H)# np.dot(M,H.T)
m.shape = (m.shape[0], m.shape[1]) # remove the third dimension m.shape = (m.shape[0], m.shape[1]) # remove the third dimension
V = np.einsum('ij,ajk,kl', H, P, H.T) V = np.einsum('ij,ajk,kl', H, P, H.T)
V.shape = (V.shape[0], V.shape[1]) # remove the third dimension V.shape = (V.shape[0], V.shape[1]) # remove the third dimension
# Return the posterior of the state # Return the posterior of the state
return (m, V) return (m, V)
def predict(self, Xnew=None, filteronly=False): def predict(self, Xnew=None, filteronly=False, include_likelihood=True, **kw):
# Run the Kalman filter to get the state # Run the Kalman filter to get the state
(m, V) = self._raw_predict(Xnew,filteronly=filteronly) (m, V) = self._raw_predict(Xnew,filteronly=filteronly)
# Add the noise variance to the state variance # Add the noise variance to the state variance
V += float(self.Gaussian_noise.variance) if include_likelihood:
V += float(self.likelihood.variance)
# Lower and upper bounds # Lower and upper bounds
lower = m - 2*np.sqrt(V) #lower = m - 2*np.sqrt(V)
upper = m + 2*np.sqrt(V) #upper = m + 2*np.sqrt(V)
# Return mean and variance # Return mean and variance
return (m, V, lower, upper) return m, V
def predict_quantiles(self, Xnew=None, quantiles=(2.5, 97.5)): def predict_quantiles(self, Xnew=None, quantiles=(2.5, 97.5), **kw):
mu, var = self._raw_predict(Xnew) mu, var = self._raw_predict(Xnew)
#import pdb; pdb.set_trace() #import pdb; pdb.set_trace()
return [stats.norm.ppf(q/100.)*np.sqrt(var + float(self.Gaussian_noise.variance)) + mu for q in quantiles] return [stats.norm.ppf(q/100.)*np.sqrt(var + float(self.Gaussian_noise.variance)) + mu for q in quantiles]
# def plot(self, plot_limits=None, levels=20, samples=0, fignum=None, # def plot(self, plot_limits=None, levels=20, samples=0, fignum=None,
# ax=None, resolution=None, plot_raw=False, plot_filter=False, # ax=None, resolution=None, plot_raw=False, plot_filter=False,
# linecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue']): # linecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue']):
@ -399,8 +391,8 @@ class StateSpace(Model):
# #
# # Return trajectory # # Return trajectory
# return Y # return Y
# #
# #
# def simulate(self,F,L,Qc,Pinf,X,size=1): # def simulate(self,F,L,Qc,Pinf,X,size=1):
# # Simulate a trajectory using the state space model # # Simulate a trajectory using the state space model
# #

View file

@ -52,6 +52,17 @@ def inject_plotting():
GP.plot_f = gpy_plot.gp_plots.plot_f GP.plot_f = gpy_plot.gp_plots.plot_f
GP.plot_magnification = gpy_plot.latent_plots.plot_magnification GP.plot_magnification = gpy_plot.latent_plots.plot_magnification
from ..models import StateSpace
StateSpace.plot_data = gpy_plot.data_plots.plot_data
StateSpace.plot_data_error = gpy_plot.data_plots.plot_data_error
StateSpace.plot_errorbars_trainset = gpy_plot.data_plots.plot_errorbars_trainset
StateSpace.plot_mean = gpy_plot.gp_plots.plot_mean
StateSpace.plot_confidence = gpy_plot.gp_plots.plot_confidence
StateSpace.plot_density = gpy_plot.gp_plots.plot_density
StateSpace.plot_samples = gpy_plot.gp_plots.plot_samples
StateSpace.plot = gpy_plot.gp_plots.plot
StateSpace.plot_f = gpy_plot.gp_plots.plot_f
from ..core import SparseGP from ..core import SparseGP
SparseGP.plot_inducing = gpy_plot.data_plots.plot_inducing SparseGP.plot_inducing = gpy_plot.data_plots.plot_inducing

View file

@ -190,6 +190,7 @@ def scatter_label_generator(labels, X, visible_dims, marker=None):
x = X[index, input_1] x = X[index, input_1]
y = X[index, input_2] y = X[index, input_2]
z = X[index, input_3] z = X[index, input_3]
yield x, y, z, this_label, index, m yield x, y, z, this_label, index, m
def subsample_X(X, labels, num_samples=1000): def subsample_X(X, labels, num_samples=1000):

View file

@ -131,14 +131,15 @@ class PlotlyPlots(AbstractPlottingLibrary):
#not matplotlib marker #not matplotlib marker
pass pass
marker_kwargs = marker_kwargs or {} marker_kwargs = marker_kwargs or {}
marker_kwargs.setdefault('symbol', marker) if 'symbol' not in marker_kwargs:
marker_kwargs['symbol'] = marker
if Z is not None: if Z is not None:
return Scatter3d(x=X, y=Y, z=Z, mode='markers', return Scatter3d(x=X, y=Y, z=Z, mode='markers',
showlegend=label is not None, showlegend=label is not None,
marker=Marker(color=color, colorscale=cmap, **marker_kwargs), marker=Marker(color=color, colorscale=cmap, **marker_kwargs),
name=label, **kwargs) name=label, **kwargs)
return Scatter(x=X, y=Y, mode='markers', showlegend=label is not None, return Scatter(x=X, y=Y, mode='markers', showlegend=label is not None,
marker=Marker(color=color, colorscale=cmap, **marker_kwargs or {}), marker=Marker(color=color, colorscale=cmap, **marker_kwargs),
name=label, **kwargs) name=label, **kwargs)
def plot(self, ax, X, Y, Z=None, color=None, label=None, line_kwargs=None, **kwargs): def plot(self, ax, X, Y, Z=None, color=None, label=None, line_kwargs=None, **kwargs):

Binary file not shown.

Before

Width:  |  Height:  |  Size: 186 KiB

After

Width:  |  Height:  |  Size: 191 KiB

Before After
Before After

Binary file not shown.

Before

Width:  |  Height:  |  Size: 180 KiB

After

Width:  |  Height:  |  Size: 185 KiB

Before After
Before After

View file

@ -730,6 +730,7 @@ class GradientTests(np.testing.TestCase):
self.assertTrue( np.allclose(var1, var2) ) self.assertTrue( np.allclose(var1, var2) )
def test_gp_VGPC(self): def test_gp_VGPC(self):
np.random.seed(10)
num_obs = 25 num_obs = 25
X = np.random.randint(0, 140, num_obs) X = np.random.randint(0, 140, num_obs)
X = X[:, None] X = X[:, None]
@ -737,6 +738,7 @@ class GradientTests(np.testing.TestCase):
kern = GPy.kern.Bias(1) + GPy.kern.RBF(1) kern = GPy.kern.Bias(1) + GPy.kern.RBF(1)
lik = GPy.likelihoods.Gaussian() lik = GPy.likelihoods.Gaussian()
m = GPy.models.GPVariationalGaussianApproximation(X, Y, kernel=kern, likelihood=lik) m = GPy.models.GPVariationalGaussianApproximation(X, Y, kernel=kern, likelihood=lik)
m.randomize()
self.assertTrue(m.checkgrad()) self.assertTrue(m.checkgrad())
def test_ssgplvm(self): def test_ssgplvm(self):
@ -744,12 +746,14 @@ class GradientTests(np.testing.TestCase):
from GPy.models import SSGPLVM from GPy.models import SSGPLVM
from GPy.examples.dimensionality_reduction import _simulate_matern from GPy.examples.dimensionality_reduction import _simulate_matern
np.random.seed(10)
D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9 D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9
_, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, False) _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, False)
Y = Ylist[0] Y = Ylist[0]
k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q)
# k = kern.RBF(Q, ARD=True, lengthscale=10.) # k = kern.RBF(Q, ARD=True, lengthscale=10.)
m = SSGPLVM(Y, Q, init="rand", num_inducing=num_inducing, kernel=k, group_spike=True) m = SSGPLVM(Y, Q, init="rand", num_inducing=num_inducing, kernel=k, group_spike=True)
m.randomize()
self.assertTrue(m.checkgrad()) self.assertTrue(m.checkgrad())
if __name__ == "__main__": if __name__ == "__main__":

View file

@ -89,6 +89,9 @@ def _image_directories():
cbook.mkdirs(result_dir) cbook.mkdirs(result_dir)
return baseline_dir, result_dir return baseline_dir, result_dir
baseline_dir, result_dir = _image_directories()
if not os.path.exists(baseline_dir):
raise SkipTest("Not installed from source, baseline not available. Install from source to test plotting")
def _sequenceEqual(a, b): def _sequenceEqual(a, b):
assert len(a) == len(b), "Sequences not same length" assert len(a) == len(b), "Sequences not same length"
@ -99,7 +102,6 @@ def _notFound(path):
raise IOError('File {} not in baseline') raise IOError('File {} not in baseline')
def _image_comparison(baseline_images, extensions=['pdf','svg','png'], tol=11): def _image_comparison(baseline_images, extensions=['pdf','svg','png'], tol=11):
baseline_dir, result_dir = _image_directories()
for num, base in zip(plt.get_fignums(), baseline_images): for num, base in zip(plt.get_fignums(), baseline_images):
for ext in extensions: for ext in extensions:
fig = plt.figure(num) fig = plt.figure(num)

View file

@ -17,5 +17,5 @@ recursive-include GPy *.h
recursive-include GPy *.pyx recursive-include GPy *.pyx
# Testing # Testing
include GPy/testing/baseline/*.png #include GPy/testing/baseline/*.png
#include GPy/testing/pickle_test.pickle #include GPy/testing/pickle_test.pickle

View file

@ -7,7 +7,9 @@ The Gaussian processes framework in Python.
* User [mailing-list](https://lists.shef.ac.uk/sympa/subscribe/gpy-users) * User [mailing-list](https://lists.shef.ac.uk/sympa/subscribe/gpy-users)
* Developer [documentation](http://gpy.readthedocs.org/en/devel/) * Developer [documentation](http://gpy.readthedocs.org/en/devel/)
* Travis-CI [unit-tests](https://travis-ci.org/SheffieldML/GPy) * Travis-CI [unit-tests](https://travis-ci.org/SheffieldML/GPy)
* [![licence](https://img.shields.io/badge/licence-BSD-blue.svg)](http://opensource.org/licenses/BSD-3-Clause) * [![licence](https://img.shields.io/badge/licence-BSD-blue.svg)](http://opensource.org/licenses/BSD-3-Clause)
[![develstat](https://travis-ci.org/SheffieldML/GPy.svg?branch=devel)](https://travis-ci.org/SheffieldML/GPy) [![covdevel](http://codecov.io/github/SheffieldML/GPy/coverage.svg?branch=devel)](http://codecov.io/github/SheffieldML/GPy?branch=devel) [![docdevel](https://readthedocs.org/projects/gpy/badge/?version=devel)](http://gpy.readthedocs.org/en/devel/) [![Research software impact](http://depsy.org/api/package/pypi/GPy/badge.svg)](http://depsy.org/package/python/GPy)
## Updated Structure ## Updated Structure
@ -27,20 +29,14 @@ A warning: This usually works, but sometimes `distutils/setuptools` opens a
whole can of worms here, specially when compiled extensions are involved. whole can of worms here, specially when compiled extensions are involved.
If that is the case, it is best to clean the repo and reinstall. If that is the case, it is best to clean the repo and reinstall.
## Continuous integration
| | Travis-CI | Codecov | RTFD |
| ---: | :--: | :---: | :---: |
| **devel:** | [![develstat](https://travis-ci.org/SheffieldML/GPy.svg?branch=devel)](https://travis-ci.org/SheffieldML/GPy) | [![covdevel](http://codecov.io/github/SheffieldML/GPy/coverage.svg?branch=devel)](http://codecov.io/github/SheffieldML/GPy?branch=devel) | [![docdevel](https://readthedocs.org/projects/gpy/badge/?version=devel)](http://gpy.readthedocs.org/en/devel/) |
## Supported Platforms: ## Supported Platforms:
[<img src="https://www.python.org/static/community_logos/python-logo-generic.svg" height=40px>](https://www.python.org/) [<img src="https://www.python.org/static/community_logos/python-logo-generic.svg" height=40px>](https://www.python.org/)
[<img src="https://upload.wikimedia.org/wikipedia/commons/5/5f/Windows_logo_-_2012.svg" height=40px>](http://www.microsoft.com/en-gb/windows) [<img src="https://upload.wikimedia.org/wikipedia/commons/5/5f/Windows_logo_-_2012.svg" height=40px>](http://www.microsoft.com/en-gb/windows)
[<img src="https://upload.wikimedia.org/wikipedia/commons/8/8e/OS_X-Logo.svg" height=40px>](http://www.apple.com/osx/) [<img src="https://upload.wikimedia.org/wikipedia/commons/8/8e/OS_X-Logo.svg" height=40px>](http://www.apple.com/osx/)
[<img src="https://upload.wikimedia.org/wikipedia/commons/3/35/Tux.svg" height=40px>](https://en.wikipedia.org/wiki/List_of_Linux_distributions) [<img src="https://upload.wikimedia.org/wikipedia/commons/3/35/Tux.svg" height=40px>](https://en.wikipedia.org/wiki/List_of_Linux_distributions)
Python 2.7, 3.3 and higher Python 2.7, 3.4 and higher
## Citation ## Citation
@ -51,14 +47,14 @@ Python 2.7, 3.3 and higher
year = {2012--2015} year = {2012--2015}
} }
### Pronounciation: ### Pronounciation:
We like to pronounce it 'g-pie'. We like to pronounce it 'g-pie'.
## Getting started: installing with pip ## Getting started: installing with pip
We are now requiring the newest version (0.16) of We are now requiring the newest version (0.16) of
[scipy](http://www.scipy.org/) and thus, we strongly recommend using [scipy](http://www.scipy.org/) and thus, we strongly recommend using
the [anaconda python distribution](http://continuum.io/downloads). the [anaconda python distribution](http://continuum.io/downloads).
With anaconda you can install GPy by the following: With anaconda you can install GPy by the following:
@ -105,7 +101,7 @@ or from within IPython
or using setuptools or using setuptools
python setup.py test python setup.py test
## Ubuntu hackers ## Ubuntu hackers
> Note: Right now the Ubuntu package index does not include scipy 0.16.0, and thus, cannot > Note: Right now the Ubuntu package index does not include scipy 0.16.0, and thus, cannot
@ -146,7 +142,7 @@ The HTML files are then stored in doc/build/html
## Funding Acknowledgements ## Funding Acknowledgements
Current support for the GPy software is coming through the following projects. Current support for the GPy software is coming through the following projects.
* [EU FP7-HEALTH Project Ref 305626](http://radiant-project.eu) "RADIANT: Rapid Development and Distribution of Statistical Tools for High-Throughput Sequencing Data" * [EU FP7-HEALTH Project Ref 305626](http://radiant-project.eu) "RADIANT: Rapid Development and Distribution of Statistical Tools for High-Throughput Sequencing Data"

View file

@ -1,5 +1,5 @@
[bumpversion] [bumpversion]
current_version = 1.0.1 current_version = 1.0.2
tag = False tag = False
commit = True commit = True
@ -11,3 +11,6 @@ universal = 1
[upload_docs] [upload_docs]
upload-dir = doc/build/html upload-dir = doc/build/html
[metadata]
description-file = README.rst

View file

@ -182,6 +182,8 @@ if not os.path.exists(user_file):
if os.path.exists(old_user_file): if os.path.exists(old_user_file):
# Move it to new location: # Move it to new location:
print("GPy: Found old config file, moving to new location {}".format(user_file)) print("GPy: Found old config file, moving to new location {}".format(user_file))
if not os.path.exists(os.path.dirname(user_file)):
os.makedirs(os.path.dirname(user_file))
os.rename(old_user_file, user_file) os.rename(old_user_file, user_file)
else: else:
# No config file exists, save informative stub to user config folder: # No config file exists, save informative stub to user config folder: