diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index 4f4de92e..01a9ea41 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -3,6 +3,9 @@
 
 import numpy as np
 from gp_base import GPBase
+from ..util.linalg import dtrtrs
+from ..inference.latent_function_inference import exact_gaussian_inference, expectation_propagation
+from .. import likelihoods
 
 class GP(GPBase):
     """
@@ -18,14 +21,25 @@ class GP(GPBase):
     .. Note:: Multiple independent outputs are allowed using columns of Y
 
     """
-    def __init__(self, X, likelihood, kernel, normalize_X=False):
-        super(GP, self).__init__(X, likelihood, kernel, normalize_X=normalize_X)
-        self.Posterior = self.inference_method.inference(K, likelihood, self.Y)
+    def __init__(self, X, Y, kernel, likelihood, inference_method=None, name='gp'):
+
+        if inference_method is None:
+            if isinstance(likelihood, likelihoods.Gaussian):
+                inference_method = exact_gaussian_inference.ExactGaussianInference()
+        else:
+            inference_method = expectation_propagation
+            print "defaulting to ", inference_method, "for latent function inference"
+
+        super(GP, self).__init__(X, Y, kernel, likelihood, inference_method, name)
+        self.parameters_changed()
 
     def parameters_changed(self):
         super(GP, self).parameters_changed()
         self.K = self.kern.K(self.X)
-        self.Posterior = self.inference_method.inference(K, likelihood, self.Y)
+        self.posterior = self.inference_method.inference(self.K, self.likelihood, self.Y)
+
+    def dL_dtheta_K(self):
+        return self.kern.dK_dtheta(self.posterior.dL_dK, self.X)
 
     def log_likelihood(self):
         return self.posterior.log_marginal
@@ -41,15 +55,15 @@ class GP(GPBase):
 
         """
         Kx = self.kern.K(_Xnew, self.X, which_parts=which_parts).T
-        LiKx, _ = dptrrs(self.posterior._woodbury_chol, np.asfortranarray(Kx), lower=1)
+        LiKx, _ = dtrtrs(self.posterior._woodbury_chol, np.asfortranarray(Kx), lower=1)
         mu = np.dot(Kx.T, self.posterior._woodbury_vector)
         if full_cov:
             Kxx = self.kern.K(_Xnew, which_parts=which_parts)
             var = Kxx - tdot(LiKx.T)
         else:
             Kxx = self.kern.Kdiag(_Xnew, which_parts=which_parts)
-            var = Kxx - np.sum(LiKx.T*LiKx, 0)
-            var = var.reshape(self.num_data, 1)
+            var = Kxx - np.sum(LiKx*LiKx, 0)
+            var = var.reshape(-1, 1)
         return mu, var
 
     def predict(self, Xnew, which_parts='all', full_cov=False, **likelihood_args):
diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py
index dd56e635..827de434 100644
--- a/GPy/core/gp_base.py
+++ b/GPy/core/gp_base.py
@@ -4,34 +4,36 @@ import warnings
 from .. import kern
 from ..util.plot import gpplot, Tango, x_frame1D, x_frame2D
 from model import Model
-from ..core.parameter import ObservableArray
+from parameter import ObservableArray
+from .. import likelihoods
 
 class GPBase(Model):
     """
     Gaussian process base model for holding shared behaviour between
     sparse_GP and GP models.
     """
-    def __init__(self, X, Y, kernel, normalize_X=False, inference_method=None, name=''):
+    def __init__(self, X, Y, kernel, likelihood, inference_method, name=''):
         super(GPBase, self).__init__(name)
         
         assert X.ndim == 2
         self.X = ObservableArray(X)
         self.num_data, self.input_dim = self.X.shape
 
+        assert Y.ndim == 2
+        self.Y = ObservableArray(Y)
+        assert Y.shape[0] == self.num_data
+        _, self.output_dim = self.Y.shape
+
         assert isinstance(kernel, kern.kern)
         self.kern = kernel
 
         assert isinstance(likelihood, likelihoods.Likelihood)
         self.likelihood = likelihood
 
-        if inference_method is None:
-            self.inference_method = self.likelihood.preferred_inference_method
-            print "defaulting to ", inference_method, "for latent function inference"
-        else:
-            self.inference_method = inference_method
+        self.inference_method = inference_method
 
-        assert self.X.shape[0] == Y.shape[0]
-        self.num_data, self.output_dim = self.Y.shape
+        # reinstate this later TODO
+        normalize_X = False
 
         if normalize_X:
             self._Xoffset = X.mean(0)[None, :]
@@ -41,8 +43,8 @@ class GPBase(Model):
             self._Xoffset = np.zeros((1, self.input_dim))
             self._Xscale = np.ones((1, self.input_dim))
         
-        self.add_parameter(self.kern, gradient=lambda:self.kern.dK_dtheta(self.dL_dK, self.X))
-        self.add_parameter(self.likelihood, gradient=lambda:self.likelihood._gradients(partial=np.diag(self.dL_dK)))
+        self.add_parameter(self.kern, gradient=self.dL_dtheta_K)
+        self.add_parameter(self.likelihood, gradient=lambda:self.posterior.dL_dtheta_lik)
 
     def posterior_samples_f(self,X,size=10,which_parts='all',full_cov=True):
         """
diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py
index 6445988a..a7be7334 100644
--- a/GPy/core/sparse_gp.py
+++ b/GPy/core/sparse_gp.py
@@ -5,7 +5,6 @@ import numpy as np
 import pylab as pb
 from ..util.linalg import mdot, jitchol, tdot, symmetrify, backsub_both_sides, chol_inv, dtrtrs, dpotrs, dpotri
 from scipy import linalg
-from ..likelihoods import Gaussian, EP,EP_Mixed_Noise
 from gp_base import GPBase
 from GPy.core.parameter import Param
 
diff --git a/GPy/core/svigp.py b/GPy/core/svigp.py
index c5ea9c6b..5d47611e 100644
--- a/GPy/core/svigp.py
+++ b/GPy/core/svigp.py
@@ -5,7 +5,6 @@ import numpy as np
 import pylab as pb
 from .. import kern
 from ..util.linalg import pdinv, mdot, tdot, dpotrs, dtrtrs, jitchol, backsub_both_sides
-from ..likelihoods import EP
 from gp_base import GPBase
 from model import Model
 import time
diff --git a/GPy/inference/latent_function_inference/__init__.py b/GPy/inference/latent_function_inference/__init__.py
index b843e7ef..fe9dd819 100644
--- a/GPy/inference/latent_function_inference/__init__.py
+++ b/GPy/inference/latent_function_inference/__init__.py
@@ -22,3 +22,6 @@ use this posterior object for making predictions, optimizing hyper-parameters,
 etc.
 
 """
+
+from exact_gaussian_inference import ExactGaussianInference
+expectation_propagation = 'foo' # TODO
diff --git a/GPy/inference/latent_function_inference/exact_gaussian_inference.py b/GPy/inference/latent_function_inference/exact_gaussian_inference.py
index 80e69b81..762744ed 100644
--- a/GPy/inference/latent_function_inference/exact_gaussian_inference.py
+++ b/GPy/inference/latent_function_inference/exact_gaussian_inference.py
@@ -2,11 +2,12 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 from posterior import Posterior
-from .../util.linalg import pdinv, dpotrs, tdot
+from ...util.linalg import pdinv, dpotrs, tdot
+import numpy as np
 log_2_pi = np.log(2*np.pi)
 
 
-class exact_gaussian_inference(object):
+class ExactGaussianInference(object):
     """
     An object for inference when the likelihood is Gaussian.
 
@@ -17,9 +18,9 @@ class exact_gaussian_inference(object):
 
     """
     def __init__(self):
-        self._YYTfactor_cache = caching.cache()
+        pass#self._YYTfactor_cache = caching.cache()
 
-     def get_YYTfactor(self, Y):
+    def get_YYTfactor(self, Y):
         """
         find a matrix L which satisfies LLT = YYT. 
 
@@ -38,16 +39,16 @@ class exact_gaussian_inference(object):
         """
         YYT_factor = self.get_YYTfactor(Y)
 
-        Wi, LW, LWi, W_logdet = pdinv(K + likelhood.covariance(Y, Y_metadata))
+        Wi, LW, LWi, W_logdet = pdinv(K + likelihood.covariance_matrix(Y, Y_metadata))
 
         alpha, _ = dpotrs(LW, YYT_factor, lower=1)
 
         dL_dK = 0.5 * (tdot(alpha) - Y.shape[1] * Wi)
 
-        log_marginal =  0.5*(-Y.size * log_2_pi - Y.shape[1] * W_logdet - np.sum(alpha * YYT_factor.T)
+        log_marginal =  0.5*(-Y.size * log_2_pi - Y.shape[1] * W_logdet - np.sum(alpha * YYT_factor))
 
-        dL_dtheta_lik = likelihood.dL_dtheta(np.diag(dL_dK))
+        dL_dtheta_lik = likelihood._gradients(np.diag(dL_dK))
 
-        return Posterior(log_marginal, dL_DK, dL_dtheta_lik, LW, alpha, K):
+        return Posterior(log_marginal, dL_dK, dL_dtheta_lik, LW, alpha, K)
 
 
diff --git a/GPy/inference/latent_function_inference/posterior.py b/GPy/inference/latent_function_inference/posterior.py
index 6d7f5504..79db69e5 100644
--- a/GPy/inference/latent_function_inference/posterior.py
+++ b/GPy/inference/latent_function_inference/posterior.py
@@ -6,27 +6,27 @@ import numpy as np
 class Posterior(object):
     """
     An object to represent a Gaussian posterior over latent function values.
-    this may be computed exactly for Gaussian likelihoods, or approximated for
+    This may be computed exactly for Gaussian likelihoods, or approximated for
     non-Gaussian likelihoods. 
 
-    The purpose of this clas is to serve as an interface between the inference
+    The purpose of this class is to serve as an interface between the inference
     schemes and the model classes. 
 
     """
-    def __init__(self, log_marginal, dLM_DK, dLM_dtheta_lik, woodbury_chol, woodbury_mean, K):
+    def __init__(self, log_marginal, dL_dK, dL_dtheta_lik, woodbury_chol, woodbury_vector, K):
         """
         log_marginal: log p(Y|X)
-        DLM_dK: d/dK log p(Y|X)
-        dLM_dtheta_lik : d/dtheta log p(Y|X) (where theta are the parameters of the likelihood)
+        DL_dK: d/dK log p(Y|X)
+        dL_dtheta_lik : d/dtheta log p(Y|X) (where theta are the parameters of the likelihood)
         woodbury_chol : a lower triangular matrix L that satisfies posterior_covariance = K - K L^{-T} L^{-1} K
-        woodbury_mean : a matrix (or vector, as Nx1 matrix) M which satisfies posterior_mean = K M
+        woodbury_vector : a matrix (or vector, as Nx1 matrix) M which satisfies posterior_mean = K M
         K : the proir covariance (required for lazy computation of various quantities)
         """
         self.log_marginal = log_marginal
-        self.dLM_DK = dLM_DK
-        self.dLM_dtheta_lik = _dLM_dtheta_lik
+        self.dL_dK = dL_dK
+        self.dL_dtheta_lik = dL_dtheta_lik
         self._woodbury_chol = woodbury_chol
-        self._woodbury_mean = woodbury_mean
+        self._woodbury_vector = woodbury_vector
         self._K = K
 
         #these are computed lazily below
@@ -37,13 +37,13 @@ class Posterior(object):
     @property
     def mean(self):
         if self._mean is None:
-            self._mean = np.dot(self._K, self._woodbury_mean)
+            self._mean = np.dot(self._K, self._woodbury_vector)
         return self._mean
 
     @property
     def covariance(self):
         if self._covariance is None:
-            LiK, _ = dpotrs
+            LiK, _ = dpotrs(self._woodbury_chol, self._K)
             self._covariance = self._K - tdot(LiK.T)
         return self._covariance
 
diff --git a/GPy/inference/optimization/__init__.py b/GPy/inference/optimization/__init__.py
index e69de29b..1a8f043b 100644
--- a/GPy/inference/optimization/__init__.py
+++ b/GPy/inference/optimization/__init__.py
@@ -0,0 +1,2 @@
+from scg import SCG
+from optimization import *
diff --git a/GPy/likelihoods/__init__.py b/GPy/likelihoods/__init__.py
index 0708c763..59e8fb74 100644
--- a/GPy/likelihoods/__init__.py
+++ b/GPy/likelihoods/__init__.py
@@ -3,4 +3,5 @@ from exponential import Exponential
 from gaussian import Gaussian
 from gamma import Gamma
 from poisson import Poisson
-from student_t import Student_t
+from student_t import StudentT
+from likelihood import Likelihood
diff --git a/GPy/likelihoods/bernoulli.py b/GPy/likelihoods/bernoulli.py
index 47737788..9fe96b9a 100644
--- a/GPy/likelihoods/bernoulli.py
+++ b/GPy/likelihoods/bernoulli.py
@@ -2,10 +2,10 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 import numpy as np
-from scipy import stats,special
+from scipy import stats, special
 import scipy as sp
-from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
-import gp_transformations
+from GPy.util.univariate_Gaussian import std_norm_pdf, std_norm_cdf
+import link_functions
 from likelihood import Likelihood
 
 class Bernoulli(Likelihood):
@@ -16,29 +16,29 @@ class Bernoulli(Likelihood):
         p(y_{i}|\\lambda(f_{i})) = \\lambda(f_{i})^{y_{i}}(1-f_{i})^{1-y_{i}}
 
     .. Note::
-        Y is expected to take values in {-1,1}
+        Y is expected to take values in {-1, 1}
         Probit likelihood usually used
     """
-    def __init__(self,gp_link=None,analytical_mean=False,analytical_variance=False):
-        super(Bernoulli, self).__init__(gp_link,analytical_mean,analytical_variance)
-        if isinstance(gp_link , (gp_transformations.Heaviside, gp_transformations.Probit)):
+    def __init__(self, gp_link=None, analytical_mean=False, analytical_variance=False):
+        super(Bernoulli, self).__init__(gp_link, analytical_mean, analytical_variance)
+        if isinstance(gp_link , (link_functions.Heaviside, link_functions.Probit)):
             self.log_concave = True
 
-    def _preprocess_values(self,Y):
+    def _preprocess_values(self, Y):
         """
         Check if the values of the observations correspond to the values
         assumed by the likelihood function.
 
-        ..Note:: Binary classification algorithm works better with classes {-1,1}
+        ..Note:: Binary classification algorithm works better with classes {-1, 1}
         """
         Y_prep = Y.copy()
         Y1 = Y[Y.flatten()==1].size
         Y2 = Y[Y.flatten()==0].size
-        assert Y1 + Y2 == Y.size, 'Bernoulli likelihood is meant to be used only with outputs in {0,1}.'
+        assert Y1 + Y2 == Y.size, 'Bernoulli likelihood is meant to be used only with outputs in {0, 1}.'
         Y_prep[Y.flatten() == 0] = -1
         return Y_prep
 
-    def _moments_match_analytical(self,data_i,tau_i,v_i):
+    def _moments_match_analytical(self, data_i, tau_i, v_i):
         """
         Moments match of the marginal approximation in EP algorithm
 
@@ -51,15 +51,15 @@ class Bernoulli(Likelihood):
         elif data_i == 0:
             sign = -1
         else:
-            raise ValueError("bad value for Bernouilli observation (0,1)")
-        if isinstance(self.gp_link,gp_transformations.Probit):
+            raise ValueError("bad value for Bernouilli observation (0, 1)")
+        if isinstance(self.gp_link, link_functions.Probit):
             z = sign*v_i/np.sqrt(tau_i**2 + tau_i)
             Z_hat = std_norm_cdf(z)
             phi = std_norm_pdf(z)
             mu_hat = v_i/tau_i + sign*phi/(Z_hat*np.sqrt(tau_i**2 + tau_i))
             sigma2_hat = 1./tau_i - (phi/((tau_i**2+tau_i)*Z_hat))*(z+phi/Z_hat)
 
-        elif isinstance(self.gp_link,gp_transformations.Heaviside):
+        elif isinstance(self.gp_link, link_functions.Heaviside):
             a = sign*v_i/np.sqrt(tau_i)
             Z_hat = std_norm_cdf(a)
             N = std_norm_pdf(a)
@@ -68,24 +68,24 @@ class Bernoulli(Likelihood):
             if np.any(np.isnan([Z_hat, mu_hat, sigma2_hat])):
                 stop
         else:
-            raise ValueError("Exact moment matching not available for link {}".format(self.gp_link.gp_transformations.__name__))
+            raise ValueError("Exact moment matching not available for link {}".format(self.gp_link.__name__))
 
         return Z_hat, mu_hat, sigma2_hat
 
-    def _predictive_mean_analytical(self,mu,variance):
+    def _predictive_mean_analytical(self, mu, variance):
 
-        if isinstance(self.gp_link,gp_transformations.Probit):
+        if isinstance(self.gp_link, link_functions.Probit):
             return stats.norm.cdf(mu/np.sqrt(1+variance))
 
-        elif isinstance(self.gp_link,gp_transformations.Heaviside):
+        elif isinstance(self.gp_link, link_functions.Heaviside):
             return stats.norm.cdf(mu/np.sqrt(variance))
 
         else:
             raise NotImplementedError
 
-    def _predictive_variance_analytical(self,mu,variance, pred_mean):
+    def _predictive_variance_analytical(self, mu, variance, pred_mean):
 
-        if isinstance(self.gp_link,gp_transformations.Heaviside):
+        if isinstance(self.gp_link, link_functions.Heaviside):
             return 0.
         else:
             raise NotImplementedError
@@ -106,7 +106,7 @@ class Bernoulli(Likelihood):
         :rtype: float
 
         .. Note:
-            Each y_i must be in {0,1}
+            Each y_i must be in {0, 1}
         """
         assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
         objective = (link_f**y) * ((1.-link_f)**(1.-y))
@@ -195,13 +195,13 @@ class Bernoulli(Likelihood):
         d3logpdf_dlink3 = 2*(y/(link_f**3) - (1-y)/((1-link_f)**3))
         return d3logpdf_dlink3
 
-    def _mean(self,gp):
+    def _mean(self, gp):
         """
         Mass (or density) function
         """
         return self.gp_link.transf(gp)
 
-    def _variance(self,gp):
+    def _variance(self, gp):
         """
         Mass (or density) function
         """
diff --git a/GPy/likelihoods/exponential.py b/GPy/likelihoods/exponential.py
index 16f3456d..8d2e8cdc 100644
--- a/GPy/likelihoods/exponential.py
+++ b/GPy/likelihoods/exponential.py
@@ -6,10 +6,10 @@ import numpy as np
 from scipy import stats,special
 import scipy as sp
 from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
-import gp_transformations
+import link_functions
 from likelihood import Likelihood
 
-class Exponential(NoiseDistribution):
+class Exponential(Likelihood):
     """
     Expoential likelihood
     Y is expected to take values in {0,1,2,...}
diff --git a/GPy/likelihoods/gamma.py b/GPy/likelihoods/gamma.py
index a28130ac..0ac70a9f 100644
--- a/GPy/likelihoods/gamma.py
+++ b/GPy/likelihoods/gamma.py
@@ -6,10 +6,10 @@ import numpy as np
 from scipy import stats,special
 import scipy as sp
 from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
-import gp_transformations
+import link_functions
 from likelihood import Likelihood
 
-class Gamma(NoiseDistribution):
+class Gamma(Likelihood):
     """
     Gamma likelihood
 
diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py
index 0aa6c42b..e3093125 100644
--- a/GPy/likelihoods/gaussian.py
+++ b/GPy/likelihoods/gaussian.py
@@ -2,13 +2,13 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 import numpy as np
-from scipy import stats,special
-import scipy as sp
-from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
-import gp_transformations
+from scipy import stats, special
+from GPy.util.univariate_Gaussian import std_norm_pdf, std_norm_cdf
+import link_functions
 from likelihood import Likelihood
+from ..core.parameter import Param
 
-class Gaussian(NoiseDistribution):
+class Gaussian(Likelihood):
     """
     Gaussian likelihood
 
@@ -19,40 +19,39 @@ class Gaussian(NoiseDistribution):
     :param N: Number of data points
     :type N: int
     """
-    def __init__(self,gp_link=None,analytical_mean=False,analytical_variance=False,variance=1., D=None, N=None):
-        self.variance = variance
-        self.N = N
-        self._set_params(np.asarray(variance))
-        super(Gaussian, self).__init__(gp_link,analytical_mean,analytical_variance)
-        if isinstance(gp_link , gp_transformations.Identity):
+    def __init__(self, gp_link=None, variance=1., name='Gaussian_noise'):
+        if gp_link is None:
+            gp_link = link_functions.Identity()
+
+        if isinstance(gp_link, link_functions.Identity):
+            analytical_variance = True
+            analytical_mean = True
+        else:
+            analytical_variance = False
+            analytical_mean = False
+
+        super(Gaussian, self).__init__(gp_link, analytical_mean, analytical_variance, name=name)
+
+        self.variance = Param('variance', variance)
+        self.add_parameter(self.variance)
+
+        if isinstance(gp_link , link_functions.Identity):
             self.log_concave = True
 
-    def _get_params(self):
-        return np.array([self.variance])
+    def covariance_matrix(self, Y, Y_metadata=None):
+        return np.eye(Y.shape[0]) * self.variance
 
-    def _get_param_names(self):
-        return ['noise_model_variance']
+    def _gradients(self, partial):
+        return np.sum(partial)
 
-    def _set_params(self, p):
-        self.variance = float(p)
-        self.I = np.eye(self.N)
-        self.covariance_matrix = self.I * self.variance
-        self.Ki = self.I*(1.0 / self.variance)
-        #self.ln_det_K = np.sum(np.log(np.diag(self.covariance_matrix)))
-        self.ln_det_K = self.N*np.log(self.variance)
-
-    def _gradients(self,partial):
-        return np.zeros(1)
-        #return np.sum(partial)
-
-    def _preprocess_values(self,Y):
+    def _preprocess_values(self, Y):
         """
         Check if the values of the observations correspond to the values
         assumed by the likelihood function.
         """
         return Y
 
-    def _moments_match_analytical(self,data_i,tau_i,v_i):
+    def _moments_match_analytical(self, data_i, tau_i, v_i):
         """
         Moments match of the marginal approximation in EP algorithm
 
@@ -66,36 +65,13 @@ class Gaussian(NoiseDistribution):
         Z_hat = 1./np.sqrt(2.*np.pi*sum_var)*np.exp(-.5*(data_i - v_i/tau_i)**2./sum_var)
         return Z_hat, mu_hat, sigma2_hat
 
-    def _predictive_mean_analytical(self,mu,sigma):
-        new_sigma2 = self.predictive_variance(mu,sigma)
+    def _predictive_mean_analytical(self, mu, sigma):
+        new_sigma2 = self.predictive_variance(mu, sigma)
         return new_sigma2*(mu/sigma**2 + self.gp_link.transf(mu)/self.variance)
 
-    def _predictive_variance_analytical(self,mu,sigma,predictive_mean=None):
+    def _predictive_variance_analytical(self, mu, sigma, predictive_mean=None):
         return 1./(1./self.variance + 1./sigma**2)
 
-    def _mass(self, link_f, y, extra_data=None):
-        NotImplementedError("Deprecated, now doing chain in noise_model.py for link function evaluation\
-                            Please negate your function and use pdf in noise_model.py, if implementing a likelihood\
-                            rederivate the derivative without doing the chain and put in logpdf, dlogpdf_dlink or\
-                            its derivatives")
-    def _nlog_mass(self, link_f, y, extra_data=None):
-        NotImplementedError("Deprecated, now doing chain in noise_model.py for link function evaluation\
-                            Please negate your function and use logpdf in noise_model.py, if implementing a likelihood\
-                            rederivate the derivative without doing the chain and put in logpdf, dlogpdf_dlink or\
-                            its derivatives")
-
-    def _dnlog_mass_dgp(self, link_f, y, extra_data=None):
-        NotImplementedError("Deprecated, now doing chain in noise_model.py for link function evaluation\
-                            Please negate your function and use dlogpdf_df in noise_model.py, if implementing a likelihood\
-                            rederivate the derivative without doing the chain and put in logpdf, dlogpdf_dlink or\
-                            its derivatives")
-
-    def _d2nlog_mass_dgp2(self, link_f, y, extra_data=None):
-        NotImplementedError("Deprecated, now doing chain in noise_model.py for link function evaluation\
-                            Please negate your function and use d2logpdf_df2 in noise_model.py, if implementing a likelihood\
-                            rederivate the derivative without doing the chain and put in logpdf, dlogpdf_dlink or\
-                            its derivatives")
-
     def pdf_link(self, link_f, y, extra_data=None):
         """
         Likelihood function given link(f)
@@ -270,7 +246,7 @@ class Gaussian(NoiseDistribution):
         d2logpdf_dlink2_dvar = self.d2logpdf_dlink2_dvar(f, y, extra_data=extra_data)
         return d2logpdf_dlink2_dvar
 
-    def _mean(self,gp):
+    def _mean(self, gp):
         """
         Expected value of y under the Mass (or density) function p(y|f)
 
@@ -279,7 +255,7 @@ class Gaussian(NoiseDistribution):
         """
         return self.gp_link.transf(gp)
 
-    def _variance(self,gp):
+    def _variance(self, gp):
         """
         Variance of y under the Mass (or density) function p(y|f)
 
diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py
index ffb6e60d..c1f517f0 100644
--- a/GPy/likelihoods/likelihood.py
+++ b/GPy/likelihoods/likelihood.py
@@ -1,19 +1,19 @@
 # Copyright (c) 2012, 2013 Ricardo Andrade
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
-
 import numpy as np
 from scipy import stats,special
 import scipy as sp
 import pylab as pb
-from GPy.util.plot import gpplot
-from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
-import gp_transformations
-from GPy.util.misc import chain_1, chain_2, chain_3
+from ..util.plot import gpplot
+from ..util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
+import link_functions
+from ..util.misc import chain_1, chain_2, chain_3
 from scipy.integrate import quad
 import warnings
+from ..core.parameterized import Parameterized
 
-class Likelihood(object):
+class Likelihood(Parameterized):
     """
     Likelihood base class
 
@@ -21,8 +21,12 @@ class Likelihood(object):
 
     The minimum required funciotnality is... TODO
     """
-    def __init__(self,gp_link,analytical_mean=False,analytical_variance=False):
-        assert isinstance(gp_link,gp_transformations.GPTransformation), "gp_link is not a valid GPTransformation."
+    def __init__(self,gp_link,analytical_mean=False,analytical_variance=False, name='likelihood_base'):
+        """
+        What are analytical_mean, analyitical_variance? TODO
+        """
+        super(Likelihood, self).__init__(name) 
+        assert isinstance(gp_link,link_functions.GPTransformation), "gp_link is not a valid GPTransformation."
         self.gp_link = gp_link
         self.analytical_mean = analytical_mean
         self.analytical_variance = analytical_variance
@@ -39,15 +43,6 @@ class Likelihood(object):
 
         self.log_concave = False
 
-    def _get_params(self):
-        return np.zeros(0)
-
-    def _get_param_names(self):
-        return []
-
-    def _set_params(self,p):
-        pass
-
     def _gradients(self,partial):
         return np.zeros(0)
 
diff --git a/GPy/likelihoods/poisson.py b/GPy/likelihoods/poisson.py
index 73a10570..355516bb 100644
--- a/GPy/likelihoods/poisson.py
+++ b/GPy/likelihoods/poisson.py
@@ -6,7 +6,7 @@ import numpy as np
 from scipy import stats,special
 import scipy as sp
 from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
-import gp_transformations
+import link_functions
 from likelihood import Likelihood
 
 class Poisson(Likelihood):
diff --git a/GPy/likelihoods/student_t.py b/GPy/likelihoods/student_t.py
index 3336235b..587e1b23 100644
--- a/GPy/likelihoods/student_t.py
+++ b/GPy/likelihoods/student_t.py
@@ -4,7 +4,7 @@
 import numpy as np
 from scipy import stats, special
 import scipy as sp
-import gp_transformations
+import link_functions
 from scipy import stats, integrate
 from scipy.special import gammaln, gamma
 from likelihood import Likelihood
diff --git a/GPy/models/gp_regression.py b/GPy/models/gp_regression.py
index cc58fbb1..4ebf2e25 100644
--- a/GPy/models/gp_regression.py
+++ b/GPy/models/gp_regression.py
@@ -16,28 +16,21 @@ class GPRegression(GP):
     :param X: input observations
     :param Y: observed values
     :param kernel: a GPy kernel, defaults to rbf
-    :param normalize_X:  whether to normalize the input data before computing (predictions will be in original scales)
-    :type normalize_X: False|True
-    :param normalize_Y:  whether to normalize the input data before computing (predictions will be in original scales)
-    :type normalize_Y: False|True
 
     .. Note:: Multiple independent outputs are allowed using columns of Y
 
     """
 
-    def __init__(self, X, Y, kernel=None, normalize_X=False, normalize_Y=False):
+    def __init__(self, X, Y, kernel=None):
         if kernel is None:
             kernel = kern.rbf(X.shape[1])
 
-        likelihood = likelihoods.Gaussian(Y, normalize=normalize_Y)
+        likelihood = likelihoods.Gaussian()
 
-        super(GPRegression, self).__init__(X, likelihood, kernel, normalize_X=normalize_X)
-        self.ensure_default_constraints()
+        super(GPRegression, self).__init__(X, Y, kernel, likelihood, name='gp regression')
 
     def getstate(self):
         return GP.getstate(self)
 
     def setstate(self, state):
         return GP.setstate(self, state)
-
-    pass
diff --git a/GPy/util/misc.py b/GPy/util/misc.py
index 5866ecf9..1cb4c182 100644
--- a/GPy/util/misc.py
+++ b/GPy/util/misc.py
@@ -3,6 +3,34 @@
 
 import numpy as np
 from scipy import weave
+from config import *
+
+def chain_1(df_dg, dg_dx):
+    """
+    Generic chaining function for first derivative
+
+    .. math::
+        \\frac{d(f . g)}{dx} = \\frac{df}{dg} \\frac{dg}{dx}
+    """
+    return df_dg * dg_dx
+
+def chain_2(d2f_dg2, dg_dx, df_dg, d2g_dx2):
+    """
+    Generic chaining function for second derivative
+
+    .. math::
+        \\frac{d^{2}(f . g)}{dx^{2}} = \\frac{d^{2}f}{dg^{2}}(\\frac{dg}{dx})^{2} + \\frac{df}{dg}\\frac{d^{2}g}{dx^{2}}
+    """
+    return d2f_dg2*(dg_dx**2) + df_dg*d2g_dx2
+
+def chain_3(d3f_dg3, dg_dx, d2f_dg2, d2g_dx2, df_dg, d3g_dx3):
+    """
+    Generic chaining function for third derivative
+
+    .. math::
+        \\frac{d^{3}(f . g)}{dx^{3}} = \\frac{d^{3}f}{dg^{3}}(\\frac{dg}{dx})^{3} + 3\\frac{d^{2}f}{dg^{2}}\\frac{dg}{dx}\\frac{d^{2}g}{dx^{2}} + \\frac{df}{dg}\\frac{d^{3}g}{dx^{3}}
+    """
+    return d3f_dg3*(dg_dx**3) + 3*d2f_dg2*dg_dx*d2g_dx2 + df_dg*d3g_dx3
 
 def opt_wrapper(m, **kwargs):
     """
@@ -57,11 +85,18 @@ def kmm_init(X, m = 10):
     return X[inducing]
 
 def fast_array_equal(A, B):
+
+
+    if config.getboolean('parallel', 'openmp'):
+        pragma_string = '#pragma omp parallel for private(i, j)'
+    else:
+        pragma_string = ''
+
     code2="""
     int i, j;
     return_val = 1;
 
-    // #pragma omp parallel for private(i, j)
+    %s
     for(i=0;i<N;i++){
        for(j=0;j<D;j++){
           if(A(i, j) != B(i, j)){
@@ -70,13 +105,18 @@ def fast_array_equal(A, B):
           }
        }
     }
-    """
+    """ % pragma_string
+
+    if config.getboolean('parallel', 'openmp'):
+        pragma_string = '#pragma omp parallel for private(i, j, z)'
+    else:
+        pragma_string = ''
 
     code3="""
     int i, j, z;
     return_val = 1;
 
-    // #pragma omp parallel for private(i, j, z)
+    %s
     for(i=0;i<N;i++){
        for(j=0;j<D;j++){
          for(z=0;z<Q;z++){
@@ -87,20 +127,33 @@ def fast_array_equal(A, B):
           }
        }
     }
-    """
+    """ % pragma_string
+
+    if config.getboolean('parallel', 'openmp'):
+        pragma_string = '#include <omp.h>'
+    else:
+        pragma_string = ''
 
     support_code = """
-    // #include <omp.h>
+    %s
     #include <math.h>
-    """
+    """ % pragma_string
 
-    weave_options = {'headers'           : ['<omp.h>'],
-                     'extra_compile_args': ['-fopenmp -O3'],
-                     'extra_link_args'   : ['-lgomp']}
 
+    weave_options_openmp = {'headers'           : ['<omp.h>'],
+                            'extra_compile_args': ['-fopenmp -O3'],
+                            'extra_link_args'   : ['-lgomp'],
+                            'libraries': ['gomp']}
+    weave_options_noopenmp = {'extra_compile_args': ['-O3']}
+
+    if config.getboolean('parallel', 'openmp'):
+        weave_options = weave_options_openmp
+    else:
+        weave_options = weave_options_noopenmp
 
     value = False
 
+
     if (A == None) and (B == None):
         return True
     elif ((A == None) and (B != None)) or ((A != None) and (B == None)):
@@ -110,14 +163,12 @@ def fast_array_equal(A, B):
             N, D = [int(i) for i in A.shape]
             value = weave.inline(code2, support_code=support_code,
                                  arg_names=['A', 'B', 'N', 'D'],
-                                 type_converters=weave.converters.blitz)
-            # libraries=['gomp'], **weave_options)
+                                 type_converters=weave.converters.blitz, **weave_options)
         elif A.ndim == 3:
             N, D, Q = [int(i) for i in A.shape]
             value = weave.inline(code3, support_code=support_code,
                                  arg_names=['A', 'B', 'N', 'D', 'Q'],
-                                 type_converters=weave.converters.blitz)
-            #libraries=['gomp'], **weave_options)
+                                 type_converters=weave.converters.blitz, **weave_options)
         else:
             value = np.array_equal(A,B)