function S=calc_S_new_data(X_new,A,V,Mu,ncomp,CV)

   % X_new(X_new==0) = eps;  
   M = ~isnan(X_new);
    X_new(isnan(X_new)) = 0;  
     Av=CV.A;
     
     
     [X_new,Xprobe] = SubtractMu( Mu, X_new, M, [], [], 1 );
     %for j = 1:n2
            A_j = repmat(M,1,ncomp) .* A;
            Psi = A_j' * A_j + diag( repmat(V,1,ncomp) );
            if ~isempty(Av)
                for i = find(M)'
                    Psi = Psi + Av{i};
                end
            end
            invPsi = inv(Psi);
            S = invPsi * A_j' * X_new;
            Sv{1} = V * invPsi;

           %posterior variances for A
    % end    
 
         %  [ dMu, A, Av, S, Sv ] = RotateToPCA( A, Av, S, Sv, CV.Isv, {}, 1 );
        
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find the PCA rotation: This has to be checked
function [ dMu, A, Av, S, Sv ] = ...
    RotateToPCA( A, Av, S, Sv, Isv, obscombj, update_bias );

n1 = size(A,1);
n2 = size(S,2);

if update_bias
    mS = mean(S,2);
    dMu = A*mS;
    S = S - repmat(mS,1,n2);
else
    dMu = 0;
end

covS = S*S';
if isempty(Isv)
    for j = 1:n2
        covS = covS + Sv{j};
    end
else
    nobscomb = length(obscombj);
    for j = 1:nobscomb
        covS = covS + ( length(obscombj{j})*Sv{j} );
    end
end
    
covS = covS / n2;
%covS = covS / (n2-n1);
[VS,D] = eig(covS);
RA = VS*sqrt(D);
A = A*RA;
covA = A'*A;
if ~isempty(Av)
    for i = 1:n1
        Av{i} = RA'*Av{i}*RA;
        covA = covA + Av{i};
    end
end
covA = covA / n1;
[VA,DA] = eig(covA);
[DA,I] = sort( -diag(DA) );
DA = -DA;
VA = VA(:,I);
A = A*VA;

if ~isempty(Av)
    for i = 1:n1
        Av{i} = VA'*Av{i}*VA;
    end
end
R = VA'*diag(1./sqrt(diag(D)))*VS';

S = R*S;
for j = 1:length(Sv)
    Sv{j} = R*Sv{j}*R';
end

end



%%%%%%
function [X,Xprobe] = SubtractMu( Mu, X, M, Xprobe, Mprobe, update_bias )

n2 = size(X,2);

if ~update_bias
    return
end
   
if issparse(X)
    X = subtract_mu( X, Mu );
    if ~isempty(Xprobe)
        Xprobe = subtract_mu( Xprobe, Mu );
    end
else
    X = X - repmat(Mu,1,n2).*M;
    if ~isempty(Xprobe)
        Xprobe = Xprobe - repmat( Mu, 1, n2 ).*Mprobe;
    end
end
end

end