function [gp] = init_gp_1d (gp)
% Initialise 1D Gaussian process model
% FORMAT [gp] = init_gp_1d (gp)
% 
% gp.       Data structure with fields
%           N         Number of discrete sample points in [0,1]
%           cfun      Covariance function; 'ard'

try
    gp.N;
catch
    gp.N=100;
end
gp.C=zeros(gp.N,gp.N);
gp.v=[1:gp.N]/gp.N;

% Get covariance function 
switch gp.cfun
    
    case 'ard',
        % Equation 45 in Mackay paper
        try
            gp.theta1;
        catch
            gp.theta1=5;
        end
        try
            gp.theta2;
        catch
            gp.theta2=0;
        end
        try
            gp.r;
        catch
            gp.r=1;
        end
        for i=1:gp.N,
            d=(gp.v(i)-gp.v).^2;
            d=d./(gp.r^2);
            gp.C(i,:)=gp.theta1*exp(-0.5*d)'+gp.theta2;
        end
        
    case 'sqexp',
        try 
            gp.s;
        catch
            gp.s=0.1;
        end
        try
            gp.sigma;
        catch
            gp.sigma=1;
        end
        lambda2=gp.s^2;
        for i=1:gp.N,
            d=(gp.v(i)-gp.v).^2;
            d=d/lambda2;
            gp.C(i,:)=gp.sigma*exp(-0.5*d)';
        end
        gp.sigma=gp.C(1,1);
    otherwise
        disp('Unknown covariance function');
end
    

% Get equivalent kernel
%gp.A=chol(gp.C);

% Use eig method
[v,d]=eig(gp.C);
dd=diag(d);
negd=find(dd<0);
dd(negd)=0;
gp.Capp=v*diag(dd)*v';
gp.A=v*diag(sqrt(dd))*v';