function k=bhatta(X1,X2,eta)
%
% function k=bhatta(X1,X2,eta)
%
% X1 is a set of points as rows in a matrix
% X2 is another set of points as rows in a matrix
% eta is a scalar parameter
%
% Requires specifying a function called kernel.m that accepts two
% vectors and outputs a scalar.
%
% Based on the following paper (please cite it if you use this code):
%
% R. Kondor and T. Jebara. A Kernel between Sets of Vectors.
% International Conference on Machine Learning (ICML), August 2003.
%

[n1,dummy]=size(X1);
[n2,dummy]=size(X2);
n=n1+n2;

X=[X1;X2];

K=zeros(n,n);
for i=1:n
for j=1:i
K(i,j)=kernel(X(i,:),X(j,:));
K(j,i)=K(i,j);
end;
end;

M=eye(n);
for i=1:n
for j=1:i-1
M(i,j)=-kernel(X(i,:),X(j,:));
end;
M(i,:)=M(i,:)/sqrt(M(i,:)*K*M(i,:)');
end;

V=K*M';
V1=V(1:n1,:);
V2=V(n1+1:n,:);
m1=(1/n1)*V1'*ones(n1,1);
m2=(1/n2)*V2'*ones(n2,1);

S1=V1'*(1/n1*eye(n1)-1/n1^2*ones(n1))*V1+eta*eye(n);
S2=V2'*(1/n2*eye(n2)-1/n2^2*ones(n2))*V2+eta*eye(n);

S1inv=inv(S1);
S2inv=inv(S2);

mdag=(S1inv*m1+S2inv*m2)/2;
Sdag=2*inv(S1inv+S2inv);

a=-m1'*S1inv*m1/4;
b=-m2'*S2inv*m2/4;
c=mdag'*Sdag*mdag/2;

k=det(S1)^(-1/4)*det(S2)^(-1/4)*det(Sdag)^(1/2)*exp(a+b+c);