function [means, bases]=subspace_bases(X, sampleLabels, energyLossThreshold, affine)
% Computes the means and bases of the subspaces in a mixed linear or affine dataset.
%
% inputs:
% X: the data set, columns are samples.
% sampleLabels: the group membership of each sample.  Must be positive integers, starting at 1, no gaps.
% energyLossThreshold: the allowable MSE in the data set.  (i.e. epsilon^2)
% affine: a flag that should be set to true for affine subspaces, false for linear.
%
% outputs:
% means: a cell array of group mean vectors.
% bases: a cell array of group basis vectors.  For a given cell, the columns of the contained array are the basis vectors.

[ambientDimension, sampleCount]=size(X);
subspaceCount = max(sampleLabels);

means=cell(subspaceCount, 1);
bases=cell(subspaceCount, 1);
for subspaceIndex=1:subspaceCount,
    groupX=X(:,sampleLabels==subspaceIndex);
    groupSize=size(groupX,2);
    if groupSize==0
        error(['No group can be empty!']);
    else
        if affine
            groupXMean=mean(groupX,2);
            groupX=groupX-groupXMean*ones(1,groupSize);
            means{subspaceIndex}=groupXMean;
        end;
        [U,S,V]=svds(groupX,min(ambientDimension,groupSize));
        sDiag=diag(S);
        uSize=subspace_dimension(sDiag,groupSize*energyLossThreshold);
%        distances=point_to_space_distance(groupX,U(:,1:uSize));
%        totalSE = sum(distances.*distances);         
        bases{subspaceIndex}(:,1:uSize)=U(:,1:uSize);
    end
end
