% LOOCV2D_2Mag2013.m
%
% Esercise 3.(i) of 2 May 2013
%--------------------------------------------------------------------
% Script that performs leave-one-out cross-validation
% (Rippa's method) to find a good epsilon for 1D/2D RBF interpolation
% Calls on: DistanceMatrixCSRBF_new (in the case of Wendland) 
% --------------------------------------------------------------------

  clear all
  close all
  
  % Gaussian RBF
  %rbf = @(e,r) exp(-(e*r).^2); 
  
  % Wendland C^2
  rbf = inline('r.^4.*(5*spones(r)-4*r)','e','r'); 
  
  %testfunction = @(x,y) sinc(x).*sinc(y);
  
  testfunction = @(x) 0.75*(exp(-(9*x-2).^2)/4+exp(-(9*x+1).^2)/49)+ ...
  0.5*(exp(-(9*x-7).^2)/4)-0.1*exp(-(9*x-4).^2)/4;

  % Parameters for shape parameter loop below
  mine = 0; maxe = 10; ne = 50;
  ep = linspace(mine,maxe,ne);
  
  NV=[9 25 81 289 396];
  
  MN=length(NV);
  
  for j=1:5,
  % Number and type of data points
  
  N = NV(j);
  
  % equispaced points
  dsites=linspace(0,1,N)';
  
  ctrs = dsites;   % centers coincide with data sites
  % Create right-hand side vector, i.e.,
  % evaluate the test function at the data points.
  
  %rhs = testfunction(dsites(:,1),dsites(:,2));
  rhs = testfunction(dsites);
  for i=1:length(ep)
      
     % Compute distance matrix between the data sites and centers
     DM_data = DistanceMatrixCSRBF_new(dsites,ctrs,ep(i));
     
     % Compute interpolation matrix
     IM = rbf(ep(i),DM_data);
     % Compute error function (i.e., "cost" of epsilon)
     %invIM = pinv(IM);
     invIM=inv(IM); %this in case of Wendland
     EF = (invIM*rhs)./diag(invIM);
     % Compute maximum norm of EF
     maxEF(i) = norm(EF(:),inf);
  end
  
  
   fprintf('Smallest maximum norm: %e\n', min(maxEF));
   fprintf('at epsilon = %f\n',ep(maxEF==min(maxEF)));
  
  % Plot cost function norm
  figure; 
  semilogy(ep,maxEF,'b');
  
 
  end