diff --git a/.travis.yml b/.travis.yml
index 51b9ca2b..63fa1c5e 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -16,8 +16,9 @@ addons:
 env:
   - PYTHON_VERSION=2.7
   #- PYTHON_VERSION=3.3
-  - PYTHON_VERSION=3.4
+  #- PYTHON_VERSION=3.4
   - PYTHON_VERSION=3.5
+  - PYTHON_VERSION=3.6
 
 before_install:
 - wget https://github.com/mzwiessele/travis_scripts/raw/master/download_miniconda.sh
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 51dfe03e..6b1a8a65 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,97 @@
 # Changelog
 
+## v1.7.6 (2017-06-19)
+
+### Fix
+
+* Appveyor not uploading to testpypi for now. [mzwiessele]
+
+### Other
+
+* Bump version: 1.7.5 → 1.7.6. [mzwiessele]
+
+
+## v1.7.5 (2017-06-19)
+
+### Fix
+
+* Splitting forecast tests into 3 to circumvent 10 minute stop of travis. [mzwiessele]
+
+### Other
+
+* Bump version: 1.7.4 → 1.7.5. [mzwiessele]
+
+
+## v1.7.4 (2017-06-19)
+
+### Fix
+
+* Paramz version for parallel optimization fix. [mzwiessele]
+
+### Other
+
+* Bump version: 1.7.3 → 1.7.4. [mzwiessele]
+
+
+## v1.7.3 (2017-06-19)
+
+### Fix
+
+* Appveyor build failing. [mzwiessele]
+
+### Other
+
+* Bump version: 1.7.2 → 1.7.3. [mzwiessele]
+
+
+## v1.7.2 (2017-06-17)
+
+### Fix
+
+* Appveyor build python 3.6. [mzwiessele]
+
+### Other
+
+* Bump version: 1.7.1 → 1.7.2. [mzwiessele]
+
+
+## v1.7.1 (2017-06-17)
+
+### Fix
+
+* Appveyor build python 3.6. [mzwiessele]
+
+### Other
+
+* Bump version: 1.7.0 → 1.7.1. [mzwiessele]
+
+
+## v1.7.0 (2017-06-17)
+
+### Fix
+
+* Support for 3.5 and higher now that 3.6 is out. [mzwiessele]
+
+### Other
+
+* Bump version: 1.6.3 → 1.7.0. [mzwiessele]
+
+
+## v1.6.3 (2017-06-17)
+
+### Other
+
+* Bump version: 1.6.2 → 1.6.3. [mzwiessele]
+
+* Merge pull request #504 from rmcantin/devel. [Max Zwiessele]
+
+* Fix python 2-3 compatibility. [Ruben Martinez-Cantin]
+
+* Merge pull request #511 from dirmeier/devel. [Max Zwiessele]
+
+* Added LICENSE file to MANIFEST.in. [dirmeier]
+
+
 ## v1.6.2 (2017-04-12)
 
 ### Fix
diff --git a/GPy/__version__.py b/GPy/__version__.py
index 51bbb3f2..9f0329de 100644
--- a/GPy/__version__.py
+++ b/GPy/__version__.py
@@ -1 +1 @@
-__version__ = "1.6.2"
+__version__ = "1.7.7"
diff --git a/GPy/inference/latent_function_inference/expectation_propagation.py b/GPy/inference/latent_function_inference/expectation_propagation.py
index 194ee6d6..81c020df 100644
--- a/GPy/inference/latent_function_inference/expectation_propagation.py
+++ b/GPy/inference/latent_function_inference/expectation_propagation.py
@@ -6,6 +6,7 @@ from paramz import ObsAr
 from . import ExactGaussianInference, VarDTC
 from ...util import diag
 from .posterior import PosteriorEP as Posterior
+from ...likelihoods import Gaussian
 
 log_2_pi = np.log(2*np.pi)
 
@@ -174,18 +175,18 @@ class EP(EPBase, ExactGaussianInference):
         if self.ep_mode=="nested":
             #Force EP at each step of the optimization
             self._ep_approximation = None
-            post_params, ga_approx, log_Z_tilde = self._ep_approximation = self.expectation_propagation(K, Y, likelihood, Y_metadata)
+            post_params, ga_approx, cav_params, log_Z_tilde = self._ep_approximation = self.expectation_propagation(K, Y, likelihood, Y_metadata)
         elif self.ep_mode=="alternated":
             if getattr(self, '_ep_approximation', None) is None:
                 #if we don't yet have the results of runnign EP, run EP and store the computed factors in self._ep_approximation
-                post_params, ga_approx, log_Z_tilde = self._ep_approximation = self.expectation_propagation(K, Y, likelihood, Y_metadata)
+                post_params, ga_approx, cav_params, log_Z_tilde = self._ep_approximation = self.expectation_propagation(K, Y, likelihood, Y_metadata)
             else:
                 #if we've already run EP, just use the existing approximation stored in self._ep_approximation
-                post_params, ga_approx, log_Z_tilde = self._ep_approximation
+                post_params, ga_approx, cav_params, log_Z_tilde = self._ep_approximation
         else:
             raise ValueError("ep_mode value not valid")
 
-        return self._inference(K, ga_approx, likelihood, Y_metadata=Y_metadata,  Z_tilde=log_Z_tilde)
+        return self._inference(Y, K, ga_approx, cav_params, likelihood, Y_metadata=Y_metadata,  Z_tilde=log_Z_tilde)
 
     def expectation_propagation(self, K, Y, likelihood, Y_metadata):
 
@@ -220,7 +221,7 @@ class EP(EPBase, ExactGaussianInference):
         # This terms cancel with the coreresponding terms in the marginal loglikelihood
         log_Z_tilde = self._log_Z_tilde(marg_moments, ga_approx, cav_params)
                          # - 0.5*np.log(tau_tilde) + 0.5*(v_tilde*v_tilde*1./tau_tilde)
-        return (post_params, ga_approx, log_Z_tilde)
+        return (post_params, ga_approx, cav_params, log_Z_tilde)
 
     def _init_approximations(self, K, num_data):
         #initial values - Gaussian factors
@@ -280,7 +281,7 @@ class EP(EPBase, ExactGaussianInference):
 
         return log_marginal, post_params
 
-    def _inference(self, K, ga_approx, likelihood, Z_tilde, Y_metadata=None):
+    def _inference(self, Y, K, ga_approx, cav_params, likelihood, Z_tilde, Y_metadata=None):
         log_marginal, post_params = self._ep_marginal(K, ga_approx, Z_tilde)
 
         tau_tilde_root = np.sqrt(ga_approx.tau)
@@ -293,8 +294,7 @@ class EP(EPBase, ExactGaussianInference):
         symmetrify(Wi) #(K + Sigma^(\tilde))^(-1)
 
         dL_dK = 0.5 * (tdot(alpha) - Wi)
-        dL_dthetaL = likelihood.exact_inference_gradients(np.diag(dL_dK), Y_metadata)
-
+        dL_dthetaL = likelihood.ep_gradients(Y, cav_params.tau, cav_params.v, np.diag(dL_dK), Y_metadata=Y_metadata, quad_mode='gh')
         return Posterior(woodbury_inv=Wi, woodbury_vector=alpha, K=K), log_marginal, {'dL_dK':dL_dK, 'dL_dthetaL':dL_dthetaL, 'dL_dm':alpha}
 
 
diff --git a/GPy/likelihoods/binomial.py b/GPy/likelihoods/binomial.py
index e63c009a..61440ec9 100644
--- a/GPy/likelihoods/binomial.py
+++ b/GPy/likelihoods/binomial.py
@@ -66,7 +66,14 @@ class Binomial(Likelihood):
         np.testing.assert_array_equal(N.shape, y.shape)
 
         nchoosey = special.gammaln(N+1) - special.gammaln(y+1) - special.gammaln(N-y+1)
-        return nchoosey + y*np.log(inv_link_f) + (N-y)*np.log(1.-inv_link_f)
+        
+        Ny = N-y
+        t1 = np.zeros(y.shape)
+        t2 = np.zeros(y.shape)
+        t1[y>0] = y[y>0]*np.log(inv_link_f[y>0])
+        t2[Ny>0] = Ny[Ny>0]*np.log(1.-inv_link_f[Ny>0])
+        
+        return nchoosey + t1 + t2
 
     def dlogpdf_dlink(self, inv_link_f, y, Y_metadata=None):
         """
@@ -86,7 +93,13 @@ class Binomial(Likelihood):
         N = Y_metadata['trials']
         np.testing.assert_array_equal(N.shape, y.shape)
 
-        return y/inv_link_f - (N-y)/(1.-inv_link_f)
+        Ny = N-y
+        t1 = np.zeros(y.shape)
+        t2 = np.zeros(y.shape)
+        t1[y>0] = y[y>0]/inv_link_f[y>0]
+        t2[Ny>0] = (Ny[Ny>0])/(1.-inv_link_f[Ny>0])        
+
+        return t1 - t2
 
     def d2logpdf_dlink2(self, inv_link_f, y, Y_metadata=None):
         """
@@ -111,7 +124,13 @@ class Binomial(Likelihood):
         """
         N = Y_metadata['trials']
         np.testing.assert_array_equal(N.shape, y.shape)
-        return -y/np.square(inv_link_f) - (N-y)/np.square(1.-inv_link_f)
+        Ny = N-y
+        t1 = np.zeros(y.shape)
+        t2 = np.zeros(y.shape)
+        t1[y>0] = -y[y>0]/np.square(inv_link_f[y>0])
+        t2[Ny>0] = -(Ny[Ny>0])/np.square(1.-inv_link_f[Ny>0])
+        return t1+t2
+
 
     def d3logpdf_dlink3(self, inv_link_f, y, Y_metadata=None):
         """
@@ -135,8 +154,14 @@ class Binomial(Likelihood):
         N = Y_metadata['trials']
         np.testing.assert_array_equal(N.shape, y.shape)
 
-        inv_link_f2 = np.square(inv_link_f)
-        return 2*y/inv_link_f**3 - 2*(N-y)/(1.-inv_link_f)**3
+        #inv_link_f2 = np.square(inv_link_f)  #TODO Remove. Why is this here?
+        
+        Ny = N-y
+        t1 = np.zeros(y.shape)
+        t2 = np.zeros(y.shape)
+        t1[y>0] = 2*y[y>0]/inv_link_f[y>0]**3
+        t2[Ny>0] = - 2*(Ny[Ny>0])/(1.-inv_link_f[Ny>0])**3
+        return t1 + t2
 
     def samples(self, gp, Y_metadata=None, **kw):
         """
diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py
index 533c6558..04fd6a33 100644
--- a/GPy/likelihoods/gaussian.py
+++ b/GPy/likelihoods/gaussian.py
@@ -57,7 +57,10 @@ class Gaussian(Likelihood):
     def update_gradients(self, grad):
         self.variance.gradient = grad
 
-    def exact_inference_gradients(self, dL_dKdiag,Y_metadata=None):
+    def ep_gradients(self, Y, cav_tau, cav_v, dL_dKdiag, Y_metadata=None, quad_mode='gk', boost_grad=1.):
+        return self.exact_inference_gradients(dL_dKdiag)
+
+    def exact_inference_gradients(self, dL_dKdiag, Y_metadata=None):
         return dL_dKdiag.sum()
 
     def _preprocess_values(self, Y):
diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py
index c5b2094f..308c6a76 100644
--- a/GPy/likelihoods/likelihood.py
+++ b/GPy/likelihoods/likelihood.py
@@ -6,8 +6,12 @@ from scipy import stats,special
 import scipy as sp
 from . import link_functions
 from ..util.misc import chain_1, chain_2, chain_3, blockify_dhess_dtheta, blockify_third, blockify_hessian, safe_exp
+from ..util.quad_integrate import quadgk_int
 from scipy.integrate import quad
+from functools import partial
+
 import warnings
+
 from ..core.parameterization import Parameterized
 
 class Likelihood(Parameterized):
@@ -223,6 +227,91 @@ class Likelihood(Parameterized):
             self.__gh_points = np.polynomial.hermite.hermgauss(T)
         return self.__gh_points
 
+    def ep_gradients(self, Y, cav_tau, cav_v, dL_dKdiag, Y_metadata=None, quad_mode='gk', boost_grad=1.):
+        if self.size > 0:
+            shape = Y.shape
+            tau,v,Y = cav_tau.flatten(), cav_v.flatten(),Y.flatten()
+            mu = v/tau
+            sigma2 = 1./tau
+
+            # assert Y.shape == v.shape
+            dlik_dtheta = np.empty((self.size, Y.shape[0]))
+            # for j in range(self.size):
+            Y_metadata_list = []
+            for index in range(len(Y)):
+                Y_metadata_i = {}
+                if Y_metadata is not None:
+                    for key in Y_metadata.keys():
+                        Y_metadata_i[key] = Y_metadata[key][index,:]
+                    Y_metadata_list.append(Y_metadata_i)
+
+            if quad_mode == 'gk':
+                f = partial(self.integrate_gk)
+                quads = zip(*map(f, Y.flatten(), mu.flatten(), np.sqrt(sigma2.flatten()), Y_metadata_list))
+                quads = np.vstack(quads)
+                quads.reshape(self.size, shape[0], shape[1])
+            elif quad_mode == 'gh':
+                f = partial(self.integrate_gh)
+                quads = zip(*map(f, Y.flatten(), mu.flatten(), np.sqrt(sigma2.flatten())))
+                quads = np.hstack(quads)
+                quads = quads.T
+            else:
+                raise Exception("no other quadrature mode available")
+            #     do a gaussian-hermite integration
+            dL_dtheta_avg = boost_grad * np.nanmean(quads, axis=1)
+            dL_dtheta = boost_grad * np.nansum(quads, axis=1)
+            # dL_dtheta = boost_grad * np.nansum(dlik_dtheta, axis=1)
+        else:
+            dL_dtheta = np.zeros(self.num_params)
+        return dL_dtheta
+
+
+    def integrate_gk(self, Y, mu, sigma, Y_metadata_i=None):
+        # gaussian-kronrod integration.
+        fmin = -np.inf
+        fmax = np.inf
+        SQRT_2PI = np.sqrt(2.*np.pi)
+        def generate_integral(f):
+            a = np.exp(self.logpdf_link(f, Y, Y_metadata_i)) * np.exp(-0.5 * np.square((f - mu) / sigma)) / (
+                SQRT_2PI * sigma)
+            fn1 = a * self.dlogpdf_dtheta(f, Y, Y_metadata_i)
+            fn = fn1
+            return fn
+
+        dF_dtheta_i = quadgk_int(generate_integral, fmin=fmin, fmax=fmax)
+        return dF_dtheta_i
+
+    def integrate_gh(self, Y, mu, sigma, Y_metadata_i=None, gh_points=None):
+        # gaussian-hermite quadrature.
+        # "calculate site derivatives E_f{d logp(y_i|f_i)/da} where a is a likelihood parameter
+        # and the expectation is over the exact marginal posterior, which is not gaussian- and is
+        # unnormalised product of the cavity distribution(a Gaussian) and the exact likelihood term.
+        #
+        # calculate the expectation wrt the approximate marginal posterior, which should be approximately the same.
+        # . This term is needed for evaluating the
+        # gradients of the marginal likelihood estimate Z_EP wrt likelihood parameters."
+        # "writing it explicitly "
+        # use them for gaussian-hermite quadrature
+
+        SQRT_2PI = np.sqrt(2.*np.pi)
+        if gh_points is None:
+            gh_x, gh_w = self._gh_points(32)
+        else:
+            gh_x, gh_w = gh_points
+
+        X = gh_x[None,:]*np.sqrt(2.)*sigma + mu
+
+        # Here X is a grid vector of possible fi values, while Y is just a single value which will be broadcasted.
+        a = np.exp(self.logpdf_link(X, Y, Y_metadata_i))
+        a = a.repeat(self.num_params,0)
+        b = self.dlogpdf_dtheta(X, Y, Y_metadata_i)
+        old_shape = b.shape
+        fn = np.array([i*j for i,j in zip(a.flatten(), b.flatten())])
+        fn = fn.reshape(old_shape)
+
+        dF_dtheta_i = np.dot(fn, gh_w)/np.sqrt(np.pi)
+        return dF_dtheta_i
+
     def variational_expectations(self, Y, m, v, gh_points=None, Y_metadata=None):
         """
         Use Gauss-Hermite Quadrature to compute
diff --git a/GPy/testing/ep_likelihood_tests.py b/GPy/testing/ep_likelihood_tests.py
index 2e1072fa..70efe210 100644
--- a/GPy/testing/ep_likelihood_tests.py
+++ b/GPy/testing/ep_likelihood_tests.py
@@ -28,10 +28,10 @@ class TestObservationModels(unittest.TestCase):
 
         self.Y_noisy = self.Y.copy()
         self.Y_verynoisy = self.Y.copy()
-        self.Y_noisy[75:80] += 1.3
+        self.Y_noisy[75] += 1.3
 
-        self.init_var = 0.3
-        self.deg_free = 5.
+        self.init_var = 0.15
+        self.deg_free = 4.
         censored = np.zeros_like(self.Y)
         random_inds = np.random.choice(self.N, int(self.N / 2), replace=True)
         censored[random_inds] = 1
@@ -83,7 +83,7 @@ class TestObservationModels(unittest.TestCase):
         # taking laplace predictions as the ground truth
         probs_mean_lap, probs_var_lap = m1.predict(self.X)
         probs_mean_ep_alt, probs_var_ep_alt = m2.predict(self.X)
-        probs_mean_ep_nested, probs_var_ep_nested = m2.predict(self.X)
+        probs_mean_ep_nested, probs_var_ep_nested = m3.predict(self.X)
 
         # for simple single dimension data , marginal likelihood for laplace and EP approximations should not be so far apart.
         self.assertAlmostEqual(m1.log_likelihood(), m2.log_likelihood(),delta=1)
@@ -107,12 +107,12 @@ class TestObservationModels(unittest.TestCase):
         ep_inf_nested = GPy.inference.latent_function_inference.EP(ep_mode='nested')
         ep_inf_frac = GPy.inference.latent_function_inference.EP(ep_mode='nested', eta=0.7)
 
-        m1 = GPy.core.GP(self.X, self.Y_noisy.copy(), kernel=self.kernel1, likelihood=studentT.copy(), inference_method=laplace_inf)
+        m1 = GPy.core.GP(self.X.copy(), self.Y_noisy.copy(), kernel=self.kernel1.copy(), likelihood=studentT.copy(), inference_method=laplace_inf)
         # optimize
         m1['.*white'].constrain_fixed(1e-5)
         m1.randomize()
 
-        m2 = GPy.core.GP(self.X, self.Y_noisy.copy(), kernel=self.kernel1, likelihood=studentT.copy(), inference_method=ep_inf_alt)
+        m2 = GPy.core.GP(self.X.copy(), self.Y_noisy.copy(), kernel=self.kernel1.copy(), likelihood=studentT.copy(), inference_method=ep_inf_alt)
         m2['.*white'].constrain_fixed(1e-5)
         # m2.constrain_bounded('.*t_scale2', 0.001, 10)
         m2.randomize()
@@ -124,20 +124,22 @@ class TestObservationModels(unittest.TestCase):
 
         optimizer='bfgs'
         m1.optimize(optimizer=optimizer,max_iters=400)
-        m2.optimize(optimizer=optimizer, max_iters=500)
+        m2.optimize(optimizer=optimizer, max_iters=400)
+        # m3.optimize(optimizer=optimizer, max_iters=500)
+
+        self.assertAlmostEqual(m1.log_likelihood(), m2.log_likelihood(),delta=200)
 
-        self.assertAlmostEqual(m1.log_likelihood(), m2.log_likelihood(),delta=10)
         # self.assertAlmostEqual(m1.log_likelihood(), m3.log_likelihood(), 3)
 
         preds_mean_lap, preds_var_lap = m1.predict(self.X)
         preds_mean_alt, preds_var_alt = m2.predict(self.X)
         # preds_mean_nested, preds_var_nested = m3.predict(self.X)
-        rmse_lap = self.rmse(preds_mean_lap, self.Y_noisy)
-        rmse_alt = self.rmse(preds_mean_alt, self.Y_noisy)
+        rmse_lap = self.rmse(preds_mean_lap, self.Y)
+        rmse_alt = self.rmse(preds_mean_alt, self.Y)
         # rmse_nested = self.rmse(preds_mean_nested, self.Y_noisy)
 
-        if rmse_alt > rmse_alt:
-            self.assertAlmostEqual(rmse_lap, rmse_alt, delta=1.)
+        if rmse_alt > rmse_lap:
+            self.assertAlmostEqual(rmse_lap, rmse_alt, delta=1.5)
         # m3.optimize(optimizer=optimizer, max_iters=500)
 
 
diff --git a/GPy/testing/gpy_kernels_state_space_tests.py b/GPy/testing/gpy_kernels_state_space_tests.py
index f39eb9d0..c06093dd 100644
--- a/GPy/testing/gpy_kernels_state_space_tests.py
+++ b/GPy/testing/gpy_kernels_state_space_tests.py
@@ -306,11 +306,7 @@ class StateSpaceKernelsTests(np.testing.TestCase):
                             gp_kernel=gp_kernel,
                             mean_compare_decimal=2, var_compare_decimal=2)
 
-    def test_forecast(self,):
-        """
-        Test time-series forecasting.
-        """
-
+    def test_forecast_regular(self,):
         # Generate data ->
         np.random.seed(339) # seed the random number generator
         #import pdb; pdb.set_trace()
@@ -334,37 +330,102 @@ class StateSpaceKernelsTests(np.testing.TestCase):
 
         #import pdb; pdb.set_trace()
 
-        def get_new_kernels():
-            periodic_kernel = GPy.kern.StdPeriodic(1,active_dims=[0,])
-            gp_kernel = GPy.kern.Linear(1, active_dims=[0,]) + GPy.kern.Bias(1, active_dims=[0,]) + periodic_kernel
-            gp_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
-            gp_kernel.std_periodic.period.constrain_bounded(0.15, 100)
+        periodic_kernel = GPy.kern.StdPeriodic(1,active_dims=[0,])
+        gp_kernel = GPy.kern.Linear(1, active_dims=[0,]) + GPy.kern.Bias(1, active_dims=[0,]) + periodic_kernel
+        gp_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
+        gp_kernel.std_periodic.period.constrain_bounded(0.15, 100)
 
-            periodic_kernel = GPy.kern.sde_StdPeriodic(1,active_dims=[0,])
-            ss_kernel = GPy.kern.sde_Linear(1,X,active_dims=[0,]) + \
-                GPy.kern.sde_Bias(1, active_dims=[0,]) + periodic_kernel
+        periodic_kernel = GPy.kern.sde_StdPeriodic(1,active_dims=[0,])
+        ss_kernel = GPy.kern.sde_Linear(1,X,active_dims=[0,]) + \
+            GPy.kern.sde_Bias(1, active_dims=[0,]) + periodic_kernel
 
-            ss_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
-            ss_kernel.std_periodic.period.constrain_bounded(0.15, 100)
+        ss_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
+        ss_kernel.std_periodic.period.constrain_bounded(0.15, 100)
 
-            return ss_kernel, gp_kernel
-
-        ss_kernel, gp_kernel = get_new_kernels()
         self.run_for_model(X_train, Y_train, ss_kernel, kalman_filter_type = 'regular',
                            use_cython=False, optimize_max_iters=30, check_gradients=True,
                            predict_X=X_test,
                            gp_kernel=gp_kernel,
                            mean_compare_decimal=2, var_compare_decimal=2)
 
+    def test_forecast_svd(self,):
+        # Generate data ->
+        np.random.seed(339) # seed the random number generator
+        #import pdb; pdb.set_trace()
+        (X,Y) = generate_sine_data(x_points=None, sin_period=5.0, sin_ampl=5.0, noise_var=2.0,
+                        plot = False, points_num=100, x_interval = (0, 40), random=True)
+
+        (X1,Y1) = generate_linear_data(x_points=X, tangent=1.0, add_term=20.0, noise_var=0.0,
+                    plot = False, points_num=100, x_interval = (0, 40), random=True)
+
+        Y = Y + Y1
+
+        X_train = X[X <= 20]
+        Y_train = Y[X <= 20]
+        X_test = X[X > 20]
+        Y_test = Y[X > 20]
+
+        X.shape = (X.shape[0],1); Y.shape = (Y.shape[0],1)
+        X_train.shape = (X_train.shape[0],1); Y_train.shape = (Y_train.shape[0],1)
+        X_test.shape = (X_test.shape[0],1); Y_test.shape = (Y_test.shape[0],1)
+        # Generate data <-
+
+        #import pdb; pdb.set_trace()
+
+        periodic_kernel = GPy.kern.StdPeriodic(1,active_dims=[0,])
+        gp_kernel = GPy.kern.Linear(1, active_dims=[0,]) + GPy.kern.Bias(1, active_dims=[0,]) + periodic_kernel
+        gp_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
+        gp_kernel.std_periodic.period.constrain_bounded(0.15, 100)
+
+        periodic_kernel = GPy.kern.sde_StdPeriodic(1,active_dims=[0,])
+        ss_kernel = GPy.kern.sde_Linear(1,X,active_dims=[0,]) + \
+            GPy.kern.sde_Bias(1, active_dims=[0,]) + periodic_kernel
+
+        ss_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
+        ss_kernel.std_periodic.period.constrain_bounded(0.15, 100)
 
-        ss_kernel, gp_kernel = get_new_kernels()
         self.run_for_model(X_train, Y_train, ss_kernel, kalman_filter_type = 'svd',
                            use_cython=False, optimize_max_iters=30, check_gradients=False,
                            predict_X=X_test,
                            gp_kernel=gp_kernel,
                            mean_compare_decimal=2, var_compare_decimal=2)
 
-        ss_kernel, gp_kernel = get_new_kernels()
+    def test_forecast_svd_cython(self,):
+        # Generate data ->
+        np.random.seed(339) # seed the random number generator
+        #import pdb; pdb.set_trace()
+        (X,Y) = generate_sine_data(x_points=None, sin_period=5.0, sin_ampl=5.0, noise_var=2.0,
+                        plot = False, points_num=100, x_interval = (0, 40), random=True)
+
+        (X1,Y1) = generate_linear_data(x_points=X, tangent=1.0, add_term=20.0, noise_var=0.0,
+                    plot = False, points_num=100, x_interval = (0, 40), random=True)
+
+        Y = Y + Y1
+
+        X_train = X[X <= 20]
+        Y_train = Y[X <= 20]
+        X_test = X[X > 20]
+        Y_test = Y[X > 20]
+
+        X.shape = (X.shape[0],1); Y.shape = (Y.shape[0],1)
+        X_train.shape = (X_train.shape[0],1); Y_train.shape = (Y_train.shape[0],1)
+        X_test.shape = (X_test.shape[0],1); Y_test.shape = (Y_test.shape[0],1)
+        # Generate data <-
+
+        #import pdb; pdb.set_trace()
+
+        periodic_kernel = GPy.kern.StdPeriodic(1,active_dims=[0,])
+        gp_kernel = GPy.kern.Linear(1, active_dims=[0,]) + GPy.kern.Bias(1, active_dims=[0,]) + periodic_kernel
+        gp_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
+        gp_kernel.std_periodic.period.constrain_bounded(0.15, 100)
+
+        periodic_kernel = GPy.kern.sde_StdPeriodic(1,active_dims=[0,])
+        ss_kernel = GPy.kern.sde_Linear(1,X,active_dims=[0,]) + \
+            GPy.kern.sde_Bias(1, active_dims=[0,]) + periodic_kernel
+
+        ss_kernel.std_periodic.lengthscale.constrain_bounded(0.25, 1000)
+        ss_kernel.std_periodic.period.constrain_bounded(0.15, 100)
+
         self.run_for_model(X_train, Y_train, ss_kernel, kalman_filter_type = 'svd',
                            use_cython=True, optimize_max_iters=30, check_gradients=False,
                            predict_X=X_test,
diff --git a/GPy/testing/inference_tests.py b/GPy/testing/inference_tests.py
index f5abf9b5..816d5488 100644
--- a/GPy/testing/inference_tests.py
+++ b/GPy/testing/inference_tests.py
@@ -64,13 +64,13 @@ class InferenceGPEP(unittest.TestCase):
     def genNoisyData(self):
         np.random.seed(1)
         X = np.random.rand(100,1)
-        self.real_std = 0.2
+        self.real_std = 0.1
         noise = np.random.randn(*X[:, 0].shape)*self.real_std
         Y = (np.sin(X[:, 0]*2*np.pi) + noise)[:, None]
         self.f = np.random.rand(X.shape[0],1)
         Y_extra_noisy = Y.copy()
-        Y_extra_noisy[50:53] += 4.
-        Y_extra_noisy[80:83] -= 2.
+        Y_extra_noisy[50] += 4.
+        # Y_extra_noisy[80:83] -= 2.
         return X, Y, Y_extra_noisy
 
     def test_inference_EP(self):
@@ -85,10 +85,11 @@ class InferenceGPEP(unittest.TestCase):
                         inference_method=inf,
                         likelihood=lik)
         K = self.model.kern.K(X)
-        post_params, ga_approx, log_Z_tilde = self.model.inference_method.expectation_propagation(K, ObsAr(Y), lik, None)
+
+        post_params, ga_approx, cav_params, log_Z_tilde = self.model.inference_method.expectation_propagation(K, ObsAr(Y), lik, None)
 
         mu_tilde = ga_approx.v / ga_approx.tau.astype(float)
-        p, m, d = self.model.inference_method._inference(K, ga_approx, lik, Y_metadata=None,  Z_tilde=log_Z_tilde)
+        p, m, d = self.model.inference_method._inference(Y, K, ga_approx, cav_params, lik, Y_metadata=None,  Z_tilde=log_Z_tilde)
         p0, m0, d0 = super(GPy.inference.latent_function_inference.expectation_propagation.EP, inf).inference(k, X,lik ,mu_tilde[:,None], mean_function=None, variance=1./ga_approx.tau, K=K, Z_tilde=log_Z_tilde + np.sum(- 0.5*np.log(ga_approx.tau) + 0.5*(ga_approx.v*ga_approx.v*1./ga_approx.tau)))
 
         assert (np.sum(np.array([m - m0,
@@ -109,19 +110,19 @@ class InferenceGPEP(unittest.TestCase):
     def test_inference_EP_non_classification(self):
         from paramz import ObsAr
         X, Y, Y_extra_noisy = self.genNoisyData()
-        deg_freedom = 5
-        init_noise_var = 0.4
+        deg_freedom = 5.
+        init_noise_var = 0.08
         lik_studentT = GPy.likelihoods.StudentT(deg_free=deg_freedom, sigma2=init_noise_var)
         # like_gaussian_noise = GPy.likelihoods.MixedNoise()
         k = GPy.kern.RBF(1, variance=2., lengthscale=1.1)
-        ep_inf_alt = GPy.inference.latent_function_inference.expectation_propagation.EP(max_iters=100, delta=0.5)
-        ep_inf_nested = GPy.inference.latent_function_inference.expectation_propagation.EP(ep_mode='nested', max_iters=100, delta=0.5)
+        ep_inf_alt = GPy.inference.latent_function_inference.expectation_propagation.EP(max_iters=4, delta=0.5)
+        # ep_inf_nested = GPy.inference.latent_function_inference.expectation_propagation.EP(ep_mode='nested', max_iters=100, delta=0.5)
         m = GPy.core.GP(X=X,Y=Y_extra_noisy,kernel=k,likelihood=lik_studentT,inference_method=ep_inf_alt)
         K = m.kern.K(X)
-        post_params, ga_approx, log_Z_tilde = m.inference_method.expectation_propagation(K, ObsAr(Y_extra_noisy), lik_studentT, None)
+        post_params, ga_approx, cav_params, log_Z_tilde = m.inference_method.expectation_propagation(K, ObsAr(Y_extra_noisy), lik_studentT, None)
 
         mu_tilde = ga_approx.v / ga_approx.tau.astype(float)
-        p, m, d = m.inference_method._inference(K, ga_approx, lik_studentT, Y_metadata=None,  Z_tilde=log_Z_tilde)
+        p, m, d = m.inference_method._inference(Y_extra_noisy, K, ga_approx, cav_params, lik_studentT, Y_metadata=None,  Z_tilde=log_Z_tilde)
         p0, m0, d0 = super(GPy.inference.latent_function_inference.expectation_propagation.EP, ep_inf_alt).inference(k, X,lik_studentT ,mu_tilde[:,None], mean_function=None, variance=1./ga_approx.tau, K=K, Z_tilde=log_Z_tilde + np.sum(- 0.5*np.log(ga_approx.tau) + 0.5*(ga_approx.v*ga_approx.v*1./ga_approx.tau)))
 
         assert (np.sum(np.array([m - m0,
diff --git a/GPy/testing/quadrature_tests.py b/GPy/testing/quadrature_tests.py
new file mode 100644
index 00000000..e519d87e
--- /dev/null
+++ b/GPy/testing/quadrature_tests.py
@@ -0,0 +1,39 @@
+from __future__ import print_function, division
+import numpy as np
+import GPy
+import warnings
+from  ..util.quad_integrate import quadgk_int, quadvgk
+
+
+
+class QuadTests(np.testing.TestCase):
+    """
+    test file for checking implementation of gaussian-kronrod quadrature.
+    we will take a function which can be integrated analytically and check if quadgk result is similar or not!
+    through this file we can test how numerically accurate quadrature implementation in native numpy or manual code is.
+    """
+    def setUp(self):
+        pass
+
+    def test_infinite_quad(self):
+        def f(x):
+            return np.exp(-0.5*x**2)*np.power(x,np.arange(3)[:,None])
+        quad_int_val = quadgk_int(f)
+        real_val = np.sqrt(np.pi * 2)
+        np.testing.assert_almost_equal(real_val, quad_int_val[0], decimal=7)
+
+    def test_finite_quad(self):
+        def f2(x):
+            return x**2
+        quad_int_val = quadvgk(f2, 1.,2.)
+        real_val = 7/3.
+        np.testing.assert_almost_equal(real_val, quad_int_val, decimal=5)
+
+if __name__ == '__main__':
+    def f(x):
+        return np.exp(-0.5 * x ** 2) * np.power(x, np.arange(3)[:, None])
+
+    quad_int_val = quadgk_int(f)
+    real_val = np.sqrt(np.pi*2)
+    np.testing.assert_almost_equal(real_val, quad_int_val[0], decimal=7)
+    print(quadgk_int(f))
diff --git a/GPy/util/__init__.py b/GPy/util/__init__.py
index 685551fd..4994ddcb 100644
--- a/GPy/util/__init__.py
+++ b/GPy/util/__init__.py
@@ -17,3 +17,4 @@ from . import multioutput
 from . import parallel
 from . import functions
 from . import cluster_with_offset
+from . import quad_integrate
diff --git a/GPy/util/datasets.py b/GPy/util/datasets.py
index 6cad1eed..f8fa8239 100644
--- a/GPy/util/datasets.py
+++ b/GPy/util/datasets.py
@@ -206,7 +206,10 @@ def authorize_download(dataset_name=None):
 
 def download_data(dataset_name=None):
     """Check with the user that the are happy with terms and conditions for the data set, then download it."""
-    import itertools
+    try:
+        from itertools import zip_longest
+    except ImportError:
+        from itertools import izip_longest as zip_longest
 
     dr = data_resources[dataset_name]
     if not authorize_download(dataset_name):
@@ -220,8 +223,8 @@ def download_data(dataset_name=None):
     if 'suffices' in dr: zip_urls += (dr['suffices'], )
     else: zip_urls += ([],)
 
-    for url, files, save_names, suffices in itertools.zip_longest(*zip_urls, fillvalue=[]):
-        for f, save_name, suffix in itertools.zip_longest(files, save_names, suffices, fillvalue=None):
+    for url, files, save_names, suffices in zip_longest(*zip_urls, fillvalue=[]):
+        for f, save_name, suffix in zip_longest(files, save_names, suffices, fillvalue=None):
             download_url(os.path.join(url,f), dataset_name, save_name, suffix=suffix)
 
     return True
diff --git a/GPy/util/quad_integrate.py b/GPy/util/quad_integrate.py
new file mode 100644
index 00000000..f3711ada
--- /dev/null
+++ b/GPy/util/quad_integrate.py
@@ -0,0 +1,119 @@
+"""
+The file for utilities related to integration by quadrature methods
+- will contain implementation for gaussian-kronrod integration.
+
+"""
+import numpy as np
+
+def getSubs(Subs, XK, NK=1):
+    M = (Subs[1, :] - Subs[0, :]) / 2
+    C = (Subs[1, :] + Subs[0, :]) / 2
+    I = XK[:, None] * M + np.ones((NK, 1)) * C
+    # A = [Subs(1,:); I]
+    A = np.vstack((Subs[0, :], I))
+    # B = [I;Subs(2,:)]
+    B = np.vstack((I, Subs[1, :]))
+    # Subs = [reshape(A, 1, []);
+    A = A.flatten()
+    # reshape(B, 1, [])];
+    B = B.flatten()
+    Subs = np.vstack((A,B))
+    # Subs = np.concatenate((A, B), axis=0)
+    return Subs
+
+def quadvgk(feval, fmin, fmax, tol1=1e-5, tol2=1e-5):
+    """
+    numpy implementation makes use of the code here: http://se.mathworks.com/matlabcentral/fileexchange/18801-quadvgk
+    We here use gaussian kronrod integration already used in gpstuff for evaluating one dimensional integrals.
+    This is vectorised quadrature which means that several functions can be evaluated at the same time over a grid of
+    points.
+    :param f:
+    :param fmin:
+    :param fmax:
+    :param difftol:
+    :return:
+    """
+
+    XK = np.array([-0.991455371120813, -0.949107912342759, -0.864864423359769, -0.741531185599394,
+                   -0.586087235467691, -0.405845151377397, -0.207784955007898, 0.,
+                   0.207784955007898, 0.405845151377397, 0.586087235467691,
+                   0.741531185599394, 0.864864423359769, 0.949107912342759, 0.991455371120813])
+    WK = np.array([0.022935322010529, 0.063092092629979, 0.104790010322250, 0.140653259715525,
+                   0.169004726639267, 0.190350578064785, 0.204432940075298, 0.209482141084728,
+                   0.204432940075298, 0.190350578064785, 0.169004726639267,
+                   0.140653259715525, 0.104790010322250, 0.063092092629979, 0.022935322010529])
+     # 7-point Gaussian weightings
+    WG = np.array([0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469,
+        0.381830050505119, 0.279705391489277, 0.129484966168870])
+
+    NK = WK.size
+    G = np.arange(2,NK,2)
+    tol1 = 1e-4
+    tol2 = 1e-4
+    Subs = np.array([[fmin],[fmax]])
+    #  number of functions to evaluate in the feval vector of functions.
+    NF = feval(np.zeros(1)).size
+    Q = np.zeros(NF)
+    neval = 0
+    while Subs.size > 0:
+        Subs = getSubs(Subs,XK)
+        M = (Subs[1,:] - Subs[0,:]) / 2
+        C = (Subs[1,:] + Subs[0,:]) / 2
+        # NM = length(M);
+        NM = M.size
+        # x = reshape(XK * M + ones(NK, 1) * C, 1, []);
+        x = XK[:,None]*M + C
+        x = x.flatten()
+        FV = feval(x)
+        # FV = FV[:,None]
+        Q1 = np.zeros((NF, NM))
+        Q2 = np.zeros((NF, NM))
+
+        # for n=1:NF
+        # F = reshape(FV(n,:), NK, []);
+        # Q1(n,:) = M. * sum((WK * ones(1, NM)). * F);
+        # Q2(n,:) = M. * sum((WG * ones(1, NM)). * F(G,:));
+        # end
+        # for i in range(NF):
+        #     F = FV
+        #     F = F.reshape((NK,-1))
+        #     temp_mat = np.sum(np.multiply(WK[:,None]*np.ones((1,NM)), F),axis=0)
+        #     Q1[i,:] = np.multiply(M, temp_mat)
+        #     temp_mat = np.sum(np.multiply(WG[:,None]*np.ones((1, NM)), F[G-1,:]), axis=0)
+        #     Q2[i,:] = np.multiply(M, temp_mat)
+        # ind = np.where(np.logical_or(np.max(np.abs(Q1 -Q2) / Q1) < tol1, (Subs[1,:] - Subs[0,:]) <= tol2) > 0)[0]
+        # Q = Q + np.sum(Q1[:,ind], axis=1)
+        # np.delete(Subs, ind,axis=1)
+
+        Q1 = np.dot(FV.reshape(NF, NK, NM).swapaxes(2,1),WK)*M
+        Q2 = np.dot(FV.reshape(NF, NK, NM).swapaxes(2,1)[:,:,1::2],WG)*M
+        #ind = np.nonzero(np.logical_or(np.max(np.abs((Q1-Q2)/Q1), 0) < difftol , M < xtol))[0]
+        ind = np.nonzero(np.logical_or(np.max(np.abs((Q1-Q2)), 0) < tol1 , (Subs[1,:] - Subs[0,:])  < tol2))[0]
+        Q = Q + np.sum(Q1[:,ind], axis=1)
+        Subs = np.delete(Subs, ind, axis=1)
+    return Q
+
+def quadgk_int(f, fmin=-np.inf, fmax=np.inf, difftol=0.1):
+    """
+    Integrate f from fmin to fmax,
+    do integration by substitution
+    x = r / (1-r**2)
+    when r goes from -1 to 1 , x goes from -inf to inf.
+    the interval for quadgk function is from -1 to +1, so we transform the space from (-inf,inf) to (-1,1)
+    :param f:
+    :param fmin:
+    :param fmax:
+    :param difftol:
+    :return:
+    """
+    difftol = 1e-4
+    def trans_func(r):
+        r2 = np.square(r)
+        x = r / (1-r2)
+        dx_dr = (1 + r2)/(1-r2)**2
+        return f(x)*dx_dr
+
+    integrand = quadvgk(trans_func, -1., 1., difftol, difftol)
+    return integrand
+
+
diff --git a/MANIFEST.in b/MANIFEST.in
index 8e665256..cf220f31 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -16,6 +16,9 @@ recursive-include GPy *.c
 recursive-include GPy *.h
 recursive-include GPy *.pyx
 
+# LICENSE
+include LICENSE.txt
+
 # Testing
 #include GPy/testing/baseline/*.png
 #include GPy/testing/pickle_test.pickle
diff --git a/README.md b/README.md
index 5a771e1b..ffbf6a34 100644
--- a/README.md
+++ b/README.md
@@ -76,7 +76,7 @@ If that is the case, it is best to clean the repo and reinstall.
 [<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)
 
-Python 2.7, 3.4 and higher
+Python 2.7, 3.5 and higher
 
 ## Citation
 
diff --git a/appveyor.yml b/appveyor.yml
index ba454487..73e13280 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -3,12 +3,14 @@ environment:
     secure: 8/ZjXFwtd1S7ixd7PJOpptupKKEDhm2da/q3unabJ00=
   COVERALLS_REPO_TOKEN:
     secure: d3Luic/ESkGaWnZrvWZTKrzO+xaVwJWaRCEP0F+K/9DQGPSRZsJ/Du5g3s4XF+tS
-  gpy_version: 1.6.2
+  gpy_version: 1.7.7
   matrix:
     - PYTHON_VERSION: 2.7
       MINICONDA: C:\Miniconda-x64
     - PYTHON_VERSION: 3.5
       MINICONDA: C:\Miniconda35-x64
+    - PYTHON_VERSION: 3.6
+      MINICONDA: C:\Miniconda36-x64
 
 #configuration:
 #  - Debug
@@ -62,21 +64,21 @@ deploy_script:
 - echo     test >> %USERPROFILE%\\.pypirc
 - echo[
 - echo [pypi] >> %USERPROFILE%\\.pypirc
-- echo username:maxz >> %USERPROFILE%\\.pypirc
-- echo password:%pip_access% >> %USERPROFILE%\\.pypirc
+- echo username = maxz >> %USERPROFILE%\\.pypirc
+- echo password = %pip_access% >> %USERPROFILE%\\.pypirc
 - echo[
 - echo [test] >> %USERPROFILE%\\.pypirc
-- echo repository:https://test.pypi.org/legacy/ >> %USERPROFILE%\\.pypirc
-- echo username:maxz >> %USERPROFILE%\\.pypirc
-- echo password:%pip_access% >> %USERPROFILE%\\.pypirc
+- echo repository = https://testpypi.python.org/pypi >> %USERPROFILE%\\.pypirc
+- echo username = maxz >> %USERPROFILE%\\.pypirc
+- echo password = %pip_access% >> %USERPROFILE%\\.pypirc
 - ps: >-
-    if ($env:APPVEYOR_REPO_BRANCH -eq 'devel') { 
-        twine upload --skip-existing -r test dist/*
+    If ($env:APPVEYOR_REPO_BRANCH -eq 'devel') { 
+        echo not deploying on devel # twine upload --skip-existing -r test dist/*
     }
-    elseif ($env:APPVEYOR_REPO_BRANCH -eq 'deploy') {  
+    ElseIf ($env:APPVEYOR_REPO_BRANCH -eq 'deploy') {  
         twine upload --skip-existing dist/*
     }
-    else {
+    Else {
         echo not deploying on other branches
     }
  
diff --git a/setup.cfg b/setup.cfg
index a52521d3..15ead644 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 1.6.2
+current_version = 1.7.7
 tag = True
 commit = True
 
diff --git a/setup.py b/setup.py
index 82bb5fc2..55f81762 100644
--- a/setup.py
+++ b/setup.py
@@ -150,7 +150,7 @@ setup(name = 'GPy',
       py_modules = ['GPy.__init__'],
       test_suite = 'GPy.testing',
       setup_requires = ['numpy>=1.7'],
-      install_requires = ['numpy>=1.7', 'scipy>=0.16', 'six', 'paramz>=0.6.9'],
+      install_requires = ['numpy>=1.7', 'scipy>=0.16', 'six', 'paramz>=0.7.4'],
       extras_require = {'docs':['sphinx'],
                         'optional':['mpi4py',
                                     'ipython>=4.0.0',
@@ -169,8 +169,8 @@ setup(name = 'GPy',
                    'Operating System :: Microsoft :: Windows',
                    'Operating System :: POSIX :: Linux',
                    'Programming Language :: Python :: 2.7',
-                   'Programming Language :: Python :: 3.3',
                    'Programming Language :: Python :: 3.5',
+                   'Programming Language :: Python :: 3.6',
                    'Framework :: IPython',
                    'Intended Audience :: Science/Research',
                    'Intended Audience :: Developers',