function P = PPLS_sift_gen(X, pub, pri, ncomp, lambda)
% Privacy Partial Least Squares (PPLS)
%
% Input:
    % X: nxd matrix, n is the number of samples and d is the dimensionality
    % pub: NxD1 matrix, indicating the presence/absence of the D1 public attributes (via binary labels)
    % pri:  NxD1 matrix, indicating the presence/absence of the D1 private attributes (via binary labels)
    % ncomp: the maximum sift dimensionality desired (default is 10)
    % lamda: tradeoff parameter between privacy/utility used in the objective function (default is 1)
% Output:
    % P: sift (transformation) matrix
% Notes: 	
% This code is code modified from the SIMPLS subfunction in PLSREGRESS
%     http://www.mathworks.com/help/stats/plsregress.html
%
% If you use this code please cite the following paper: 
% SensorSift: Balancing Sensor Data Privacy and Utility 
%             in Automated Face Understanding,
%    Miro Enev, Jaeyeon Jung, Liefeng Bo, Xiaofeng Ren, Tadayoshi Kohno 
%    Annual Computer Security Applications Conference (ACSAC), 2012 

% Copyright (c) 2013, Miro Enev, Jaeyeon Jung, Liefeng Bo, Xiaofeng Ren, Tadayoshi Kohno 
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
% 1. Redistributions of source code must retain the above copyright
%    notice, this list of conditions and the following disclaimer.
% 2. Redistributions in binary form must reproduce the above copyright
%    notice, this list of conditions and the following disclaimer in the
%    documentation and/or other materials provided with the distribution.
% 3. All advertising materials mentioning features or use of this software
%    must display the following acknowledgement:
%    This product includes software developed by the Intel Science and 
%       Technology Center.
% 4. Neither the name of the Intel Science and Technology Center nor the
%    names of its contributors may be used to endorse or promote products
%    derived from this software without specific prior written permission.
% 
% THIS SOFTWARE IS PROVIDED BY US ''AS IS'' AND ANY
% EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
% DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
% DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
% (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
% LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
% ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

if nargin < 4, ncomp = size(pub, 2); end
if nargin < 5, lambda = 10; end

meanX = mean(X,1);
X0 = bsxfun(@minus, X, meanX);
P = simpls(X0,pub,pri,ncomp,lambda);

function P = simpls(X0,pub,pri,ncomp,lambda)

[n,dx] = size(X0);
dy = size(pub,2);

% Preallocate outputs
outClass = superiorfloat(X0,pub);
Xloadings = zeros(dx,ncomp,outClass);

% An orthonormal basis for the span of the X loadings, to make the successive
% deflation X0'*pub simple - each new basis vector can be removed from Cov
% separately.
V = zeros(dx,ncomp);

Cov = X0'*pub;
priCovsquare = (X0'*pri)*pri'*X0;
P = [];

for i = 1:ncomp
    % Find unit length ti=X0*ri and ui=pub*ci whose covariance, ri'*X0'*pub*ci, is
    % jointly maximized, subject to ti'*tj=0 for j=1:(i-1).
    [ri,si,ci] = svd(Cov,'econ'); ri = ri(:,1); ci = ci(:,1); si = si(1);
    
    % find transformation matrix    
    CovSquare = ( Cov*Cov' - lambda*priCovsquare );
    [ri, beta] = eigs( CovSquare, 1, 'LR');
    
    % sort eigenvectors (sift components)
    ri = real( ri );
    beta = real( diag( beta ) );
    [~, order] = sort( beta, 'descend' );
    
    % append h to sift matrix
    ri = ri( :, order(1) );
    P = [P ri];
    
    ti = X0*ri;
    normti = norm(ti); 
        ti = ti ./ normti; % ti'*ti == 1
    Xloadings(:,i) = X0'*ti;

    % Update the orthonormal basis with modified Gram Schmidt (more stable),
    vi = Xloadings(:,i);
    for repeat = 1:2
        for j = 1:i-1
            vj = V(:,j);
            vi = vi - (vi'*vj)*vj;
        end
    end
    vi = vi ./ norm(vi);
    V(:,i) = vi;

    % Deflate Cov, i.e. project onto the ortho-complement of the X loadings.
    % First remove projections along the current basis vector, then remove any
    % component along previous basis vectors that's crept in as noise from
    % previous deflations.
    Cov = Cov - vi*(vi'*Cov);
    Vi = V(:,1:i);
    Cov = Cov - Vi*(Vi'*Cov);
    disp([i ncomp]);
end
