% ========================================================================
% Relaxed Collaborative Representation, Version 1.0
% Copyright(c) 2012 Meng Yang, Lei Zhang, and David Zhang
% All Rights Reserved.
%
% ----------------------------------------------------------------------
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is here
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
%----------------------------------------------------------------------
%
% This is an implementation of the algorithm for Face Recognition on AR
% with disguise.
%
% Please refer to the following paper
%
% Meng Yang, Lei Zhang, David Zhang, and Shenlong Wang, "Relaxed Collaborative
% Representation for Patter Classification," In Proc. CVPR 2012.
% 
% Contact: {csmyang,cslzhang}@comp.polyu.edu.hk
%----------------------------------------------------------------------

clear all;
clc;

%----RCR's parameter setting---- %
lambda            =   0.0005;
tau               =   0.005;
gamma             =   0.001;
Iteration         =   3;
simple_v          =   0; % 1: simple version; 0:original version
disguise          =   1; % 1: sunglass; 0: scarf

%Other parameters;
par.nClass        =   100;
par.height        =   83;
par.weight        =   60;
par.nDim          =   par.height*par.weight;
Block_num_h       =   4;
Block_num_w       =   2;
Block_num         =    Block_num_h*Block_num_w;

%loading data
load('AR_database');
Tr_DAT   =   []; trls     =   []; 
for ci = 1:100
    Tr_DAT = [Tr_DAT Tr_dataMatrix(:,1+7*(ci-1):4+7*(ci-1))];
    Tr_DAT = [Tr_DAT Tt_dataMatrix(:,1+7*(ci-1):4+7*(ci-1))];
    trls   = [trls repmat(ci,[1 8])];
end
clear Tr_dataMatrix Tr_sampleLabels Tt_dataMatrix Tt_sampleLabels

load('AR_database_Occlusion.mat');
Tt_DAT_sunglass = []; ttls_sunglass = [];
Tt_DAT_scarf    = []; ttls_scarf    = [];
for ci = 1:100
    Tt_DAT_sunglass = [Tt_DAT_sunglass Tr_dataMatrix(:,1+6*(ci-1):1+6*(ci-1))];
    Tt_DAT_sunglass = [Tt_DAT_sunglass Tt_dataMatrix(:,1+6*(ci-1):1+6*(ci-1))];
    ttls_sunglass   = [ttls_sunglass repmat(ci,[1 2])];
    
    Tt_DAT_scarf    = [Tt_DAT_scarf Tr_dataMatrix(:,4+6*(ci-1):4+6*(ci-1))];
    Tt_DAT_scarf    = [Tt_DAT_scarf Tt_dataMatrix(:,4+6*(ci-1):4+6*(ci-1))];
    ttls_scarf      = [ttls_scarf repmat(ci,[1 2])];
end
clear Tr_dataMatrix Tr_sampleLabels Tt_dataMatrix Tt_sampleLabels
if disguise
Tt_DAT = Tt_DAT_sunglass;
ttls   = ttls_sunglass;
else
Tt_DAT = Tt_DAT_scarf;
ttls   = ttls_scarf;
end

% partition training set
Tem_DAT = [];
for i = 1:size(Tr_DAT,2)
    tem = imresize(reshape(Tr_DAT(:,i),[165 120]),[83 60]);
    Tem_DAT(:,i) = tem(:);
end
Tr_DAT = Tem_DAT;
Tem_DAT = [];
for i = 1:size(Tt_DAT,2)
    tem = imresize(reshape(Tt_DAT(:,i),[165 120]),[83 60]);
    Tem_DAT(:,i) = tem(:);
end
Tt_DAT = Tem_DAT;

tr_dat  =  double(Tr_DAT);
tt_dat  =  double(Tt_DAT);
Tr_DAT  =  [];
step_h    =  floor(par.height/Block_num_h);
step_w    =  floor(par.weight/Block_num_w);

n_tem  = 1;
num_tr = size(tr_dat,2);
tr_dat = reshape(tr_dat,[par.height par.weight num_tr]);
for h_tem = 1:Block_num_h
for w_tem = 1:Block_num_w
    tem_dat = tr_dat((h_tem-1)*step_h+1:h_tem*step_h,(w_tem-1)*step_w+1:w_tem*step_w,:);
    Tr_DAT(:,:,n_tem) = reshape(tem_dat,[step_h*step_w,num_tr]);
    n_tem = n_tem + 1;
end
end
for bi = 1:Block_num
    tem = Tr_DAT(:,:,bi);
    Tr_DAT(:,:,bi)  =  tem./(repmat(sqrt(sum(tem.*tem)), [size(tem,1),1]) );
end
tr_dat  =  single(Tr_DAT); clear Tr_DAT Tem_DAT Tt_DAT tem_dat;

par.Dlabels       =   trls;
[hh,ww,tem] = size(tr_dat);
ID = [];

% offline computing  
BlkD        = [];
Sum_inv     = single(zeros(ww));
for bi = 1:Block_num
     BlkD(bi).TT = single(tr_dat(:,:,bi)'*tr_dat(:,:,bi)) ;
     tem = inv(BlkD(bi).TT + (lambda+tau)*eye(ww));
     BlkD(bi).P = single(tem);
     Sum_inv    = single(Sum_inv + tau/Block_num*tem);
end
Mat_inv = inv(eye(ww) -Sum_inv);
    
%% RCR Testing  
for indTest = 1:size(tt_dat,2)
    fprintf(['Totalnum:' num2str(size(tt_dat,2)) 'Nowprocess:' num2str(indTest) '\n']);
    
    % partition testing data
    n_tem  = 1;
    tt_y   = reshape(tt_dat(:,indTest),[par.height par.weight 1]);
    tty    = [];
    for h_tem = 1:Block_num_h
    for w_tem = 1:Block_num_w
        tem_dat = tt_y((h_tem-1)*step_h+1:h_tem*step_h,(w_tem-1)*step_w+1:w_tem*step_w,:);
        tty(:,n_tem) = tem_dat(:)./norm(tem_dat(:),2);
        n_tem = n_tem + 1;
    end
    end

    tic
    % initialize the coefficint
    for bi = 1:Block_num
        blkd(bi).P = BlkD(bi).P;
    end
    mat_inv = Mat_inv;
    for bi = 1:Block_num
        tem_v1(bi).v = tr_dat(:,:,bi)'*tty(:,bi);
        s(:,bi) = blkd(bi).P*tem_v1(bi).v;
        w(bi)   =  1;
        wb(bi)  =  0;
    end
    coef_m  =  mean(s,2);
    Diff_w   = norm(w-wb,2)./norm(w,2);
    
    % main loop
    iter     = 0;
    while Diff_w > 1e-3 & iter < Iteration
    % computing the weight
    for bi = 1:Block_num
        wb(bi) = w(bi);
        er(bi) = norm(s(:,bi)-coef_m,2)^2;
        w(bi) = exp(-1-tau*er(bi)/gamma);
    end
    sum_w = sum(w);
    w     = w*Block_num/sum_w;
    sum_w = Block_num;

    if ~simple_v % original version: computing projection matrix
    Sum_inv     = single(zeros(ww));
    for bi = 1:Block_num
        tem = inv(BlkD(bi).TT + (lambda + tau*w(bi))*eye(ww)); blkd(bi).P = single(tem);
        tem_ww(bi) = tau*w(bi)*w(bi)/sum_w;
        Sum_inv    = single(Sum_inv + tem*tem_ww(bi));
    end
    mat_inv = inv(eye(ww) - Sum_inv);
    end
    
    % computing the coefficients
    coef_mm = zeros(size(coef_m));
    for bi = 1:Block_num
        s(:,bi) = blkd(bi).P*tem_v1(bi).v;
        coef_mm  = coef_mm + w(bi)/sum_w*s(:,bi);
    end
    coef_m = zeros(size(coef_m));
    for bi = 1:Block_num
        tem_v = mat_inv*coef_mm;
        s(:,bi) = s(:,bi)+tau*w(bi)/sum_w*blkd(bi).P*tem_v;
        coef_m  = coef_m + w(bi)/sum_w*s(:,bi);
    end

    iter = iter + 1;
    Diff_w   = norm(w-wb,2)./norm(w,2);
    end
    
    % computing representation residual
    gap = zeros(1,par.nClass);
    for w_tem = 1:Block_num
    for indClass  =  1:par.nClass
        sc        =  s(:,w_tem);
        temp_s = sc(indClass==trls);
        zz                =  tty(:,w_tem)-tr_dat(:,indClass==trls,w_tem)*temp_s;
        gap(indClass)     =  gap(indClass)+zz(:)'*zz(:)*w(w_tem);      
    end
    end
    toc
    index = find(gap==min(gap));
    ID = [ID index(1)];
end

cornum = sum(ID==ttls);
if disguise
    fprintf('sunglass experiment:  ');
else
    fprintf('scarf experiment:  ');
end
fprintf('%s%f\n', ' Recogition rate = ', cornum/length(ttls));