function [model] = SmoothPCA(X,time,LV,lambda,deriv,group)
% [model] = SmooothPCA(X,T,LV,lambda,deriv,group);
% 
% X         data
% time      column with times of measurement
% LV        number of components
% lambda    Can be a scalar or column vector with lambda values
% deriv     Order of derivative for smooting (only 1 and 2 can be chosen)
% group     vector containing object-identification-information
%
% model-structure contains per value of lambda:
% Z         Final result for Z (smooth scores)
% P         Final result for P (loadings)
% RES       Sum of squared errors after each component (total Var(x) in 1)
% DD        Overall difference matrix
% ExpVar    %Explained variance per component
%
% Q = ||X - ZP'||^2 + lambda||DZ||^2
%
% JAW UVA-BDA, 14 February, 2008

% Here each component is separately estimated. Then after convergence, the
% component is deflated from the original data. I'm not sure whether this
% is the best option, maybe it is possible to simultaneously calculate all
% components, but then one should define a smooth surface. At the moment I
% don't know how to do that.

% The loadings are forced to be orthogonal by a projection of a new loading
% on the existing ones.

% smooth_pca_jw2 is different from smooth_pca_jw in that now the scores are
% not forced to come from the X-space.

% Note that this mfile can deal with irregular time sequence.
% In this case the differential matrix is multiplied by the inverse of the 
% time difference between the measurements. 

% Adapted:
% 020506 CMR %ExpVar toegevoegd
% 020506 CMR ongelijk aantal timepoints per object (D matrix per object)
%            (check: DD is gelijk aan D uit smooth_pca2_missing als group
%            is een kolom met enkel 1-en)
% 030506 CMR lambda als vector toegevoegd


I = size(X,1);
conv_c = [1E-6,500];

totaal = 0;
for g = 1:max(group)
    D=[];
    aantal = length(group(find(group == g)));
    D=zeros(aantal-2,aantal);
    totaal = totaal+aantal;
    Tg     = time(find(group == g));
    if deriv == 1;
        % Make 1st order diff matrix
        for i = 2:aantal-1
            h=Tg(i)-Tg(i-1);
            k=Tg(i+1)-Tg(i);
            D(i-1,i-1)=-1/(h+k);
            D(i-1,i+1)=1/(h+k);
        end
        
    end
    if deriv == 2;
        % Make diff matrix in case of 2nd order diff penalty for unequidistant
        % sampling
        for i = 2:aantal-1
            h = Tg(i) - Tg(i-1);
            k = Tg(i+1) - Tg(i);
            D(i-1,i-1) = 2*h / (k*h*h + h*k*k);
            D(i-1,i) = -2*(k+h) / (k*h*h + h*k*k);
            D(i-1,i+1) = 2*k / (k*h*h + h*k*k);
        end
    end
    
    % DD opvullen
    if deriv==1
        if g==1
            DD=D;
        else
            [DD_I,DD_J]=size(DD);
            DD = [DD zeros(DD_I,aantal);zeros(aantal-2,DD_J) D];
        end
    end

    if deriv==2
        if g==1
            DD=D;
        else
            [DD_I,DD_J]=size(DD);
            DD = [DD zeros(DD_I,aantal);zeros(aantal-2,DD_J) D];
        end
    end
end
% opslaan diff-matrix
model.D=DD;

for l=1:length(lambda);
    Z = [];P = [];RES = []; ExpVar=[];
    XX = X;
    for lv = 1:LV
        [U,S,V] = svds(XX,1);
        z = U*S;
        p = V;
        kount = 0;
        sse=1E6;
        XXadj = pinv(eye(I) + lambda(l)*DD'*DD)*XX;

        while (kount==0 || (SSEdiff>conv_c(1) && kount<conv_c(2)))

            kount = kount+1;

            % Estimate p
            p = (z\XX)';
            if lv > 1
                p=p-P*P'*p;
            end
            p = p/sqrt(p'*p);

            % Estimate z
            z = XXadj*p;

            % Check convergence
            E = XX - z*p';
            SSEdiff = abs(sse - trace(E'*E))/sse;
            sse = trace(E'*E);

        end
        XX = XX - z*p';
        Z = [Z z];P = [P p];
        RES = [RES sse];
    end
    % info per lambda wegschrijven)
    a=genvarname(['lambda' num2str(l)]);
    eval(['model.' a '.lambda = lambda(l);']);
    eval(['model.' a '.Z = Z;'])
    eval(['model.' a '.P = P;'])

    RES = [trace(X'*X) RES];
    for i =1:size(RES,2)-1;
        expvar = (RES(i)-RES(i+1))/RES(1);
        ExpVar=[ExpVar;expvar];
    end
    eval(['model.' a '.RES = RES;'])
    eval(['model.' a '.ExpVar = ExpVar;'])
end