From 7a0b3946e4d57bcc319ac62731b15e298398b66d Mon Sep 17 00:00:00 2001
From: Alexander Grigorievskiy <alex.grigorievskiy@gmail.com>
Date: Tue, 10 Mar 2015 13:10:57 +0200
Subject: [PATCH 1/3] ENH: Adding standard periodic kernel

The standard periodic kernel is the one which is mentioned
in the Rasmussen and Williams book about Gaussian Processes.
---
 GPy/kern/__init__.py               |   1 +
 GPy/kern/_src/standard_periodic.py | 166 +++++++++++++++++++++++++++++
 2 files changed, 167 insertions(+)
 create mode 100644 GPy/kern/_src/standard_periodic.py

diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py
index f479c387..ac10f498 100644
--- a/GPy/kern/__init__.py
+++ b/GPy/kern/__init__.py
@@ -6,6 +6,7 @@ from ._src.brownian import Brownian
 from ._src.stationary import Exponential, OU, Matern32, Matern52, ExpQuad, RatQuad, Cosine
 from ._src.mlp import MLP
 from ._src.periodic import PeriodicExponential, PeriodicMatern32, PeriodicMatern52
+from ._src.standard_periodic import StdPeriodic
 from ._src.independent_outputs import IndependentOutputs, Hierarchical
 from ._src.coregionalize import Coregionalize
 from ._src.ODE_UY import ODE_UY
diff --git a/GPy/kern/_src/standard_periodic.py b/GPy/kern/_src/standard_periodic.py
new file mode 100644
index 00000000..59baa6e3
--- /dev/null
+++ b/GPy/kern/_src/standard_periodic.py
@@ -0,0 +1,166 @@
+# -*- coding: utf-8 -*-
+
+# Copyright (c) 2014, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+"""
+The standard periodic kernel which mentioned in:
+
+[1] Gaussian Processes for Machine Learning, C. E. Rasmussen, C. K. I. Williams.
+The MIT Press, 2005.
+
+
+[2] Introduction to Gaussian processes. D. J. C. MacKay. In C. M. Bishop, editor, 
+Neural Networks and Machine Learning, pages 133-165. Springer, 1998.
+"""
+
+from kern import Kern
+from ...core.parameterization import Param
+from ...core.parameterization.transformations import Logexp
+
+import numpy as np
+
+class StdPeriodic(Kern):
+    """
+    Standart periodic kernel
+
+    .. math::
+
+       k(x,y) = \theta_1 \exp \left[  - \frac{1}{2} {}\sum_{i=1}^{input\_dim}  
+       \left( \frac{\sin(\frac{\pi}{\lambda_i} (x_i - y_i) )}{l_i} \right)^2 \right] }
+
+    :param input_dim: the number of input dimensions
+    :type input_dim: int
+    :param variance: the variance :math:`\theta_1` in the formula above
+    :type variance: float
+    :param wavelength: the vector of wavelengths :math:`\lambda_i`. If None then 1.0 is assumed.
+    :type wavelength: array or list of the appropriate size (or float if there is only one wavelength parameter)
+    :param lengthscale: the vector of lengthscale :math:`\l_i`. If None then 1.0 is assumed.
+    :type lengthscale: array or list of the appropriate size (or float if there is only one lengthscale parameter)
+    :param ARD1: Auto Relevance Determination with respect to wavelength. 
+        If equal to "False" one single wavelength parameter :math:`\lambda_i` for 
+        each dimension is assumed, otherwise there is one lengthscale 
+        parameter per dimension.
+    :type ARD1: Boolean
+    :param ARD2: Auto Relevance Determination with respect to lengthscale. 
+        If equal to "False" one single wavelength parameter :math:`l_i` for 
+        each dimension is assumed, otherwise there is one lengthscale 
+        parameter per dimension.
+    :type ARD2: Boolean
+    :param active_dims: indices of dimensions which are used in the computation of the kernel
+    :type wavelength: array or list of the appropriate size
+    :param name: Name of the kernel for output
+    :type String
+    :param useGPU: whether of not use GPU
+    :type Boolean
+    """
+    
+    def __init__(self, input_dim, variance=1., wavelength=None, lengthscale=None, ARD1=False, ARD2=False, active_dims=None, name='std_periodic',useGPU=False):
+        super(StdPeriodic, self).__init__(input_dim, active_dims, name, useGPU=useGPU)
+        self.input_dim = input_dim
+        self.ARD1 = ARD1 # correspond to wavelengths        
+        self.ARD2 = ARD2 # correspond to lengthscales
+        
+        self.name = name
+        
+        if self.ARD1 == False:
+            if wavelength is not None:
+                wavelength = np.asarray(wavelength)
+                assert wavelength.size == 1, "Only one wavelength needed for non-ARD kernel"
+            else:
+                wavelength = np.ones(1)
+        else:
+            if wavelength is not None:
+                wavelength = np.asarray(wavelength)
+                assert wavelength.size == input_dim, "bad number of wavelengths"
+            else:
+                wavelength = np.ones(input_dim)
+        
+        if self.ARD2 == False:
+            if lengthscale is not None:
+                lengthscale = np.asarray(lengthscale)
+                assert lengthscale.size == 1, "Only one lengthscale needed for non-ARD kernel"
+            else:
+                lengthscale = np.ones(1)
+        else:
+            if lengthscale is not None:
+                lengthscale = np.asarray(lengthscale)
+                assert lengthscale.size == input_dim, "bad number of lengthscales"
+            else:
+                lengthscale = np.ones(input_dim)
+        
+        self.variance = Param('variance', variance, Logexp())
+        assert self.variance.size==1, "Variance size must be one"
+        self.wavelengths =  Param('wavelengths', wavelength, Logexp())
+        self.lengthscales =  Param('lengthscales', lengthscale, Logexp())
+        
+        self.link_parameters(self.variance,  self.wavelengths, self.lengthscales)
+
+    def parameters_changed(self):
+        """
+        This functions deals as a callback for each optimization iteration. 
+        If one optimization step was successfull and the parameters
+        this callback function will be called to be able to update any 
+        precomputations for the kernel.
+        """
+        
+        pass
+        
+        
+    def K(self, X, X2=None):
+        """Compute the covariance matrix between X and X2."""
+        if X2 is None: 
+            X2 = X
+            
+        base = np.pi * (X[:, None, :] - X2[None, :, :]) / self.wavelengths
+        exp_dist = np.exp( -0.5* np.sum( np.square(  np.sin( base ) / self.lengthscales ), axis = -1 ) ) 
+            
+        return self.variance * exp_dist
+
+
+    def Kdiag(self, X):
+        """Compute the diagonal of the covariance matrix associated to X."""
+        ret = np.empty(X.shape[0])
+        ret[:] = self.variance
+        return ret
+
+    def update_gradients_full(self, dL_dK, X, X2=None):
+        """derivative of the covariance matrix with respect to the parameters."""
+        if X2 is None: 
+            X2 = X
+        
+        base = np.pi * (X[:, None, :] - X2[None, :, :]) / self.wavelengths
+        
+        sin_base = np.sin( base )         
+        exp_dist = np.exp( -0.5* np.sum( np.square(  sin_base / self.lengthscales ), axis = -1 ) ) 
+        
+        dwl = self.variance * (1.0/np.square(self.lengthscales)) * sin_base*np.cos(base) * (base / self.wavelengths)
+        
+        dl = self.variance * np.square( sin_base) / np.power( self.lengthscales, 3) 
+        
+        self.variance.gradient = np.sum(exp_dist * dL_dK)    
+        #target[0] += np.sum( exp_dist * dL_dK)        
+        
+        if self.ARD1: # different wavelengths
+            self.wavelengths.gradient = (dwl * exp_dist[:,:,None] * dL_dK[:, :, None]).sum(0).sum(0)
+        else:  # same wavelengths
+            self.wavelengths.gradient = np.sum(dwl.sum(-1) * exp_dist * dL_dK)
+            
+        if self.ARD2: # different lengthscales
+            self.lengthscales.gradient = (dl * exp_dist[:,:,None] * dL_dK[:, :, None]).sum(0).sum(0)
+        else: # same lengthscales
+            self.lengthscales.gradient = np.sum(dl.sum(-1) * exp_dist * dL_dK)
+        
+    def update_gradients_diag(self, dL_dKdiag, X):
+        """derivative of the diagonal of the covariance matrix with respect to the parameters."""
+        self.variance.gradient = np.sum(dL_dKdiag)
+        self.wavelengths.gradient = 0
+        self.lengthscales.gradient = 0
+
+    def gradients_X(self, dL_dK, X, X2=None):
+        """derivative of the covariance matrix with respect to X."""
+    
+        raise NotImplemented("Periodic kernel: dK_dX not implemented")
+
+    def gradients_X_diag(self, dL_dKdiag, X):
+        
+        raise NotImplemented("Periodic kernel: dKdiag_dX not implemented")
\ No newline at end of file

From 0e2f235fc5e591c06ff7ec15debd547c09f3d533 Mon Sep 17 00:00:00 2001
From: Alexander Grigorievskiy <alex.grigorievskiy@gmail.com>
Date: Fri, 17 Apr 2015 18:50:30 +0300
Subject: [PATCH 2/3] TEST: Adding test for periodic kernel

---
 GPy/testing/periodic_kernel_tests.py | 220 +++++++++++++++++++++++++++
 1 file changed, 220 insertions(+)
 create mode 100644 GPy/testing/periodic_kernel_tests.py

diff --git a/GPy/testing/periodic_kernel_tests.py b/GPy/testing/periodic_kernel_tests.py
new file mode 100644
index 00000000..d95330a7
--- /dev/null
+++ b/GPy/testing/periodic_kernel_tests.py
@@ -0,0 +1,220 @@
+# -*- coding: utf-8 -*-
+"""
+Alexander Grigorevskiy, 2015
+"""
+
+import numpy as np
+import scipy as sp
+from scipy import linalg
+import GPy
+import matplotlib.pyplot as plt
+
+def normalize(D):
+    """
+    Function which mormalizes the data to zero mean unit variance.
+    D - column vector, or column-wise matrix
+    """
+    
+    
+    n_rows = D.shape[0]    
+    means = np.nanmean(D, axis= 0)    
+    tmp = D - np.tile( means, (n_rows,1) ) # temporary result. Data with substracted mean                                                                         
+  
+    stds = np.nanstd(tmp,axis=0, ddof=1 ) # one degree of freadom as matlab default
+  
+    result = np.divide( tmp, stds ) 
+    
+    return (result,means,stds)
+    
+def generate_data(func_type_num, gd_plot=False):
+    """
+    Function generates data.
+    """
+    
+    func_types = [ 'linear', 'sin' ]
+    func_type = func_types[func_type_num]   
+
+    
+    time_int = np.array((0, 20)) # time interval
+    points_num = 300
+    
+    t_points = np.random.rand(points_num) * ( time_int[1] - time_int[0] ) + time_int[0]
+    t_points = np.sort( t_points )
+    
+    t_regular_points = np.linspace(time_int[0], time_int[1], num=points_num )    
+    
+    noise_sigma = 1
+    # linear function ->
+    def linear_function(tt):
+        a = 2; b= 1
+        ff = lambda tt:  a*tt + b
+        return ff(tt)
+    # linear function <-
+
+    # sin function ->
+    def sin_function(tt):
+        a = 2; b= 10; # a-period
+        ff = lambda tt:  b * np.sin( 2*np.pi/a * tt )
+        return ff(tt)
+    # sin function <-
+
+    if func_type == 'linear':
+        xx_func = linear_function( t_points ) + np.random.randn( len(t_points) ) * noise_sigma
+        xx_reg_func = linear_function( t_regular_points ) 
+       
+    elif func_type == 'sin':
+        xx_func = sin_function( t_points ) + np.random.randn( len(t_points) ) * noise_sigma
+        xx_reg_func = sin_function( t_regular_points )
+         
+    if gd_plot:
+        plt.figure(1)  
+        plt.plot( t_points, xx_func, 'b.' )
+        plt.plot( t_regular_points, xx_reg_func, '-b' )
+        plt.title('Generated Samples and Underlying Function')
+        plt.show()
+        #plt.close()
+    return  (t_points,  xx_func)
+
+def generate_data2D(func_type_num, gd_plot=False):
+    """
+    Function generates data.
+    """
+    
+    func_type = 'sin'
+
+    
+    time_int = np.array((0, 20)) # time interval
+    points_num = 300
+    
+    t_points_x = np.random.rand(points_num) * ( time_int[1] - time_int[0] ) + time_int[0]
+    t_points_y = np.random.rand(points_num) * ( time_int[1] - time_int[0] ) + time_int[0]
+    
+    t_regular_points_x = np.linspace(time_int[0], time_int[1], num=points_num )    
+    t_regular_points_y = np.linspace(time_int[0], time_int[1], num=points_num )
+    
+    noise_sigma = 1
+   
+    # sin function ->
+    def sin_function(tt1,tt2):
+        a1 = 2; a2=8; b= 10; # a - wavelengths
+        ff = lambda tt1, tt2:  b * np.sin( 2*np.pi/a1 * tt1 + 2*np.pi/a2 * tt2)
+        return ff(tt1,tt2)
+    # sin function <-
+
+   
+      
+    if func_type == 'sin':
+        xx_func = sin_function( t_points_x, t_points_y) + np.random.randn( len(t_points_x) ) * noise_sigma
+        xx_reg_func = sin_function( t_regular_points_x, t_regular_points_y)
+         
+#    if gd_plot:
+#        plt.figure(1)
+#        xx,yy = np.meshgrid(t_regular_points_x, t_regular_points_x)
+#        
+#        plt.title('Generated Samples and Underlying Function')
+#        plt.show()
+#        #plt.close()
+    return  (t_points_x, t_points_y, xx_func)
+    
+def test_simple():
+    """
+    Simple test
+    """
+    
+    (xx, yy)  =  generate_data( 1, gd_plot=True )
+    xx = xx[:, np.newaxis ]; 
+    yy = yy[:, np.newaxis ];  
+    
+    input_dim = 1
+    xx1 = xx # np.hstack( (xx, xx**2))    # possible input data transformation
+    
+    (xx_norm, xx_mean, xx_std) = normalize( xx1) # normalize in place   
+    (yy_norm, yy_mean, yy_std) = normalize( yy ) # normalize in place
+   
+
+    #periodic_kernel = GPy.kern.Linear(input_dim)
+    periodic_kernel = GPy.kern.StdPeriodic(input_dim)
+    
+    KK = periodic_kernel
+ 
+   # Show sample paths ->    
+    Cov = KK.K( xx_norm, xx_norm); # KK.plot();GPy.kern.periodic_Matern52
+    #plt.figure(1)
+    plt.matshow(Cov); plt.colorbar()
+    plt.title('Covariance BEFORE optimization')
+    #plt.show()        
+    
+    zz = np.random.multivariate_normal( np.zeros(xx.shape[0]), Cov, 5 )
+    plt.figure(2)  
+    for i in range(5):
+        plt.plot(xx_norm,zz[i,:])
+    plt.title('Covariance sample paths')
+    plt.show()
+    # Show sample paths <-
+    
+    gpy_reg = GPy.models.GPRegression( xx_norm, yy, KK )
+    print "Before Optimization ->\n"
+    print gpy_reg
+    gpy_reg.plot()
+    
+    gpy_reg.optimize_restarts(num_restarts = 3, robust=True, parallel=False, num_processes=4)    
+    
+    print "After Optimization ->\n"
+    print gpy_reg
+    gpy_reg.plot()
+     
+    return gpy_reg
+    
+def test_2D_sine():
+    """
+    Test function of 2D argument.
+    """
+    
+    (xx, yy,zz)  =  generate_data2D( 1, gd_plot=True )
+    xx = xx[:, np.newaxis ]; 
+    yy = yy[:, np.newaxis ];  
+    zz = zz[:, np.newaxis ];  
+    
+    input_dim = 1
+    xx1 = np.hstack( (xx, yy))    # possible input data transformation
+    
+    (xx_norm, xx_mean, xx_std) = normalize( xx1) # normalize in place   
+    (yy_norm, yy_mean, yy_std) = normalize( yy ) # normalize in place
+   
+
+    #periodic_kernel = GPy.kern.Linear(input_dim)
+    periodic_kernel = GPy.kern.StdPeriodic(input_dim,active_dims=[0,])
+    
+    KK = periodic_kernel
+ 
+   # Show sample paths ->    
+    Cov = KK.K( xx_norm, xx_norm); # KK.plot();GPy.kern.periodic_Matern52
+    #plt.figure(1)
+    plt.matshow(Cov); plt.colorbar()
+    plt.title('Covariance BEFORE optimization')
+    #plt.show()        
+    
+    zz = np.random.multivariate_normal( np.zeros(xx.shape[0]), Cov, 5 )
+    plt.figure(2)  
+    for i in range(5):
+        plt.plot(xx_norm,zz[i,:])
+    plt.title('Covariance sample paths')
+    plt.show()
+    # Show sample paths <-
+    
+    gpy_reg = GPy.models.GPRegression( xx_norm, yy, KK )
+    print "Before Optimization ->\n"
+    print gpy_reg
+    gpy_reg.plot()
+    gpy_reg.optimize_restarts(num_restarts = 3, robust=True, parallel=False, num_processes=4)    
+    
+    print "After Optimization ->\n"
+    print gpy_reg
+    gpy_reg.plot()
+    
+    return gpy_reg
+    
+    
+if __name__ == '__main__':   
+    reg = test_simple()
+    #reg = test_2D_sine()
\ No newline at end of file

From acc15e59d2f41b6a306b410cfe63f2c55ea66df3 Mon Sep 17 00:00:00 2001
From: Alexander Grigorievskiy <alex.grigorievskiy@gmail.com>
Date: Wed, 15 Jul 2015 13:06:27 +0300
Subject: [PATCH 3/3] TEST: Making test for the Standard Periodic Kernel
 similar to other kernel tests.

---
 GPy/kern/_src/standard_periodic.py   |  16 +-
 GPy/testing/kernel_tests.py          |   7 +-
 GPy/testing/periodic_kernel_tests.py | 220 ---------------------------
 3 files changed, 14 insertions(+), 229 deletions(-)
 delete mode 100644 GPy/testing/periodic_kernel_tests.py

diff --git a/GPy/kern/_src/standard_periodic.py b/GPy/kern/_src/standard_periodic.py
index 59baa6e3..8e2966b2 100644
--- a/GPy/kern/_src/standard_periodic.py
+++ b/GPy/kern/_src/standard_periodic.py
@@ -156,11 +156,11 @@ class StdPeriodic(Kern):
         self.wavelengths.gradient = 0
         self.lengthscales.gradient = 0
 
-    def gradients_X(self, dL_dK, X, X2=None):
-        """derivative of the covariance matrix with respect to X."""
-    
-        raise NotImplemented("Periodic kernel: dK_dX not implemented")
-
-    def gradients_X_diag(self, dL_dKdiag, X):
-        
-        raise NotImplemented("Periodic kernel: dKdiag_dX not implemented")
\ No newline at end of file
+#    def gradients_X(self, dL_dK, X, X2=None):
+#        """derivative of the covariance matrix with respect to X."""
+#    
+#        raise NotImplemented("Periodic kernel: dK_dX not implemented")
+#
+#    def gradients_X_diag(self, dL_dKdiag, X):
+#        
+#        raise NotImplemented("Periodic kernel: dKdiag_dX not implemented")
\ No newline at end of file
diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py
index 043d5e9a..e8ce6bb0 100644
--- a/GPy/testing/kernel_tests.py
+++ b/GPy/testing/kernel_tests.py
@@ -312,7 +312,12 @@ class KernelGradientTestsContinuous(unittest.TestCase):
         k = GPy.kern.LinearFull(self.D, self.D-1)
         k.randomize()
         self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
-
+        
+    def test_standard_periodic(self):
+        k = GPy.kern.StdPeriodic(self.D, self.D-1)
+        k.randomize()
+        self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
+        
 class KernelTestsMiscellaneous(unittest.TestCase):
     def setUp(self):
         N, D = 100, 10
diff --git a/GPy/testing/periodic_kernel_tests.py b/GPy/testing/periodic_kernel_tests.py
deleted file mode 100644
index d95330a7..00000000
--- a/GPy/testing/periodic_kernel_tests.py
+++ /dev/null
@@ -1,220 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Alexander Grigorevskiy, 2015
-"""
-
-import numpy as np
-import scipy as sp
-from scipy import linalg
-import GPy
-import matplotlib.pyplot as plt
-
-def normalize(D):
-    """
-    Function which mormalizes the data to zero mean unit variance.
-    D - column vector, or column-wise matrix
-    """
-    
-    
-    n_rows = D.shape[0]    
-    means = np.nanmean(D, axis= 0)    
-    tmp = D - np.tile( means, (n_rows,1) ) # temporary result. Data with substracted mean                                                                         
-  
-    stds = np.nanstd(tmp,axis=0, ddof=1 ) # one degree of freadom as matlab default
-  
-    result = np.divide( tmp, stds ) 
-    
-    return (result,means,stds)
-    
-def generate_data(func_type_num, gd_plot=False):
-    """
-    Function generates data.
-    """
-    
-    func_types = [ 'linear', 'sin' ]
-    func_type = func_types[func_type_num]   
-
-    
-    time_int = np.array((0, 20)) # time interval
-    points_num = 300
-    
-    t_points = np.random.rand(points_num) * ( time_int[1] - time_int[0] ) + time_int[0]
-    t_points = np.sort( t_points )
-    
-    t_regular_points = np.linspace(time_int[0], time_int[1], num=points_num )    
-    
-    noise_sigma = 1
-    # linear function ->
-    def linear_function(tt):
-        a = 2; b= 1
-        ff = lambda tt:  a*tt + b
-        return ff(tt)
-    # linear function <-
-
-    # sin function ->
-    def sin_function(tt):
-        a = 2; b= 10; # a-period
-        ff = lambda tt:  b * np.sin( 2*np.pi/a * tt )
-        return ff(tt)
-    # sin function <-
-
-    if func_type == 'linear':
-        xx_func = linear_function( t_points ) + np.random.randn( len(t_points) ) * noise_sigma
-        xx_reg_func = linear_function( t_regular_points ) 
-       
-    elif func_type == 'sin':
-        xx_func = sin_function( t_points ) + np.random.randn( len(t_points) ) * noise_sigma
-        xx_reg_func = sin_function( t_regular_points )
-         
-    if gd_plot:
-        plt.figure(1)  
-        plt.plot( t_points, xx_func, 'b.' )
-        plt.plot( t_regular_points, xx_reg_func, '-b' )
-        plt.title('Generated Samples and Underlying Function')
-        plt.show()
-        #plt.close()
-    return  (t_points,  xx_func)
-
-def generate_data2D(func_type_num, gd_plot=False):
-    """
-    Function generates data.
-    """
-    
-    func_type = 'sin'
-
-    
-    time_int = np.array((0, 20)) # time interval
-    points_num = 300
-    
-    t_points_x = np.random.rand(points_num) * ( time_int[1] - time_int[0] ) + time_int[0]
-    t_points_y = np.random.rand(points_num) * ( time_int[1] - time_int[0] ) + time_int[0]
-    
-    t_regular_points_x = np.linspace(time_int[0], time_int[1], num=points_num )    
-    t_regular_points_y = np.linspace(time_int[0], time_int[1], num=points_num )
-    
-    noise_sigma = 1
-   
-    # sin function ->
-    def sin_function(tt1,tt2):
-        a1 = 2; a2=8; b= 10; # a - wavelengths
-        ff = lambda tt1, tt2:  b * np.sin( 2*np.pi/a1 * tt1 + 2*np.pi/a2 * tt2)
-        return ff(tt1,tt2)
-    # sin function <-
-
-   
-      
-    if func_type == 'sin':
-        xx_func = sin_function( t_points_x, t_points_y) + np.random.randn( len(t_points_x) ) * noise_sigma
-        xx_reg_func = sin_function( t_regular_points_x, t_regular_points_y)
-         
-#    if gd_plot:
-#        plt.figure(1)
-#        xx,yy = np.meshgrid(t_regular_points_x, t_regular_points_x)
-#        
-#        plt.title('Generated Samples and Underlying Function')
-#        plt.show()
-#        #plt.close()
-    return  (t_points_x, t_points_y, xx_func)
-    
-def test_simple():
-    """
-    Simple test
-    """
-    
-    (xx, yy)  =  generate_data( 1, gd_plot=True )
-    xx = xx[:, np.newaxis ]; 
-    yy = yy[:, np.newaxis ];  
-    
-    input_dim = 1
-    xx1 = xx # np.hstack( (xx, xx**2))    # possible input data transformation
-    
-    (xx_norm, xx_mean, xx_std) = normalize( xx1) # normalize in place   
-    (yy_norm, yy_mean, yy_std) = normalize( yy ) # normalize in place
-   
-
-    #periodic_kernel = GPy.kern.Linear(input_dim)
-    periodic_kernel = GPy.kern.StdPeriodic(input_dim)
-    
-    KK = periodic_kernel
- 
-   # Show sample paths ->    
-    Cov = KK.K( xx_norm, xx_norm); # KK.plot();GPy.kern.periodic_Matern52
-    #plt.figure(1)
-    plt.matshow(Cov); plt.colorbar()
-    plt.title('Covariance BEFORE optimization')
-    #plt.show()        
-    
-    zz = np.random.multivariate_normal( np.zeros(xx.shape[0]), Cov, 5 )
-    plt.figure(2)  
-    for i in range(5):
-        plt.plot(xx_norm,zz[i,:])
-    plt.title('Covariance sample paths')
-    plt.show()
-    # Show sample paths <-
-    
-    gpy_reg = GPy.models.GPRegression( xx_norm, yy, KK )
-    print "Before Optimization ->\n"
-    print gpy_reg
-    gpy_reg.plot()
-    
-    gpy_reg.optimize_restarts(num_restarts = 3, robust=True, parallel=False, num_processes=4)    
-    
-    print "After Optimization ->\n"
-    print gpy_reg
-    gpy_reg.plot()
-     
-    return gpy_reg
-    
-def test_2D_sine():
-    """
-    Test function of 2D argument.
-    """
-    
-    (xx, yy,zz)  =  generate_data2D( 1, gd_plot=True )
-    xx = xx[:, np.newaxis ]; 
-    yy = yy[:, np.newaxis ];  
-    zz = zz[:, np.newaxis ];  
-    
-    input_dim = 1
-    xx1 = np.hstack( (xx, yy))    # possible input data transformation
-    
-    (xx_norm, xx_mean, xx_std) = normalize( xx1) # normalize in place   
-    (yy_norm, yy_mean, yy_std) = normalize( yy ) # normalize in place
-   
-
-    #periodic_kernel = GPy.kern.Linear(input_dim)
-    periodic_kernel = GPy.kern.StdPeriodic(input_dim,active_dims=[0,])
-    
-    KK = periodic_kernel
- 
-   # Show sample paths ->    
-    Cov = KK.K( xx_norm, xx_norm); # KK.plot();GPy.kern.periodic_Matern52
-    #plt.figure(1)
-    plt.matshow(Cov); plt.colorbar()
-    plt.title('Covariance BEFORE optimization')
-    #plt.show()        
-    
-    zz = np.random.multivariate_normal( np.zeros(xx.shape[0]), Cov, 5 )
-    plt.figure(2)  
-    for i in range(5):
-        plt.plot(xx_norm,zz[i,:])
-    plt.title('Covariance sample paths')
-    plt.show()
-    # Show sample paths <-
-    
-    gpy_reg = GPy.models.GPRegression( xx_norm, yy, KK )
-    print "Before Optimization ->\n"
-    print gpy_reg
-    gpy_reg.plot()
-    gpy_reg.optimize_restarts(num_restarts = 3, robust=True, parallel=False, num_processes=4)    
-    
-    print "After Optimization ->\n"
-    print gpy_reg
-    gpy_reg.plot()
-    
-    return gpy_reg
-    
-    
-if __name__ == '__main__':   
-    reg = test_simple()
-    #reg = test_2D_sine()
\ No newline at end of file