function K=WDMR(Xe,Xv,lambda,k,r)
%WDMR computes the Weight Determination by Manifold Regularization kernel 
%matrix. 
%
%  K=WDMR(Xe,Xv,lambda,k,r) computes the kernel matrix K for Xe, an 
%  nx-by-Ne matrix, containing manifold constrained estimation regressors 
%  and Xv, an nx-by-Nv matrix, containing the validation regressors. lambda
%  and k are scalar design variables.


%% Compute M
[De,Ne] = size(Xe);
[Dv,Nv] = size(Xv);
X=[Xe Xv];
[D,N] = size(X);
% STEP1: COMPUTE PAIRWISE DISTANCES & FIND NEIGHBORS
X2= sum(X.^2,1);
distance = repmat(X2,N,1)+repmat(X2',1,N)-2*X'*X;
[sorted,index] = sort(distance);
neighborhood = index(2:(1+k),:);
% STEP2: SOLVE FOR RECONSTRUCTION WEIGHTS
W = zeros(k,N);
for ii=1:N
    zl = (X(:,neighborhood(:,ii))-repmat(X(:,ii),1,k));% shift ith pt to origin
    C = zl'*zl;                                        % local covariance
    C = C + eye(k,k)*r*trace(C)/k;                   % regularlization (K>D)
    W(:,ii) = C\ones(k,1);                           % solve Cw=1
    W(:,ii) = W(:,ii)/sum(W(:,ii));                  % enforce sum(w)=1
end
% STEP 3: COMPUTE EMBEDDING FROM EIGENVECTS OF COST MATRIX M=(I-W)'(I-W)
M = sparse(1:N,1:N,ones(1,N),N,N,4*k*N);
for ii=1:N
    w = W(:,ii);
    jj = neighborhood(:,ii);
    M(ii,jj) = M(ii,jj) - w';
    M(jj,ii) = M(jj,ii) - w;
    M(jj,jj) = M(jj,jj) + w*w';
end;

%% Compute Kernel
K=(1-lambda)*inv(lambda*M+(1-lambda)*[eye(Ne,Ne) zeros(Ne,Nv);zeros(Nv,Ne+Nv)]);
