function [E,V,eigen]=extension_prepare(G,H,kernel,sigma,condnumber,optarg);
% function extension_prepare: computes extension matrices for Gaussian and
%                               bandlimited Nystrom extension
%Syntax:
%   [E,V,eigen]=extension_prepare(G,H,kernel,sigma,condnumber,optarg);
%   IN:
%       G: is a m by d matrix of the small set
%       H: is a n by d matrix of the set where one wants to extend
%       kernel: choice of kernel. Must be 'gaussian' of 'bandlimited'
%       sigma: width of the kernel
%       condnumber: condition number for the extension procedure
%       optarg (optional): for bandlimited extension, is the dimension of
%                           the Bessel kernel
%   OUT:
%       E: the extension matrix, each column is the extension of a
%           geometric harmonic
%       V: the geometric harmonics on the set (rows)
%       eigen: the distribution of eigenvalues
%   Example:
%       [E,V,eigen]=extension_prepare(G,H,kernel,sigma,condnumber,optarg);
%       c=V*f; %decompose f on the inner basis
%       F=E*c; %extend f as F

[n,d]=size(H);
[m,tmp]=size(G);
if(tmp~=d)
    disp('Check dimension of G and H!')
    return;
end
if(nargin<6)
    optarg=d;
elseif(optarg<d)
    disp('Beware! You''re choosing a dimension parameter that is lower than the dimension of the ambiant space')
end
if (condnumber<1)
    disp('Condition number too small, must be at least equal to 1')
    return;
end
DG2=distances(G,G);
DHG2=distances(G,H);
switch(lower(kernel))
    case 'gaussian'
        K1=exp(-DG2/sigma^2);
        K2=exp(-DHG2/sigma^2);
    case 'bandlimited'
        DG2=DG2+1e-15;
        DHG2=DHG2+1e-15;
        K1=besselj(d/2,2*pi*sqrt(DG2))./(DG2.^(d/4));
        K2=besselj(d/2,2*pi*sqrt(DHG2))./(DHG2.^(d/4));
    otherwise
        disp('Unknown kernel')
        return;
end

[U,S,V]=svd(K1);
eigen=diag(S);
ind=find(condnumber*eigen>=eigen(1));
U=U(:,ind)*diag(1./eigen(ind));
E=K2*U;
V=V(:,ind)';
return;