diff --git a/.gitchangelog.rc b/.gitchangelog.rc
index 58d5e788b..8132cf2a7 100644
--- a/.gitchangelog.rc
+++ b/.gitchangelog.rc
@@ -76,7 +76,7 @@ ignore_regexps = [
 ##
 section_regexps = [
     ('New', [
-    r'^[nN]ew\s*:\s*((dev|use?r|pkg|test|doc)\s*:\s*)?([^\n]*)$',
+	r'^[nN]ew\s*:\s*((dev|use?r|pkg|test|doc)\s*:\s*)?([^\n]*)$',
      ]),
     ('Changes', [
         r'^[cC]hg\s*:\s*((dev|use?r|pkg|test|doc)\s*:\s*)?([^\n]*)$',
@@ -87,7 +87,6 @@ section_regexps = [
 
     ('Other', None ## Match all lines
      ),
-
 ]
 
 
@@ -147,7 +146,7 @@ tag_filter_regexp = r'^v[0-9]+\.[0-9]+(\.[0-9]+)?$'
 ##
 ## This label will be used as the changelog Title of the last set of changes
 ## between last valid tag and HEAD if any.
-unreleased_version_label = "%%__version__%% (unreleased)"
+unreleased_version_label = "Unreleased"
 
 
 ## ``output_engine`` is a callable
@@ -178,7 +177,6 @@ unreleased_version_label = "%%__version__%% (unreleased)"
 ##        Examples:
 ##           - makotemplate("restructuredtext")
 ##
-
 #output_engine = rest_py
 #output_engine = mustache("restructuredtext")
 output_engine = mustache("markdown")
@@ -189,4 +187,4 @@ output_engine = mustache("markdown")
 ##
 ## This option tells git-log whether to include merge commits in the log.
 ## The default is to include them.
-include_merge = True
\ No newline at end of file
+include_merge = True
diff --git a/.travis.yml b/.travis.yml
index 7897b0098..51b9ca2b4 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -28,6 +28,10 @@ before_install:
 install:
 - echo $PATH
 - source install_retry.sh
+- if [[ "$TRAVIS_OS_NAME" == "osx" ]];
+  then
+    conda install --yes pandoc;
+  fi;
 - pip install codecov
 - pip install coveralls
 - pip install pypandoc
diff --git a/CHANGELOG.md b/CHANGELOG.md
index e553ffb60..167b4ddca 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,19 +1,86 @@
 # Changelog
 
-## v1.5.5 (2016-10-03)
+## v1.5.6 (2016-11-07)
+
+### New
+
+* Added ploy basis kernel tests and import. [mzwiessele]
+
+* Gitchangelogrc. [mzwiessele]
+
+### Changes
+
+* Added polynomial basis func kernel. [mzwiessele]
+
+### Fix
+
+* Installation #451. [Max Zwiessele]
+
+* Pandoc install under travis osx. [mzwiessele]
+
+* Pandoc install under travis osx. [mzwiessele]
+
+* Pypi changing to pypi.org. [mzwiessele]
 
 ### Other
 
-* Bump version: 1.5.4 → 1.5.5. [Max Zwiessele]
+* Bump version: 1.5.5 → 1.5.6. [mzwiessele]
 
+* Merge pull request #448 from thangbui/devel. [Max Zwiessele]
 
-## v1.5.4 (2016-10-03)
+  Added pep.py -- Sparse Gaussian processes using Power Expectation Propagation
 
-### New
+* Renamed pep test scripts. [Thang Bui]
 
-* Added deployment pull request instructions for developers. [mzwiessele]
+* Fixed seed in pep test script #448. [Thang Bui]
 
-* Using gitchangelog to keep track of changes and log new features. [mzwiessele]
+* Added tests. [Thang Bui]
+
+* Added pep.py -- Sparse Gaussian processes using Power Expectation Propagation. [Thang Bui]
+
+  This allows interpolation between FITC (EP or alpha = 1), and Titsias's variational (VarDTC, VFE when alpha = 0).
+
+* Merge pull request #452 from SheffieldML/setupreq. [Max Zwiessele]
+
+  fix: Installation #451
+
+* Merge pull request #447 from SheffieldML/polinomial. [Max Zwiessele]
+
+  Polynomial
+
+* Merge branch 'devel' into polinomial. [mzwiessele]
+
+* Merge pull request #449 from SheffieldML/deploy. [Max Zwiessele]
+
+  Deploy
+
+* Update setup.py. [Mike Croucher]
+
+* Merge pull request #446 from SheffieldML/devel. [Max Zwiessele]
+
+  newest patch fixing some issues
+
+* Merge branch 'devel' of github.com:SheffieldML/GPy into devel. [mzwiessele]
+
+* Merge branch 'deploy' into devel. [Max Zwiessele]
+
+* Merge pull request #442 from SheffieldML/devel. [Max Zwiessele]
+
+  New Major for GPy
+
+* Merge pull request #426 from SheffieldML/devel. [Max Zwiessele]
+
+  some fixes from issues and beckdaniels warped gp improvements
+
+
+## v1.5.5 (2016-10-03)
+
+### Other
+
+* Bump version: 1.5.4 → 1.5.5. [Max Zwiessele]
+
+
+## v1.5.4 (2016-10-03)
 
 ### Changes
 
@@ -21,14 +88,10 @@
 
 * Fixed naming in variational priors : [Max Zwiessele]
 
-* Changelog update. [mzwiessele]
-
 ### Fix
 
 * Bug in dataset (in fn download_url) which wrongly interprets the Content-Length meta data, and just takes first character. [Michael T Smith]
 
-* What's new update fix #425 in changelog. [mzwiessele]
-
 ### Other
 
 * Bump version: 1.5.3 → 1.5.4. [Max Zwiessele]
@@ -39,25 +102,14 @@
 
 * Merge branch 'kurtCutajar-devel' into devel. [mzwiessele]
 
-* Bump version: 1.5.2 → 1.5.3. [mzwiessele]
-
-* Merge branch 'devel' into kurtCutajar-devel. [mzwiessele]
-
-* Bump version: 1.5.1 → 1.5.2. [mzwiessele]
-
-* Minor readme changes. [mzwiessele]
-
-* Bump version: 1.5.0 → 1.5.1. [mzwiessele]
-
-* Bump version: 1.4.3 → 1.5.0. [mzwiessele]
 
-* Bump version: 1.4.2 → 1.4.3. [mzwiessele]
+## v1.5.3 (2016-09-06)
 
-* Bump version: 1.4.1 → 1.4.2. [mzwiessele]
+### Other
 
-* Merge branch 'devel' of github.com:SheffieldML/GPy into devel. [mzwiessele]
+* Bump version: 1.5.2 → 1.5.3. [mzwiessele]
 
-* [kern] fix #440. [mzwiessele]
+* Merge branch 'devel' into kurtCutajar-devel. [mzwiessele]
 
 * [doc] cleanup. [mzwiessele]
 
@@ -92,6 +144,63 @@
 * Added core code for GpSSM and GpGrid. [kcutajar]
 
 
+## v1.5.2 (2016-09-06)
+
+### New
+
+* Added deployment pull request instructions for developers. [mzwiessele]
+
+### Other
+
+* Bump version: 1.5.1 → 1.5.2. [mzwiessele]
+
+* Minor readme changes. [mzwiessele]
+
+
+## v1.5.1 (2016-09-06)
+
+### Fix
+
+* What's new update fix #425 in changelog. [mzwiessele]
+
+### Other
+
+* Bump version: 1.5.0 → 1.5.1. [mzwiessele]
+
+
+## v1.5.0 (2016-09-06)
+
+### New
+
+* Using gitchangelog to keep track of changes and log new features. [mzwiessele]
+
+### Other
+
+* Bump version: 1.4.3 → 1.5.0. [mzwiessele]
+
+
+## v1.4.3 (2016-09-06)
+
+### Changes
+
+* Changelog update. [mzwiessele]
+
+### Other
+
+* Bump version: 1.4.2 → 1.4.3. [mzwiessele]
+
+
+## v1.4.2 (2016-09-06)
+
+### Other
+
+* Bump version: 1.4.1 → 1.4.2. [mzwiessele]
+
+* Merge branch 'devel' of github.com:SheffieldML/GPy into devel. [mzwiessele]
+
+* [kern] fix #440. [mzwiessele]
+
+
 ## v1.4.1 (2016-09-06)
 
 ### Other
diff --git a/GPy/__version__.py b/GPy/__version__.py
index 8c118e616..d0113b355 100644
--- a/GPy/__version__.py
+++ b/GPy/__version__.py
@@ -1 +1 @@
-__version__ = "1.5.5"
+__version__ = "1.5.6"
diff --git a/GPy/inference/latent_function_inference/__init__.py b/GPy/inference/latent_function_inference/__init__.py
index 90fbf0f1a..3938a6a48 100644
--- a/GPy/inference/latent_function_inference/__init__.py
+++ b/GPy/inference/latent_function_inference/__init__.py
@@ -67,6 +67,7 @@ def __setstate__(self, state):
 from .expectation_propagation import EP, EPDTC
 from .dtc import DTC
 from .fitc import FITC
+from .pep import PEP
 from .var_dtc_parallel import VarDTC_minibatch
 from .var_gauss import VarGauss
 from .gaussian_grid_inference import GaussianGridInference
diff --git a/GPy/inference/latent_function_inference/pep.py b/GPy/inference/latent_function_inference/pep.py
new file mode 100644
index 000000000..79706292e
--- /dev/null
+++ b/GPy/inference/latent_function_inference/pep.py
@@ -0,0 +1,93 @@
+from .posterior import Posterior
+from ...util.linalg import jitchol, tdot, dtrtrs, dtrtri, pdinv
+from ...util import diag
+import numpy as np
+from . import LatentFunctionInference
+log_2_pi = np.log(2*np.pi)
+
+class PEP(LatentFunctionInference):
+    '''
+    Sparse Gaussian processes using Power-Expectation Propagation
+    for regression: alpha \approx 0 gives VarDTC and alpha = 1 gives FITC
+    
+    Reference: A Unifying Framework for Sparse Gaussian Process Approximation using 
+    Power Expectation Propagation, https://arxiv.org/abs/1605.07066
+    
+    '''
+    const_jitter = 1e-6
+    
+    def __init__(self, alpha):
+        super(PEP, self).__init__()
+        self.alpha = alpha
+
+    def inference(self, kern, X, Z, likelihood, Y, mean_function=None, Y_metadata=None):
+        assert mean_function is None, "inference with a mean function not implemented"
+
+        num_inducing, _ = Z.shape
+        num_data, output_dim = Y.shape
+
+        #make sure the noise is not hetero
+        sigma_n = likelihood.gaussian_variance(Y_metadata)
+        if sigma_n.size >1:
+            raise NotImplementedError("no hetero noise with this implementation of PEP")
+
+        Kmm = kern.K(Z)
+        Knn = kern.Kdiag(X)
+        Knm = kern.K(X, Z)
+        U = Knm
+
+        #factor Kmm
+        diag.add(Kmm, self.const_jitter)
+        Kmmi, L, Li, _ = pdinv(Kmm)
+
+        #compute beta_star, the effective noise precision
+        LiUT = np.dot(Li, U.T)
+        sigma_star = sigma_n + self.alpha * (Knn - np.sum(np.square(LiUT),0))
+        beta_star = 1./sigma_star
+
+        # Compute and factor A
+        A = tdot(LiUT*np.sqrt(beta_star)) + np.eye(num_inducing)
+        LA = jitchol(A)
+
+        # back substitute to get b, P, v
+        URiy = np.dot(U.T*beta_star,Y)
+        tmp, _ = dtrtrs(L, URiy, lower=1)
+        b, _ = dtrtrs(LA, tmp, lower=1)
+        tmp, _ = dtrtrs(LA, b, lower=1, trans=1)
+        v, _ = dtrtrs(L, tmp, lower=1, trans=1)
+        tmp, _ = dtrtrs(LA, Li, lower=1, trans=0)
+        P = tdot(tmp.T)
+
+        alpha_const_term = (1.0-self.alpha) / self.alpha
+
+        #compute log marginal
+        log_marginal = -0.5*num_data*output_dim*np.log(2*np.pi) + \
+                       -np.sum(np.log(np.diag(LA)))*output_dim + \
+                       0.5*output_dim*(1+alpha_const_term)*np.sum(np.log(beta_star)) + \
+                       -0.5*np.sum(np.square(Y.T*np.sqrt(beta_star))) + \
+                       0.5*np.sum(np.square(b)) + 0.5*alpha_const_term*num_data*np.log(sigma_n)
+        #compute dL_dR
+        Uv = np.dot(U, v)
+        dL_dR = 0.5*(np.sum(U*np.dot(U,P), 1) - (1.0+alpha_const_term)/beta_star + np.sum(np.square(Y), 1) - 2.*np.sum(Uv*Y, 1) \
+            + np.sum(np.square(Uv), 1))*beta_star**2 
+
+        # Compute dL_dKmm
+        vvT_P = tdot(v.reshape(-1,1)) + P
+        dL_dK = 0.5*(Kmmi - vvT_P)
+        KiU = np.dot(Kmmi, U.T)
+        dL_dK += self.alpha * np.dot(KiU*dL_dR, KiU.T)
+
+        # Compute dL_dU
+        vY = np.dot(v.reshape(-1,1),Y.T)
+        dL_dU = vY - np.dot(vvT_P, U.T)
+        dL_dU *= beta_star
+        dL_dU -= self.alpha * 2.*KiU*dL_dR
+
+        dL_dthetaL = likelihood.exact_inference_gradients(dL_dR)
+        dL_dthetaL += 0.5*alpha_const_term*num_data / sigma_n
+        grad_dict = {'dL_dKmm': dL_dK, 'dL_dKdiag':dL_dR * self.alpha, 'dL_dKnm':dL_dU.T, 'dL_dthetaL':dL_dthetaL}
+
+        #construct a posterior object
+        post = Posterior(woodbury_inv=Kmmi-P, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=L)
+
+        return post, log_marginal, grad_dict
diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py
index 164b66b7d..d82399109 100644
--- a/GPy/kern/__init__.py
+++ b/GPy/kern/__init__.py
@@ -32,7 +32,7 @@
 from .src.splitKern import SplitKern,DEtime
 from .src.splitKern import DEtime as DiffGenomeKern
 from .src.spline import Spline
-from .src.basis_funcs import LogisticBasisFuncKernel, LinearSlopeBasisFuncKernel, BasisFuncKernel, ChangePointBasisFuncKernel, DomainKernel
+from .src.basis_funcs import LogisticBasisFuncKernel, LinearSlopeBasisFuncKernel, BasisFuncKernel, ChangePointBasisFuncKernel, DomainKernel, PolynomialBasisFuncKernel
 from .src.grid_kerns import GridRBF
 
 from .src.sde_matern import sde_Matern32
diff --git a/GPy/kern/src/basis_funcs.py b/GPy/kern/src/basis_funcs.py
index 5d589aa61..569a12f19 100644
--- a/GPy/kern/src/basis_funcs.py
+++ b/GPy/kern/src/basis_funcs.py
@@ -102,6 +102,26 @@ def _K(self, X, X2):
                 phi2 = phi2[:, None]
             return phi1.dot(phi2.T)
 
+class PolynomialBasisFuncKernel(BasisFuncKernel):
+    def __init__(self, input_dim, degree, variance=1., active_dims=None, ARD=True, name='polynomial_basis'):
+        """
+        A linear segment transformation. The segments start at start, \
+        are then linear to stop and constant again. The segments are
+        normalized, so that they have exactly as much mass above
+        as below the origin.
+
+        Start and stop can be tuples or lists of starts and stops.
+        Behaviour of start stop is as np.where(X<start) would do.
+        """
+        self.degree = degree
+        super(PolynomialBasisFuncKernel, self).__init__(input_dim, variance, active_dims, ARD, name)
+
+    @Cache_this(limit=3, ignore_args=())
+    def _phi(self, X):
+        phi = np.empty((X.shape[0], self.degree+1))
+        for i in range(self.degree+1):
+            phi[:, [i]] = X**i
+        return phi
 
 class LinearSlopeBasisFuncKernel(BasisFuncKernel):
     def __init__(self, input_dim, start, stop, variance=1., active_dims=None, ARD=False, name='linear_segment'):
diff --git a/GPy/kern/src/kern.py b/GPy/kern/src/kern.py
index 537d0f7c0..931591fd5 100644
--- a/GPy/kern/src/kern.py
+++ b/GPy/kern/src/kern.py
@@ -3,7 +3,6 @@
 import sys
 import numpy as np
 from ...core.parameterization.parameterized import Parameterized
-from paramz.core.observable_array import ObsAr
 from paramz.caching import Cache_this
 from .kernel_slice_operations import KernCallsViaSlicerMeta
 from functools import reduce
diff --git a/GPy/kern/src/rbf.py b/GPy/kern/src/rbf.py
index 5f3fc97ac..9059ccc03 100644
--- a/GPy/kern/src/rbf.py
+++ b/GPy/kern/src/rbf.py
@@ -42,6 +42,7 @@ def dK2_drdr(self, r):
 
     def dK2_drdr_diag(self):
         return -self.variance # as the diagonal of r is always filled with zeros
+    
     def __getstate__(self):
         dc = super(RBF, self).__getstate__()
         if self.useGPU:
diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py
index 5babeb414..c2725c8e9 100644
--- a/GPy/testing/kernel_tests.py
+++ b/GPy/testing/kernel_tests.py
@@ -495,6 +495,13 @@ def test_basis_func_changepoint(self):
         k = GPy.kern.Add(ks)
         self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
 
+    def test_basis_func_poly(self):
+        ks = []
+        for i in range(self.X.shape[1]):
+            ks.append(GPy.kern.PolynomialBasisFuncKernel(1, 5, ARD=i%2==0, active_dims=[i]))
+        k = GPy.kern.Add(ks)
+        self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
+
     def test_basis_func_domain(self):
         start_stop = np.random.uniform(self.X.min(0), self.X.max(0), (4, self.X.shape[1])).T
         start_stop.sort(axis=1)
diff --git a/GPy/testing/pep_tests.py b/GPy/testing/pep_tests.py
new file mode 100644
index 000000000..2aa6a7848
--- /dev/null
+++ b/GPy/testing/pep_tests.py
@@ -0,0 +1,94 @@
+# Copyright (c) 2014, James Hensman, 2016, Thang Bui
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+import unittest
+import numpy as np
+import GPy
+
+class PEPgradienttest(unittest.TestCase):
+    def setUp(self):
+        ######################################
+        # # 1 dimensional example
+        np.random.seed(10)
+
+        N = 20
+        # sample inputs and outputs
+        self.X1D = np.random.uniform(-3., 3., (N, 1))
+        self.Y1D = np.sin(self.X1D) + np.random.randn(N, 1) * 0.05
+
+        ######################################
+        # # 2 dimensional example
+
+        # sample inputs and outputs
+        self.X2D = np.random.uniform(-3., 3., (N, 2))
+        self.Y2D = np.sin(self.X2D[:, 0:1]) * np.sin(self.X2D[:, 1:2]) + np.random.randn(N, 1) * 0.05
+
+        #######################################
+        # # more datapoints, check in alpha limits, the log marginal likelihood
+        # # is consistent with FITC and VFE/Var_DTC
+        M = 5
+        np.random.seed(42)
+        self.X1 = np.c_[np.linspace(-1., 1., N)]
+        self.Y1 = np.sin(self.X1) + np.random.randn(N, 1) * 0.05
+        self.kernel = GPy.kern.RBF(input_dim=1, lengthscale=0.5, variance=1)
+        self.Z = np.random.uniform(-1, 1, (M, 1))
+        self.lik_noise_var = 0.01
+
+    def test_pep_1d_gradients(self):
+        m = GPy.models.SparseGPRegression(self.X1D, self.Y1D)
+        m.inference_method = GPy.inference.latent_function_inference.PEP(alpha=np.random.rand())
+        self.assertTrue(m.checkgrad())
+
+    def test_pep_2d_gradients(self):
+        m = GPy.models.SparseGPRegression(self.X2D, self.Y2D)
+        m.inference_method = GPy.inference.latent_function_inference.PEP(alpha=np.random.rand())
+        self.assertTrue(m.checkgrad())
+
+    def test_pep_vfe_consistency(self):
+        vfe_model = GPy.models.SparseGPRegression(
+            self.X1, 
+            self.Y1, 
+            kernel=self.kernel, 
+            Z=self.Z
+        )
+        vfe_model.inference_method = GPy.inference.latent_function_inference.VarDTC()
+        vfe_model.Gaussian_noise.variance = self.lik_noise_var
+        vfe_lml = vfe_model.log_likelihood()
+
+        pep_model = GPy.models.SparseGPRegression(
+            self.X1, 
+            self.Y1, 
+            kernel=self.kernel, 
+            Z=self.Z
+        )
+        pep_model.inference_method = GPy.inference.latent_function_inference.PEP(alpha=1e-5)
+        pep_model.Gaussian_noise.variance = self.lik_noise_var
+        pep_lml = pep_model.log_likelihood()
+
+        self.assertAlmostEqual(vfe_lml[0, 0], pep_lml[0], delta=abs(0.01*pep_lml[0]))
+
+    def test_pep_fitc_consistency(self):
+        fitc_model = GPy.models.SparseGPRegression(
+            self.X1D, 
+            self.Y1D, 
+            kernel=self.kernel, 
+            Z=self.Z
+        )
+        fitc_model.inference_method = GPy.inference.latent_function_inference.FITC()
+        fitc_model.Gaussian_noise.variance = self.lik_noise_var
+        fitc_lml = fitc_model.log_likelihood()
+
+        pep_model = GPy.models.SparseGPRegression(
+            self.X1D, 
+            self.Y1D, 
+            kernel=self.kernel, 
+            Z=self.Z
+        )
+        pep_model.inference_method = GPy.inference.latent_function_inference.PEP(alpha=1)
+        pep_model.Gaussian_noise.variance = self.lik_noise_var
+        pep_lml = pep_model.log_likelihood()
+
+        self.assertAlmostEqual(fitc_lml, pep_lml[0], delta=abs(0.001*pep_lml[0]))
+
+
+
diff --git a/appveyor.yml b/appveyor.yml
index 9e7351c17..3738ea1d6 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -3,7 +3,7 @@ environment:
     secure: 8/ZjXFwtd1S7ixd7PJOpptupKKEDhm2da/q3unabJ00=
   COVERALLS_REPO_TOKEN:
     secure: d3Luic/ESkGaWnZrvWZTKrzO+xaVwJWaRCEP0F+K/9DQGPSRZsJ/Du5g3s4XF+tS
-  gpy_version: 1.5.5
+  gpy_version: 1.5.6
   matrix:
     - PYTHON_VERSION: 2.7
       MINICONDA: C:\Miniconda-x64
@@ -66,7 +66,7 @@ deploy_script:
 - echo password:%pip_access% >> %USERPROFILE%\\.pypirc
 - echo[
 - echo [test] >> %USERPROFILE%\\.pypirc
-- echo repository:https://testpypi.python.org/pypi >> %USERPROFILE%\\.pypirc
+- echo repository:https://test.pypi.org/legacy/ >> %USERPROFILE%\\.pypirc
 - echo username:maxz >> %USERPROFILE%\\.pypirc
 - echo password:%pip_access% >> %USERPROFILE%\\.pypirc
 - ps: >-
diff --git a/setup.cfg b/setup.cfg
index 11dee055e..77a46f72e 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 1.5.5
+current_version = 1.5.6
 tag = True
 commit = True
 
diff --git a/setup.py b/setup.py
index 2f1dc25f9..ec18c3387 100644
--- a/setup.py
+++ b/setup.py
@@ -149,6 +149,7 @@ def ismac():
       include_package_data = True,
       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'],
       extras_require = {'docs':['sphinx'],
                         'optional':['mpi4py',