function model = spca(X,k,sumabsv,options)
% model = spca(X,k,sumabsv,options)
%
% Calculates a sparse PCA model (with penalty on loadings).
% The components are either nested (i.e. based on current residuals), or
% calculated simultanousely for all components.
% The vectors/matrixes (scores and loadings) are calculated alternating
% (i.e. fix one calculate the other - switch)
%
% input:    X - data (the algorithm do not center or scale - so do that in
%               advance if necessary).
%           k - number of components/clusters
%           sumabsv - sum of the absolute values of the (normalized)
%           loading vector (1 < sumabst < sqrt(size(X,2)))
%           options - optional struct.
%
% output:   model.T - normalized scores. to get variance incorporated in
%           the scores do: T = model.T*diag(model.D);
%           model.V - normalized loadings
%           model.D - diagonal elements
%           model.details - details on model
% options:
%           algorithm:  [{'asls'} | 'defl']    Algorithm used for
%                       estimation of components. 'asls': Alternating
%                       Schrunken Least Squares, calculates the entire set
%                       of components simultanously. 'defl': calculates one
%                       component at a time based on the current residuals:
%                       E = X - tp';
%           orthscores: [{0} | 1] Force the scores to be mutually
%                       orthogonal. A model with orthogonal scores can not
%                       be used for an unique prediction of new samples.
%           v_initial:  [{'svd'} | V]   Pass a loading matrix (V) as
%                       initals. If nothing is specified SVD loadings are
%                       used.
%
%I/O: model     = spca(X,ncomp,sumabsv,options);    % calibration
%I/O: pred      = spca(X,model);                    % projection of new X data onto existing model
%I/O: options   = spca('options');                  % returns default options structure
%
% Morten Rasmussen, mortenr@life.ku.dk, Copenhagen University 2011

if ischar(X)
    options = [];
    options.name          = 'options';
    options.algorithm     = 'asls';
    options.v_initial = 'svd';
    options.orthscores = 0;
    model = options;
    return;
end

if nargin==1;
    error('At least 3 inputs are required')
end

if nargin==2; %prediction
    if isstruct(k);
        disp('Performing prediction');
        model = spca_pred(X,k);
        return
    else
        error('A double and a (model) struct is required for prediction, otherwise minimum 3 inputs for calibration')
    end
end

if nargin<4;
    options = spca('options');
end

algorithm = options.algorithm;

idnan = isnan(X);
if any(idnan(:));
    error('The current version is not able to handle missing values... The quick and dirtly trick is to impute some PCA guesses')
end

switch algorithm
    case {'asls'}
        disp('The entire set of components are estimates simultanously')
        model = spca_asls(X,k,sumabsv,options);
    case {'defl'}
        disp('Components are estimates via deflation')
        model = spca_defl(X,k,sumabsv,options);
end

%%%%%%%%%%%%
function model = spca_asls(X,k,sumabsv,options)
% model = spca(X,k,sumabsv,options)
%
% Calculates a sparse PCA model (with penalty on loadings).
% The components calcualted simultanously.
% The matrices (scores and loadings) are calculated alternating (i.e.
% fix on calculate the other - switch)
%
% input:    X - data (the algorithm do not center or scale - so do that in advance if necessary)
%           k - number of components/clusters
%           sumabsv - sum of the absolute values of the (normalized)
%           loading vector (1 < sumabst < sqrt(size(X,2)))
%           options - optional struct.
%
% output:   model.T - normalized scores
%           model.V - normalized loadings
%           model.D - diagonal elements
%           model.details - details on model


if ischar(options.v_initial)
    [u s V] = svds(X,k);
else
    V = options.v_initial;
end

orthT = options.orthscores;

% calculate all together without deflation.
% initial seettings
n = size(X,1);
vold = V; v = vold; 
told = zeros(n,k);
convth = 1e-10;
difft = 1;
diffv = 1;
count = 0;
maxit = 10000;
oldfit = 0;
fit = 0;
diff = 1;
E = X;
sse0 = sum(E(:).^2);


% while (difft>convth | diffv>convth) & count<maxit;
while diff>convth & count<maxit;
    count = count +1;
    
    % calculate score - LS estimate;
    t = [];
    if orthT==0;
        t = X*pinv(vold)';
    elseif orthT==1;
        for i=1:k;
            if i==1;
                t(:,i) = X*vold(:,1);
            else
                t(:,i) = (eye(n) - t*t')*X*vold(:,i);
            end
            t(:,i) = t(:,i)./sqrt(t(:,i)'*t(:,i));
        end
    end
    
    % normalize
    t = t*diag(1./sqrt(diag(t'*t)));
    
    % calculate loads
    [v Dv] = SMRupdate3(X',v,t,0,sumabsv,0); 
    d = diag(t'*X*v);
    
    % check convergence
    E = X-t*diag(d)*v';
    sse = sum(E(:).^2);
    fit = 1 - sse/sse0;
    diff = fit - oldfit;
    
    vold = v;
    told = t;
    oldfit = fit;
end

d = diag(t'*X*v);

model.T = t;
model.V = v;
model.D = d;
ev = calcexplainedvar(X,model);
model.explvar.tot = ev.tot;
model.explvar.prcomp = ev.prcomp;
model.details.orthscores = orthT;
model.details.sumabsv = sumabsv;
model.details.algorithm = 'asls';
model.details.numiter = count;

function model = spca_defl(X,k,sumabsv,options)
% model = spca(X,k,sumabsv,options)
%
% Calculates a sparse PCA model (with penalty on loadings).
% The components are nested (i.e. based on current residuals).
% The vectors (scores and loadings) are calculated alternating (i.e.
% fix on calculate the other - switch)
%
% input:    X - data (the algorithm do not center or scale - so do that in advance if necessary)
%           k - number of components/clusters
%           sumabsv - sum of the absolute values of the (normalized)
%           loading vector (1 < sumabst < sqrt(size(X,2)))
%           options - optional struct.
%
% output:   model.T - normalized scores
%           model.V - normalized loadings
%           model.D - diagonal elements
%           model.details - details on model

if ischar(options.v_initial)
    [u s V] = svd(X);
else
    V = options.v_initial;
end

orthT = options.orthscores;

Xres = X;
for i=1:k;
    vinit = V(:,i);
    % calculate component
    if i>1 & orthT==1;
        out = singleL1orth(Xres,sumabsv,T,vinit);
    else
        out = singleL1(Xres,sumabsv,vinit);
    end
    % collect
    T(:,i) = out.t;
    VV(:,i) = out.v;
    D(i) = out.d;
    %deflate
    Xm = out.t*out.d*out.v';
    Xres = Xres - Xm;
end

model.T = T;
model.V = VV;
model.D = D';
ev = calcexplainedvar(X,model);
model.explvar.tot = ev.tot;
model.explvar.prcomp = ev.prcomp;
model.details.orthscores = orthT;
model.details.sumabsv = sumabsv;
model.details.algorithm = 'defl';

function out = singleL1(X,sumabsv,v)
% calculate single component with L1 constraint in loading mode.
% input: X - data, sumabsv - constraint, v - initial guess on loadings
% (optional).

n = size(X,1);
if nargin<3;
    % calculate initials if not provided
    [u s v] = svds(X,1);
end

vold = v;
told = zeros(n,1);

% some settings
convth = 1e-10;
difft = 1;
diffv = 1;
count = 0;
maxit = 10000;
nonneg = 0; 

while (difft>convth | diffv>convth) & count<maxit;
    count = count +1;
    
    % calculate score - LS estimate;
    t = X*v;
    t = t./sqrt(t'*t);
    
    % calculate loads
    % LS estimate;
    v = t'*X; v = v';
    % soft thresshold v
    % find lambda fullfilling sum(abs(v))<=sumabsv
    l = searchL(v,sumabsv,nonneg);
    v = softth(v,l);
    v = v./sqrt(v'*v);
    
    % check convergence
    difft = told - t; difft = difft'*difft;
    diffv = vold - v; diffv = diffv'*diffv;
    
    vold = v;
    told = t;
end

out.t = t;
out.v = v;
out.d = t'*X*v;

function out = singleL1orth(X,sumabsv,T,v)
% calculate single component with L1 constraint in loading mode.
% input: X - data, sumabsv - L1 constrain on loadings, T - the k-1 scores,
% v - initial guess on loadings (optional).

n = size(X,1);
if nargin<4;
    % calculate initials if not provided
    [u s v] = svds(X,1);
end

vold = v;
told = zeros(n,1);

% some settings
convth = 1e-10;
difft = 1;
diffv = 1;
count = 0;
maxit = 10000;

while (difft>convth | diffv>convth) & count<maxit;
    count = count +1;
    
    % calculate score - LS estimate;
    t = (eye(n) - T*T')*X*v;
    t = t./sqrt(t'*t);
    
    % calculate loads
    % LS estimate;
    v = t'*X; v = v';
    % soft thresshold v
    % find lambda fullfilling sum(abs(v))<=sumabsv
    l = searchL(v,sumabsv);
    v = softth(v,l);
    v = v./sqrt(v'*v);
    
    % check convergence
    difft = told - t; difft = difft'*difft;
    diffv = vold - v; diffv = diffv'*diffv;
    
    vold = v;
    told = t;
end

out.t = t;
out.v = v;
out.d = t'*X*v;


function xst = softth(x,l,nonneg)
% soft threshold of x by l;
if nargin==2;
    nonneg=0;
end
xst=sign(x).*max(0, abs(x)-l);
if nonneg==1;
    xst(xst<0) = 0;
end

function l=searchL(t,sumabst,nonneg)

if norm(t,2)==0 || sum(abs(t./norm(t,2)))<=sumabst
    l=0;
    return
end
l1 = 0;
l2 = max(abs(t))-1e-4;

iter = 1;
maxit = 1000;
while iter < maxit
    % make a new guess
    lnew = (l1+l2)/2;
    stnew = softth(t,lnew,nonneg);
    % check it
    if sum(abs(stnew/norm(stnew,2)))<sumabst
        l2 = lnew;
    else
        l1= lnew;
    end
    if (l2-l1)<1e-5
        l=lnew;
        return
    end
    iter = iter+1;
end




function pred = spca_pred(X,model)
% predicts scores for a new sample (X) based on a calibration model (model)

if model.details.orthscores == 1;
    error('There is no unique solution for new sample in the case where the calibration model has orthogonal scores');
end

t = X*pinv(model.V)';
% normalize
t = t*diag(1./model.D);

pred.t = t;
pred.d = model.D;

function ev = calcexplainedvar(X,model)
% calculates the explained variance (in %) for the entire model (ev.tot) and for
% the individual components. OBS: Beaware that the component wise ev's are
% NOT additive as in ordinary PCA.

[n p] = size(X);

E = X - model.T*diag(model.D)*model.V';
if n<=p;
    SStot = trace(X'*X);
    SSe = trace(E'*E);
else
    SStot = trace(X*X');
    SSe = trace(E*E');
end

ev.tot = 1 - SSe/SStot;

ncomp = size(model.T,2);

for i=1:ncomp;
    E = X - model.T(:,i)*model.D(i)*model.V(:,i)';
    if n<=p;
        SSe = trace(E'*E);
    else
        SSe = trace(E*E');
    end
    evprcomp(i) = 1 - SSe/SStot;
end

ev.prcomp = evprcomp;


function [A Da] = SMRupdate3(X,A,B,nonneg,sumabsa,orth)
% Sparse Matrix Regression: [A Da] = SMRupdate3(X,A,B,nonneg,sumabsa,orth)
% X (I,JK) is the unfolded matrix of x (I,J,K);
% B (J,k) is loads of second mode
% A is the loads of the first mode, and the one that is estimated.
% nonneg is {0,1} for nonnegativity - OBS: do not work with active sparsity
% and nonnegativevity constraints
% sumabsa is the maximum L1 norm of of L1 normalized vectors of A = [a1,..,ak];
% 1 < sumabs < sqrt(I);

if nargin<5;
    orth=0;
end

k = size(B,2);
[n p] = size(X);
idnan = isnan(X);
E = X(~idnan);
sse0 = sum(E(:).^2);
convth = 1e-9;
diff = 1;
maxit = 1000;
oldfit = 0;
c = 0;
while diff>convth & c<maxit;
    c = c+1;
    for i=1:k;
        Bi = B(:,i);
        ic = true(k,1);
        ic(i) = false;
        Anoti = A(:,ic);
        Bnoti = B(:,ic);
        if any(idnan(:))
            for j=1:n;
                data = X(j,:);
                data(idnan(j,:)) = 0;
                D = diag(~idnan(j,:));
                Ails(j) = (data*Bi - Anoti(j,:)*Bnoti'*D*Bi)*pinv(Bi'*D*Bi);
            end
        else
            Ails = (X*Bi - Anoti*Bnoti'*Bi)*pinv(Bi'*Bi);
        end
        % find lambda fullfilling sum(abs(a))<=sumabsa
        l = searchL(Ails,sumabsa,0);
        Ails = softth(Ails,l,0);
        A(:,i) = Ails;
    end
    
    % check convergence
    E = X - A*B';
    E = E(~idnan);
    sse = sum(E(:).^2);
    fit = 1-sse/sse0    ;
    diff = fit - oldfit;
    oldfit = fit;
end

% normalize
Da = diag(sqrt(diag(A'*A)));
Dainv = pinv(Da);
A = A*Dainv;

%%%%%%%%%%%%%%%%%%%%

