function [  ] = test()
load('wiki_data.mat');
nBit = 32;
lambda = 1e-6;

XTrain = XAll(ITrain, :); YTrain = YAll(ITrain, :); LTrain = LAll(ITrain, :);
XTest = XAll(ITest, :); YTest = YAll(ITest, :); LTest = LAll(ITest, :);
dimX = size(XTrain, 2);
dimY = size(YTrain, 2);

XMean = mean(XTrain);
XTrain = bsxfun(@minus, XTrain, XMean);
XTest = bsxfun(@minus, XTest, XMean);

YMean = mean(YTrain);
YTrain = bsxfun(@minus, YTrain, YMean);
YTest = bsxfun(@minus, YTest, YMean);

Cxx = XTrain'*XTrain + lambda*eye(dimX);
Cyy = YTrain'*YTrain + lambda*eye(dimY);
Cxy = XTrain'*YTrain;

% evaluation
WTrue = (LTest * LTrain' >= 1);
MAPXY = zeros(3, 1);

[CCAWx, CCAWy, ~] = trainCCA(Cxx, Cyy, Cxy, nBit);
BxTest = compactbit(XTest * CCAWx >= 0);
ByTrain = compactbit(YTrain * CCAWy >= 0);
DHamm = hammingDist(BxTest, ByTrain);
[~, orderH] = sort(DHamm, 2);
MAPXY(1, 1) = calcMAP(orderH, WTrue);

[SCMsWx, SCMsWy] = trainSCM_seq(XTrain, YTrain, nBit, LTrain);
BxTest = compactbit(XTest * SCMsWx >= 0);
ByTrain = compactbit(YTrain * SCMsWy >= 0);
DHamm = hammingDist(BxTest, ByTrain);
[~, orderH] = sort(DHamm, 2);
MAPXY(2, 1) = calcMAP(orderH, WTrue);

[SCMoWx, SCMoWy] = trainSCM_orth(XTrain, YTrain, nBit, LTrain);
BxTest = compactbit(XTest * SCMoWx >= 0);
ByTrain = compactbit(YTrain * SCMoWy >= 0);
DHamm = hammingDist(BxTest, ByTrain);
[~, orderH] = sort(DHamm, 2);
MAPXY(3, 1) = calcMAP(orderH, WTrue);

MAPXY
end


% compute mean average precision (MAP)
function [MAP, succRate] = calcMAP (orderH, neighbor, topR)
[Q, N] = size(neighbor);
if (~exist('topR', 'var'))
    topR = N;
end
if (size(orderH, 1)~=Q || size(orderH, 2)~=N)
    disp(sprintf('size not match: %d=%d %d=%d', size(orderH, 1), Q, size(orderH, 2), N));
    MAP = 0;
    succRate = 0;
    return;
end
pos = [1: topR];
MAP = 0;
numSucc = 0;
for i = 1: Q
    ngb = neighbor(i, orderH(i, 1:topR));
    nRel = sum(ngb);
    if nRel > 0
        prec = cumsum(ngb) ./ pos;
        ap = mean(prec(ngb));
        MAP = MAP + ap;
        numSucc = numSucc + 1;
    end
end
MAP = MAP / numSucc;
succRate = numSucc / Q;
end

function cb = compactbit(b)
%
% b = bits array
% cb = compacted string of bits (using words of 'word' bits)

[nSamples nbits] = size(b);
nwords = ceil(nbits/8);
cb = zeros([nSamples nwords], 'uint8');

for j = 1:nbits
    w = ceil(j/8);
    cb(:,w) = bitset(cb(:,w), mod(j-1,8)+1, b(:,j));
end
end

function Dh=hammingDist(B1, B2)
%
% Compute hamming distance between two sets of samples (B1, B2)
%
% Dh=hammingDist(B1, B2);
%
% Input
%    B1, B2: compact bit vectors. Each datapoint is one row.
%    size(B1) = [ndatapoints1, nwords]
%    size(B2) = [ndatapoints2, nwords]
%    It is faster if ndatapoints1 < ndatapoints2
% 
% Output
%    Dh = hamming distance. 
%    size(Dh) = [ndatapoints1, ndatapoints2]

% example query
% Dhamm = hammingDist(B2, B1);
% this will give the same result than:
%    Dhamm = distMat(U2>0, U1>0).^2;
% the size of the distance matrix is:
%    size(Dhamm) = [Ntest x Ntraining]

% loop-up table:
bit_in_char = uint16([...
    0 1 1 2 1 2 2 3 1 2 2 3 2 3 3 4 1 2 2 3 2 3 ...
    3 4 2 3 3 4 3 4 4 5 1 2 2 3 2 3 3 4 2 3 3 4 ...
    3 4 4 5 2 3 3 4 3 4 4 5 3 4 4 5 4 5 5 6 1 2 ...
    2 3 2 3 3 4 2 3 3 4 3 4 4 5 2 3 3 4 3 4 4 5 ...
    3 4 4 5 4 5 5 6 2 3 3 4 3 4 4 5 3 4 4 5 4 5 ...
    5 6 3 4 4 5 4 5 5 6 4 5 5 6 5 6 6 7 1 2 2 3 ...
    2 3 3 4 2 3 3 4 3 4 4 5 2 3 3 4 3 4 4 5 3 4 ...
    4 5 4 5 5 6 2 3 3 4 3 4 4 5 3 4 4 5 4 5 5 6 ...
    3 4 4 5 4 5 5 6 4 5 5 6 5 6 6 7 2 3 3 4 3 4 ...
    4 5 3 4 4 5 4 5 5 6 3 4 4 5 4 5 5 6 4 5 5 6 ...
    5 6 6 7 3 4 4 5 4 5 5 6 4 5 5 6 5 6 6 7 4 5 ...
    5 6 5 6 6 7 5 6 6 7 6 7 7 8]);

n1 = size(B1,1);
[n2, nwords] = size(B2);

Dh = zeros([n1 n2], 'uint16');
for j = 1:n1
    for n=1:nwords
        y = bitxor(B1(j,n),B2(:,n));
        Dh(j,:) = Dh(j,:) + bit_in_char(y+1);
    end
end
end
