classdef ndSparse
%ndSparse - A class of N-dimensional sparse arrays.
%
% by Matt Jacobson
%
% Copyright, Xoran Technologies, Inc. 2010
%
%
% USAGE:
%
% S=ndSparse(X) where X is an ordinary MATLAB sparse matrix converts X into
% an ndSparse object. S can be reshaped into an N-dimensional sparse array using
% its RESHAPE method, for arbitrary N.
%
% S=ndSparse(X,[M,N,P,...]) is equivalent to reshape(ndSparse(X),[M,N,P,...]).
%
%The class also has a variety of static methods that can be used to construct instances
%of the class. For example,
%
% S=ndSparse.build(Coordinates,Values,[m,n,p,...],nzmax)
%
%lets you generate an N-dimensional sparse array from a table of explicit entries.
%This is a generalization to N dimensions of S=sparse(i,j,s,m,n,nzmax).
%
%Other such methods include:
% ndSparse.accumarray
% ndSparse.sprand
% ndSparse.sprandn
% ndSparse.spalloc
%
%EXAMPLES:
%
% >> A=ndSparse.build( [1 1 1; 2 1 1;2 2 2] , [50,60 70]) %Builds a 2x2x2 sparse
% %array from table
%
% A(:,:,1) =
%
% (1,1) 50
% (2,1) 60
%
% A(:,:,2) =
%
% (2,2) 70
%
%Many of the same manipulations common to ordinary multidimensional MATLAB full arrays
%can be performed on the sparse 3D array A generated above. It can be permuted, summed,
%concatentated, and so forth e.g.,
%
% >> B=sum( permute([A,A+10], [3,2,1]) ,2)
%
% B(:,:,1) =
%
% (1,1) 120
% (2,1) 20
%
% B(:,:,2) =
%
% (1,1) 140
% (2,1) 160
%
%Other overloaded methods include BSXFUN, REPMAT, CIRCSHIFT, CONVN, FLIPDIMS, SQUEEZE,
%SHIFTDIM and many more. Type "methods ndSparse" for a full list and use
%"help ndSparse.methodname" to get details of usage.
%
%When browsing the list of methods, note that certain common operations have
%different implementations, optimized for different situations. Specifically,
%SUM, ANY,ALL, MIN, MAX... have alternative implementations SUMML, ANYML, ALLML,
%MINML, MAXML which are optimized for "low-dimensional" ndSparse objects OBJ.
%Here, low-dimensional means that a normal N-column MATLAB sparse matrix won't
%consume too much memory on your platform for N=MAX(NUMEL(OBJ)./SIZE(OBJ)).
%
%
%Another feature of the class is that bi-operand operations are allowed between ndSparse
%objects and MATLAB objects of any numeric type (single, uint16, etc...). This is not true
%of ordinary MATLAB sparse matrices, as of R2010b.
%
% >> C=eye(2,'single')*B(:,:,2)
%
% C =
%
% (1,1) 140
% (2,1) 160
%
% >> whos A B C
%
% Name Size Bytes Class Attributes
%
% A 2x2x2 136 ndSparse
% B 2x1x2 140 ndSparse
% C 2x1 104 ndSparse
%
%To convert back to an ordinary n-D full array, use the class' overloaded FULL method.
%To convert to a normal 2D sparse matrix, use the methods SPARSE or SPARSE2D. For example,
%SPARSE2D will convert an MxNxPx...xQ ndSparse array to the two dimensional
%(M*N*P*...)xQ sparse matrix in native MATLAB form.
properties (Constant)
oldstyle=false; %Set this to true so that sum, mean, any, all etc...
%use memory-liberal algorithms
end
properties (Access=private)
data;
ndShape;
end
methods
function obj=ndSparse(data,ndShape)
%ndSparse class constructor
%
% OBJ=ndSparse(A) where A is an ordinary MATLAB sparse matrix converts A into
% an ndSparse object. OBJ can be reshaped into an N-dimensional sparse array
% using its RESHAPE method, for arbitrary N.
%
% OBJ=ndSparse(A,[M,N,P,...]) is equivalent to reshape(ndSparse(A),[M,N,P,...])
if nargin<2,
ndShape=size(data);
elseif prod(ndShape)~=numel(data)
error 'The number of DATA elements is not consistent with the given ndShape.'
end
[mm,nn,ndShape]=ndShape2MN(ndShape);
if isa(data,'ndSparse')
obj=data;
if nargin>1
obj.data=reshape(obj.data,mm,nn);
obj.ndShape=ndShape;
end
return
end
if ~issparse(data),
data=reshape(data,mm,nn);
data=sparse(mkCompat(data));
end
data=reshape(data,mm,nn);
obj.data=data;
obj.ndShape=ndShape;
end
function obj=set.ndShape(obj,ndShape)
obj.ndShape=untrail1s(ndShape(:).');
end
function varargout=size(obj,varargin)
%size - same syntax as for full arrays
varargout=parseSize(obj.ndShape,nargout,varargin{:});
end
function L=length(obj)
%length - same behavior as for full arrays
if isempty(obj)
L=0;
else
L=max(obj.ndShape);
end
end
function out=end(obj,K,N)
%end - sameindexing rules as for full arrays
sz= obj;
if N==1,
out=numel(sz); return
end
[c{1:N}]=size(obj);
c=[c{:}];
out=c(K);
end
function N=ndims(obj)
%ndims - same behavior as for full arrays
N=length(size(obj));
end
function N=nnz(obj)
%nnz - same behavior as for full arrays
N=nnz(obj.data);
end
function N=nzmax(obj)
%nzmax - same behavior as for ordinary 2D sparse arrays
N=nzmax(obj.data);
end
function s=nonzeros(obj)
%nnonzers - same behavior as for ordinary 2D sparse arrays
s=nonzeros(obj.data);
end
function N=numel(obj)
%numel - same behavior as for full arrays
N=numel(obj.data);
end
function bool=isempty(obj)
%isempty - same behavior as for full arrays
bool=isempty(obj.data);
end
function bool=islogical(obj)
%islogical - returns true of data is stored as type logical
bool=islogical(obj.data);
end
function bool=issparse(obj)
%issparse - returns true for ndSparse objects
bool=true;
end
function bool=isnumeric(obj)
%isnumeric - returns true for non-logical ndSparse objects
bool=~islogical(obj);
end
function bool=isfloat(obj)
%isfloat - returns true for non-logical ndSparse objects
bool=~islogical(obj);
end
function bool=isreal(obj)
%isreal - works as for normal arrays
bool=isreal(obj.data);
end
function bool=isequal(varargin)
%isequal - works as for normal arrays
N=length(varargin);
sizeCell=varargin;
bool=false;
for ii=1:N%for-loop instead of cellfun is deliberate
if ~isnumeric(varargin{ii}) return; end
sizeCell{ii}=size(varargin{ii});
end
if ~isequal(sizeCell{:}) return; end
sz=sizeCell{1};
for ii=1:N%for-loop instead of cellfun is deliberate
varargin{ii}=mkCompat(varargin{ii},sz);
end
bool=isequal(varargin{:});
end
function bool=isequalwithequalnans(varargin)
%isequalwithequalnans - works as for normal arrays
N=length(varargin);
sizeCell=varargin;
bool=false;
for ii=1:N%for-loop instead of cellfun is deliberate
if ~isnumeric(varargin{ii}) return; end
sizeCell{ii}=size(varargin{ii});
end
if ~isequal(sizeCell{:}) return; end
sz=sizeCell{1};
for ii=1:N%for-loop instead of cellfun is deliberate
varargin{ii}=mkCompat(varargin{ii},sz);
end
bool=isequalwithequalnans(varargin{:});
end
function obj=isnan(obj)
%isnan - works as for normal arrays
obj.data=isnan(obj.data);
end
function obj=isinf(obj)
%isinf - works as for normal arrays
obj.data=isinf(obj.data);
end
function out=isfinite(obj)
%isinf - works as for normal arrays
out=reshape( isfinite(obj.data) , obj.ndShape);
end
function obj=real(obj)
%real - works as for normal arrays
obj.data=real(obj.data);
end
function obj=imag(obj)
%imag - works as for normal arrays
obj.data=imag(obj.data);
end
function obj=conj(obj)
%conj - works as for normal arrays
obj.data=conj(obj.data);
end
function obj=abs(obj)
%abs - works as for normal arrays
obj.data=abs(obj.data);
end
function obj=sqrt(obj)
%sqrt - works as for normal arrays
obj.data=sqrt(obj.data);
end
function c=underlyingClass(obj)
%underlyingClass(obj) will return the storage type of the data (double/logical)
c=class(obj.data);
end
function out=logical(obj)
%logical - converts storage type to double
out=logical(obj.data);
end
function varargout=find(obj,varargin)
%find - same syntax as for full arrays
[varargout{1:nargout}]=find(obj.data,varargin{:});
if nargout>1,
II=varargout{1};
JJ=varargout{2};
ndSubs=IJ2ndCoords(II,JJ,obj.ndShape);
II=ndSubs{1};
if length(ndSubs)>2
JJ=sub2ind(obj.ndShape(2:end),ndSubs{2:end});
else
JJ=ndSubs{2};
end
varargout(1:2)={II,JJ};
end
end
function out=double(obj)
%double - converts storage type to double
out=double(obj.data);
end
function out=sparse(obj)
%sparse(obj) - converts ndSparse object of dimension MxNxPxQx....
%into an ordinary 2D sparse matrix of dimensions Mx(N*P*Q...)
ndShape=obj.ndShape;
out=reshape(obj.data,ndShape(1),prod(ndShape(2:end)));
end
function data=sparse2d(obj)
%sparse2d(obj) - converts ndSparse object of dimension MxNxPxQx...YxZ
%into an ordinary 2D sparse matrix of dimensions (M*N*P*Q*...*Y)xZ
data=sparse(obj.data);
end
function data=full(obj)
%full - same syntax as for 2D sparse matrices
data=reshape( full(obj.data), obj.ndShape);
end
function obj=spfun(fun,obj)
%spfun - same syntax as for 2D sparse matrices
obj.data=spfun(fun,obj.data);
end
function obj=spones(obj)
%spones - same syntax as for 2D sparse matrices
obj.data=spones(obj.data);
end
function obj=reshape(obj,varargin)
%reshape - same syntax as for full arrays
map=cellfun('isempty',varargin);
switch nnz(map)
case 0
case 1
missing=prod(obj.ndShape)/prod([varargin{:}]);
if mod(missing,1), error 'Bad reshape args'; end
varargin(map)={missing};
otherwise
error 'Too many empty args'
end
newshape=[varargin{:}];
obj=ndSparse(obj,newshape);
end
function objnew=subsref(obj,S)
%subsref - same indexing rules as for full arrays
if isempty(S.subs)
objnew=obj; return
end
[E,targetshape,~,other]=nd2matrixIndex(S.subs,obj.ndShape,'subsref');
if other.boolConsol
obj=reshape(obj,other.quasiShape);
end
data=builtin('subsref',obj.data,E);
objnew=ndSparse(data,targetshape);
end
function objnew=subsasgn(obj,S,rhs)
%subsasgn - same indexing rules as for full arrays
if isempty(S.subs)
error ' An indexing expression on the left side of an assignment must have at least one subscript.'
end
if isequal(rhs,[])%NULL ASSIGNMENT
sz=size(obj);
nn=length(sz);
N=length(S.subs);
dim=find(~strcmp(':',S.subs));
qq=length(dim);
if qq>1,
error 'A null assignment can have only one non-colon index.';
elseif qq==0
dim=1;
idx=':';
else
idx=S.subs{dim};
end
if 1<N && N<nn %Consolidated indexing of trailing dims
outcell=parseSize(size(obj),N);
sz=[outcell{:}];
obj=reshape(obj,sz);
nn=N;
end
if dim>nn
T0=1;
else
T0=1:sz(dim);
end
try
idx=T0(idx); %force idx to be numeric nonlogical
catch ME
strpres=@(t,p) ~isempty(strfind(t,p));
if strpres(ME.identifier, 'badsubscript')
if strpres(ME.message,'index must be a positive integer or logical')
error 'Subscript indices must either be real positive integers or logicals.'
elseif strpres(ME.message,'index out of bounds')
error 'Index exceeds matrix dimensions.'
end
else
error 'Bad subscript'
end
end
if dim>nn,
sz=trail1s(sz,dim);
sz(end)=0;
objnew=ndSparse.build(sz);
return;
end
if ndSparse.oldstyle%memory liberal style
data=ExtractDim1Reshaping(obj.data,dim,sz);
data(idx,:)=[];
sz(dim)=size(data,1);
data=invExtractDim1Reshaping(data,dim,sz);
objnew=ndSparse(data,sz);
else
if isequal(idx,':') || ( sz(dim)==1 && ~isempty(idx) )
sz(dim)=0;
objnew=ndSparse.build(sz);
else
[ndSubs,vals]=getEntryTable(obj,length(sz));
todelete=ismember(ndSubs{dim},idx);
ndSubs=[ndSubs{:}];
ndSubs(todelete,:)=[];
vals(todelete)=[];
T=T0;
T(idx)=[];
T0(T)=1:length(T);
ndSubs(:,dim)=T0(ndSubs(:,dim));
sz(dim)=length(T);
objnew=ndSparse.build(ndSubs , vals, sz);
end
end
elseif isempty(rhs)
error 'Improper null assignment. Right hand side must be [].'
else%NON-NULL,
[E,targetshape,boolGrow,other]=nd2matrixIndex(S.subs,obj.ndShape,'subsasgn');
if boolGrow
N = length(targetshape);
[ndSubs,vals]=getEntryTable(obj,N);
objnew=ndSparse.build(ndSubs,vals,targetshape,nzmax(obj));
elseif other.boolConsol
objnew=reshape(obj,other.quasiShape);
else
objnew=obj;
end
data=objnew.data;
if ~isscalar(rhs)
if prod(other.subshapeND)~=numel(rhs)
errorFlag=true;
else
errorFlag=false;
end
if other.LinearIndexing
if errorFlag, error ' In an assignment A(:) = B, the number of elements in A and B must be the same'; end
else
if errorFlag
error 'Subscripted assignment dimension mismatch.'
end
rhs=reshape(rhs,other.subshape);
end
end
data=builtin('subsasgn',data,E,rhs);
objnew=ndSparse(data , targetshape);
end
end
function idx=subsindex(obj)
%subsindex - same indexing rules as for full arrays
switch class(obj.data);
case 'double'
idx=obj.data-1;
case 'logical'
idx=find(obj.data)-1;
end
end
function objnew=sum(obj,varargin)
%sum - same syntax as for full arrays
objnew=sumEngine(0,obj,varargin{:});
end
function objnew=summl(obj,varargin)
%summl - An alternative implementation of SUM. It
%makes more liberal use of memory allocation, but is
%sometimes faster. Best suited for low dimensional arrays.
objnew=sumEngine(1,obj,varargin{:});
end
function objnew=any(obj,varargin)
%any - same syntax as for full arrays
objnew=anyEngine(0,obj,varargin{:});
end
function objnew=anyml(obj,varargin)
%anyml - An alternative implementation of ANY. It
%makes more liberal use of memory allocation, but is
%sometimes faster. Best suited for low dimensional arrays.
objnew=anyEngine(1,obj,varargin{:});
end
function objnew=all(obj,varargin)
%all - same syntax as for full arrays
objnew=allEngine(0,obj,varargin{:});
end
function objnew=allml(obj,varargin)
%allml - An alternative implementation of ALL. It
%makes more liberal use of memory allocation, but is
%sometimes faster. Best suited for low dimensional arrays.
objnew=allEngine(1,obj,varargin{:});
end
function objnew=mean(obj,varargin)
%mean - same syntax as for full arrays
objnew=meanEngine(0,obj,varargin{:});
end
function objnew=meanml(obj,varargin)
%meanml - An alternative implementation of MEAN. It
%makes more liberal use of memory allocation, but is
%sometimes faster. Best suited for low dimensional arrays.
objnew=meanEngine(1,obj,varargin{:});
end
function objnew=circshift(obj,varargin)
%circshift - same syntax as for full arrays
objnew=circshiftEngine(0,obj,varargin{:});
end
function objnew=circshiftml(obj,varargin)
%circshiftml - An alternative implementation of CIRCSHIFT. It
%makes more liberal use of memory allocation, but is
%sometimes faster. Best suited for low dimensional arrays.
objnew=circshiftEngine(1,obj,varargin{:});
end
function objnew=cat(obj,varargin)
%cat - same syntax as for full arrays
objnew=catEngine(0,obj,varargin{:});
end
function objnew=catml(obj,varargin)
%catml - An alternative implementation of CAT. It
%makes more liberal use of memory allocation, but is
%sometimes faster. Best suited for low dimensional arrays.
objnew=catEngine(1,obj,varargin{:});
end
function varargout=max(obj,varargin)
%max - same syntax as for full arrays
[varargout{1:nargout}]=maxEngine(0,obj,varargin{:});
end
function varargout=maxml(obj,varargin)
%maxml - An alternative implementation of MAX. It
%makes more liberal use of memory allocation, but is
%sometimes faster. Best suited for low dimensional arrays.
[varargout{1:nargout}]=maxEngine(1,obj,varargin{:});
end
function varargout=min(obj,varargin)
%min - same syntax as for full arrays
[varargout{1:nargout}]=minEngine(0,obj,varargin{:});
end
function varargout=minml(obj,varargin)
%minml - An alternative implementation of MIN. It
%makes more liberal use of memory allocation, but is
%sometimes faster. Best suited for low dimensional arrays.
[varargout{1:nargout}]=minEngine(1,obj,varargin{:});
end
function obj=horzcat(varargin)
%horzcat - same syntax as for full arrays
obj=cat(2,varargin{:});
end
function obj=vertcat(varargin)
%vertcat - same syntax as for full arrays
obj=cat(1,varargin{:});
end
function objnew=permute(obj,order)
%permute - same syntax as for full arrays
ndShape=obj.ndShape;
N=length(ndShape);
if length(order)<N,
error 'ORDER must have at least N elements for an N-D array'
else
[ndShape,N]=trail1s(ndShape,max(order));
end
if issorted(order), objnew=obj; return; end
[ndSubs,vals]=getEntryTable(obj,N);
ndSubs=ndSubs(:,order);
newshape=ndShape(order);
objnew=ndSparse.build(ndSubs,vals,newshape,nzmax(obj));
end
function objnew=ipermute(obj,order)
%ipermute - same syntax as for full arrays
[null,order]=sort(order);
objnew=permute(obj,order);
end
%%Unary ops
function obj=uplus(obj)
%uplus - same behavior as for full arrays
end
function obj=uminus(obj)
%uminus - same behavior as for full arrays
obj.data=-obj.data;
end
function obj=not(obj)
%not - same behavior as for full arrays
obj.data = ~obj.data ;
end
%%Binary ops
function obj=plus(L,R)
%plus - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) + mkCompat(R,s),s);
end
function obj=minus(L,R)
%minus - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) - mkCompat(R,s),s);
end
function obj=times(L,R)
%times - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) .* mkCompat(R,s),s);
end
function obj=rdivide(L,R)
%rdivide - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) ./ mkCompat(R,s), s);
end
function obj=ldivide(L,R)
%ldivide - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) .\ mkCompat(R,s),s);
end
function obj=power(L,R)
%power - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) .^ mkCompat(R,s),s);
end
function obj=lt(L,R)
%lt - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) < mkCompat(R,s),s);
end
function obj=le(L,R)
%le - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) <= mkCompat(R,s),s);
end
function obj=gt(L,R)
%gt - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) > mkCompat(R,s),s);
end
function obj=ge(L,R)
%ge - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) >= mkCompat(R,s),s);
end
function obj=eq(L,R)
%eq - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) == mkCompat(R,s),s);
end
function obj=ne(L,R)
%ne - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) ~= mkCompat(R,s),s);
end
function obj=and(L,R)
%and - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) & mkCompat(R,s),s);
end
function obj=or(L,R)
%or - same behavior as for full arrays
s=getshape(L,R);
obj = finalObject( mkCompat(L,s) | mkCompat(R,s),s);
end
%%Methods for 2D Matrices
function obj=inv(obj)
%inv - same behavior as for full arrays
if ndims(obj)>2, error 'Operation transpose defined for 2D arrays only'; end
obj = ndSparse(inv(sparse2d(obj)));
end
function obj=triu(obj,varargin)
%triu
if ndims(obj)>2, error 'Operation transpose defined for 2D arrays only'; end
obj = ndSparse( triu( sparse2d(obj) , varargin{:} ) );
end
function obj=tril(obj,varargin)
%tril - same behavior as for full arrays
if ndims(obj)>2, error 'Operation transpose defined for 2D arrays only'; end
obj = ndSparse( tril( sparse2d(obj) , varargin{:} ));
end
function obj=transpose(obj)
%transpose - same behavior as for full arrays
if ndims(obj)>2, error 'Operation transpose defined for 2D arrays only'; end
obj = ndSparse(sparse2d(obj).');
end
function obj=ctranspose(obj)
%ctranspose - same behavior as for full arrays
if ndims(obj)>2, error 'Operation ctranspose defined for 2D arrays only'; end
obj = ndSparse(sparse2d(obj)');
end
function obj=mtimes(L,R)
%mtimes - same behavior as for full arrays
if ( isscalar(R) ||isscalar(L))
obj=L.*R; return
end
if ~all([ndims(L),ndims(R)]==2) ,
error 'Operation mtimes defined for 2D arrays only'; end
obj = finalObject( mkCompat(L)*mkCompat(R));
end
function obj=mrdivide(L,R)
%mrdivide - same behavior as for full arrays
if ( isscalar(R))
obj=L./R; return
end
if ~all([ndims(L),ndims(R)]==2),
error 'Operation mrdivide defined for 2D arrays only'; end
obj = finalObject( mkCompat(L)/mkCompat(R));
end
function obj=mldivide(L,R)
%mldivide - same behavior as for full arrays
if ( isscalar(L))
obj=L.\R; return
end
if ~all([ndims(L),ndims(R)]==2) && ~isscalar(L),
error 'Operation mldivide defined for 2D arrays only'; end
obj = finalObject( mkCompat(L)\mkCompat(R));
end
function obj=mpower(L,R)
%mpower - same behavior as for full arrays
if ~all([ndims(L),ndims(R)]==2),
error 'Operation mpower defined for 2D arrays only'; end
obj = finalObject( mkCompat(L)^mkCompat(R));
end
function display(obj)
%display
l=inputname(1);
nn=ndims(obj);
simplecase=false;
if ~prod(obj.ndShape) %display empty array
T=evalc('full(obj)'); %only way to get at the builtin display method
simplecase=true;
elseif nn<=2
T=evalc('sparse2d(obj)'); %only way to get at the builtin display method
simplecase=true;
elseif ~nnz(obj)
dimchar=repmat({'-by-'},2,ndims(obj));
dimchar(1,:)=arrayfun(@(c) num2str(c), size(obj), 'Unif',0);
dimchar=[dimchar{1:end-1}];
T=evalc( ' ['' All zero sparse: '' dimchar] ' );
simplecase=true;
end
if simplecase
T=strrep(T,'ans =',[l ' =']);
jj=find(T~=sprintf('\n'),1,'last');
T=T(1:jj);
disp(T), disp ' '
return
end
[ndSubs,vals,ndShape]=getEntryTable(obj,ndims(obj));
trailingdims=ndShape(3:end);
data=sparse2d(obj);
sz=obj.ndShape(1:2);
M=prod(sz);
sheets=unique([ndSubs{:,3:end}],'rows');
sheetCorners=num2cell([ones(size(sheets,1),2), sheets], 1);
cornerRows=sub2ind(ndShape(1:end-1),sheetCorners{1:end-1});
cornerCols=sheetCorners{end};
rangecell=num2cell(sheets);
numSheets=size(rangecell,1);
commas=rangecell(1,:); commas(:)={','};
rangecell_aschar=cellfun(@(x)num2str(x),rangecell,'uni',0);
loopctr=0;
rg=0:M-1;
for ii=1:numSheets
loopctr=loopctr+1;
bookNum=cornerCols(ii);
rr=cornerRows(ii)+rg;
T=evalc('reshape(data(rr,bookNum),sz)'); %only way to get at the builtin display method
if ~isempty(l),
lbl=l;
%if numBooks>1
cc=[rangecell_aschar(loopctr,:); commas];
lbl=[l '(:,:,' cc{1:end-1} ')'];
%end
T=strrep(T,'ans =',[lbl ' =']);
end
jj=find(T~=sprintf('\n'),1,'last');
T=T(1:jj);
disp(T), disp ' '
end
end
function objnew=repmat(obj,varargin)
%repmat - same syntax as for full arrays
repDims=[varargin{:}];
ndShape=size(obj);
N=max(length(repDims),length(ndShape));
repDims=trail1s(repDims,N);
ndShape=trail1s(ndShape,N);
newshape=ndShape.*repDims;
if ~all(round(repDims)==repDims)
error 'Array repititions must be integer-valued'
end
if ~all(repDims)%some repDims are 0
objnew=ndSparse([],newshape); return
end
if ~any(repDims~=1),%No actual repeated data
objnew=obj;return
end
[ndSubs,vals,ndShape]=getEntryTable(obj,N);
if isempty(vals)
objnew=ndSparse.build(newshape); return
end
entrytable=[ndSubs{:},vals(:)];
for ii=1:N %loop over dims
thisRep=repDims(ii);
if thisRep==1, continue; end
repCoords=bsxfun(@plus, entrytable(:,ii), (0:thisRep-1)*ndShape(ii) );
entrytable=repmat(entrytable,thisRep,1);
entrytable(:,ii)=repCoords(:);
end
objnew=ndSparse.build(entrytable(:,1:end-1),...
entrytable(:,end),...
newshape);
end
function out=bsxfun(fun,L,R)
%bsxfun - same syntax as for full arrays
isL=isa(L,'ndSparse');
isR=isa(R,'ndSparse');
isfun=isa(fun,'function_handle');
if ~isfun, error 'First Argument must be a function handle.'; end
N=max(ndims(L),ndims(R));
[Lshape{1:N}]=size(L); Lshape=[Lshape{:}];
[Rshape{1:N}]=size(R); Rshape=[Rshape{:}];
if any(Lshape~=Rshape & Lshape~=1 & Rshape~=1)
error 'Non-singleton dimensions of the two input arrays must match each other.'
end
newshape=max([Lshape;Rshape], [],1);
newshape=newshape.*logical(Lshape.*Rshape); %mask out the zero-dimensions
Lreps=ones(size(newshape));
Lreps(Lshape==1)=newshape(Lshape==1);
Rreps=ones(size(newshape));
Rreps(Rshape==1)=newshape(Rshape==1);
if isL && isR
L=repmat(L,Lreps);
R=repmat(R,Rreps);
out=ndSparse(bsxfun(fun,L.data,R.data) , newshape);
return
elseif isL
Lsamp=L.data; Rsamp=R;
if isempty(Lsamp), Lsamp=reshape(Lsamp,0,1); else Lsamp=Lsamp(1); end
if isempty(Rsamp), Rsamp=reshape(Rsamp,0,1); else Rsamp=Rsamp(1); end
isOutFull=~issparse( bsxfun(fun, Lsamp, mkCompat(Rsamp) ) );
if isOutFull%cheapest to do everything as full array operations
out=bsxfun(fun,full(L),full(R)); return
else%sparse output - permutation approach is probably cheaper
L=repmat(L,Lreps);
idx=(Rshape~=1);
order=1:N;
order=[order(idx), order(~idx)];
L=permute(L,order);
R=permute(R,order);
data=bsxfun(fun,reshape(L.data,numel(R),[]),mkCompat(R(:)));
out=ndSparse( data , newshape(order));
out=ipermute(out,order);
return
end
else%isR
Lsamp=L; Rsamp=R.data;
if isempty(Lsamp), Lsamp=reshape(Lsamp,0,1); else Lsamp=Lsamp(1); end
if isempty(Rsamp), Rsamp=reshape(Rsamp,0,1); else Rsamp=Rsamp(1); end
isOutFull=~issparse( bsxfun(fun, mkCompat(Lsamp), Rsamp ) );
if isOutFull%cheapest to do everything as full array operations
out=bsxfun(fun,full(L),full(R)); return
else%sparse output - permutation approach is probably cheaper
R=repmat(R,Rreps);
idx=(Lshape~=1);
order=1:N;
order=[order(idx), order(~idx)];
L=permute(L,order);
R=permute(R,order);
data=bsxfun(fun,mkCompat(L(:)), reshape(R.data,numel(L),[]));
out=ndSparse( data , newshape(order));
out=ipermute(out,order);
return
end
end
end
function obj=squeeze(obj)
%squeeze - same syntax as for full arrays
ndShape=obj.ndShape;
ndShape(ndShape==1)=[];
obj=reshape(obj,ndShape);
end
function [obj,N]=shiftdim(obj,N)
%shiftdim - same syntax as for full arrays
if nargin>1
order=1:ndims(obj);
if N>=0
order=circshift(order,[0,-N]);
obj=permute(obj,order);
else
obj=reshape(obj, [ones(1,-N) ,obj.ndShape]);
end
else
sz=obj.ndShape;
N=find(sz~=1,1,'first');
obj=reshape(obj,sz(N:end));
N=N-1;
end
end
function objnew=flipdim(obj,dim)
%flipdim - same behavior as for full arrays
[ndSubs,vals,ndShape]=getEntryTable(obj,max(ndims(obj), dim));
ndSubs{dim}=ndShape(dim)+1 - ndSubs{dim};
nzm=nzmax(obj);
objnew=ndSparse.build([ndSubs{:}], vals, ndShape,nzm);
end
function objnew=flipud(obj)
%flipud - same behavior as for full arrays, except that it works
%for 3D and higher-dimensional arrays as well. That is,
%flipud(obj)=flipdim(obj,1)
objnew=flipdim(obj,1);
end
function objnew=fliplr(obj)
%fliplr - same behavior as for full arrays, except that it works
%for 3D and higher-dimensional arrays as well. That is,
%fliplr(obj)=flipdim(obj,2)
objnew=flipdim(obj,2);
end
function objnew=rot90(obj,K)
%rot90 - same behavior as for full arrays, except that it works
%for 3D and higher-dimensional arrays as well. The rotation will
%be applied to the first two dimensions.
if nargin<2, K=1; else K=mod(K,4); end
switch K
case 0
objnew=obj;
case 1
order=1:ndims(obj);
order([1 2])=order([2 1]);
objnew=permute(obj,order);
objnew=flipud(objnew);
case 2
objnew=flipud(fliplr(obj));
case 3
order=1:ndims(obj);
order([1 2])=order([2 1]);
objnew=permute(obj,order);
objnew=fliplr(objnew);
end
end
function obj=convn(varargin)
%convn method for ndSparse
%
%SYNTAXES:
%
%
%(1) OBJ=CONVN(A,B,SHAPE) performs the N-dimensional convolution of
% arrays A and B. By MATLAB dispatching rules, at least one of which
% A and B will be ndSparse.
%
%(2) OBJ=CONVN(K_1,K_2,...K_N,A,SHAPE) performs separable convolution of
% array A with kernels K_1,...,K_N. Each K_i must be a vector and will
% be applied along dimension i of A. Whether K_i is a row or column
% vector has no impact here on the effect of the convolution.
% For this syntax to be triggered, it is required that N>2 and that
% N=ndims(A)-isvector(A).
%
%The SHAPE input argument is optional and works like MATLAB's usual
%convolution functions:
%
%
% 'full' - (default) returns the full N-D convolution
% 'same' - returns the central part of the convolution that
% is the same size as A.
% 'valid' - returns only the part of the result that can be
% computed without assuming zero-padded arrays.
%
%
%CAUTION: The output of this routine will always be ndSparse and is intended
%for situations where the output data is expected to be sparse
%intrinsically, e.g., when sparse data is convolved with small or
%sparse kernels. If there is reason to think the output will not be sparse,
%it would be advisable to pre-cast to full arrays and use the native
%MATLAB convn method.
if ischar(varargin{end})
truncType=varargin{end}; varargin(end)=[];
else
truncType='full';
end
B=varargin{end};
varargin(end)=[];
B=ndSparse(B);
sizeB=size(B);
ndimsB=length(sizeB);
A=varargin{1};
nv=length(varargin);
sepFlag=(nv>1);%separable input kernel
if sepFlag
nvReq=ndimsB-isvector(B);
if nv~=nvReq
error('The number N of separable kernels K_i in CONVN(K_1,K_2,...,K_N,A) must equal N=ndims(A)-isvector(A)');
end
if ~all(cellfun(@(c) isvector(c) , varargin) );
error 'Separable kernel data appears to have been input. All the 1D kernels need to be vectors'
end
varargin=cellfun(@(c) c(:), varargin,'uni',0);
for ii=2:nv
A=kron(varargin{ii},A);
end
Lengths=cellfun(@(c) length(c),varargin);
A=reshape(A,Lengths);
A=ndSparse(A);
[A,B]=deal(B,A);
else
A=ndSparse(A);
end
ndimsA=ndims(A);
N=max(ndimsA,ndimsB);
[Acoords,Avals,Ashape]=getEntryTable(A,N);
[Bcoords,Bvals,Bshape]=getEntryTable(B,N);
convShape=Ashape+Bshape-1;
switch truncType
case 'full'
%Nothing to do
lb=[]; ub=[];
case 'same'
lb=floor(Bshape/2)+1;
ub=Ashape+lb-1;
case 'valid'
lb=Bshape;
ub=convShape+1-Bshape;
if any(lb>ub) %result will be empty
convShape=ub-lb+1;
convShape(convShape<0)=0;
obj=ndSparse.build([],[],convShape );
return
end
otherwise
error 'SHAPE must be ''full'', ''same'', or ''valid''.'
end
Bcoords=[Bcoords{:}];
Corner=cellfun(@(c) min(c), Acoords,'uni',0);
Acoords=[Acoords{:}];
S.type='()'; S.subs=Corner;
cornerval=full(subsref(A,S));
if ~isequal(cornerval,0)
Acoords=[[Corner{:}];Acoords];
Avals=[0;Avals];
end
convVals=Avals(:)*Bvals(:).';
convCoords=bsxfun(@plus, reshape(Bcoords.',N,1,[]) ,Acoords.'-1);
convCoords=reshape( convCoords,N,[]).';
if ~isempty(lb)%Only necessary for 'same' and 'valid'
idx = all( bsxfun(@ge,convCoords,lb) & bsxfun(@le,convCoords,ub) ,2);
convCoords=bsxfun(@minus, convCoords(idx,:) , lb-1);
convVals=convVals(idx);
convShape=ub-lb+1;
end
obj=ndSparse.build(convCoords,convVals,convShape);
end
end%NONSTATIC NON HIDDEN METHODS
methods (Hidden)
function data=getData(obj)
data=obj.data;
end
function ndShape=getndShape(obj)
ndShape=obj.ndShape;
end
function objnew=sumEngine(memliberal, obj,varargin)
%sumEngine
memliberal=memliberal | ndSparse.oldstyle;
if ndims(obj)<=2 %Matrix case
objnew=ndSparse( sum(obj.data, varargin{:}) );
return;
end
if isempty(varargin) || ischar(varargin{1})
dim=find(obj.ndShape,1,'first');
else
dim=varargin{1};
end
sz=size(obj);
nn=length(sz);
if dim>nn || sz(dim)==1,objnew=obj; return; end
if memliberal
data=ExtractDim1Reshaping(obj.data, dim, sz);
data=sum(data,1);
sz(dim)=1;
data=invExtractDim1Reshaping(data,dim,sz);
objnew=ndSparse(data,sz);
else
[ndSubs,vals]=getEntryTable(obj,nn);
ndSubs{dim}(:)=1;
sz(dim)=1;
objnew=ndSparse.accumarray([ndSubs{:}],vals,sz);
end
end
function objnew=anyEngine(memliberal,obj,varargin)
%anyEngine
memliberal=memliberal | ndSparse.oldstyle;
if ndims(obj)<=2 %Matrix case
objnew=ndSparse( any(obj.data, varargin{:}) );
return;
end
if isempty(varargin)
dim=find(obj.ndShape,1,'first');
else
dim=varargin{1};
end
sz=size(obj);
nn=length(sz);
if dim>nn || sz(dim)==1,
objnew=obj;
if ~islogical(objnew.data)
objnew.data=logical(objnew.data);
end
return;
end
if memliberal
data=ExtractDim1Reshaping(obj.data,dim,sz);
data=any(data,1);
sz(dim)=1;
data=invExtractDim1Reshaping(data,dim,sz);
objnew=ndSparse(data,sz);
else
objnew = sumEngine(0, obj~=0 ,dim)>0;
end
end
function objnew=allEngine(memliberal,obj,varargin)
%allEngine
memliberal=memliberal | ndSparse.oldstyle;
if ndims(obj)<=2 %Matrix case
objnew=ndSparse( all(obj.data, varargin{:}) );
return;
end
if isempty(varargin)
dim=find(obj.ndShape,1,'first');
else
dim=varargin{1};
end
sz=size(obj);
nn=length(sz);
if dim>nn || sz(dim)==1,
objnew=obj;
if ~islogical(objnew.data)
objnew.data=logical(objnew.data);
end
return;
end
if memliberal
data=ExtractDim1Reshaping(obj.data,dim,sz);
data=all(data,1);
sz(dim)=1;
data=invExtractDim1Reshaping(data,dim,sz);
objnew=ndSparse(data,sz);
else
objnew=( sumEngine(0, obj~=0 , dim) == sz(dim) );
end
end
function objnew=meanEngine(memliberal, obj,varargin)
%meanEngine
memliberal=memliberal | ndSparse.oldstyle;
if ndims(obj)<=2 %Matrix case
objnew=ndSparse( mean(obj.data, varargin{:}) );
return;
end
if isempty(varargin)
dim=find(obj.ndShape,1,'first');
else
dim=varargin{1};
end
sz=size(obj);
nn=length(sz);
if dim>nn || sz(dim)==1,objnew=obj; return; end
if memliberal
data=ExtractDim1Reshaping(obj.data,dim,sz);
data=mean(data,1);
sz(dim)=1;
data=invExtractDim1Reshaping(data,dim,sz);
objnew=ndSparse(data,sz);
else
objnew=sumEngine(0,obj,dim)/sz(dim);
end
end
function objnew=circshiftEngine(memliberal,obj,circDims)
%circshiftEngine
memliberal=memliberal | ndSparse.oldstyle;
ndShape=size(obj);
nn=length(ndShape);
if nn==2
obj.data=circshift(obj.data,circDims);
objnew=obj;
return;
end
N=max(length(circDims),nn);
circDims=[circDims,ones(1,N-length(circDims))];%BUG? should be zeros()?
sz=trail1s(ndShape,N);
initialshape=sz;
map=(circDims~=0);
firstShift=find(map,1,'first');
if isempty(firstShift) || firstShift>nn,
objnew=obj;return
end
if memliberal
%reformulate
sz=circshift(sz,[0,1-firstShift]);
circDims=circDims(firstShift:end);
data=ExtractDim1Reshaping(obj.data,firstShift,initialshape);
mm=length(circDims);
for ii=1:mm
data=reshape(data,sz(ii),[]);
data=data.';
thisShift=circDims(ii);
if thisShift
data=circshift(data,[0,thisShift]);
end
end
objnew=ndSparse(data, initialshape);
else%~memliberal
[ndSubs,vals]=getEntryTable(obj,nn);
if isempty(vals); objnew=obj; return; end %all zero array
circDims=circDims(1:nn);
for ii=1:nn
if circDims(ii)==0, continue; end
ndSubs{ii}=mod(ndSubs{ii}+circDims(ii)-1,sz(ii))+1;
end
objnew=ndSparse.build(ndSubs,vals,sz,nzmax(obj));
end
end
function objnew=catEngine(memliberal,dim,varargin)
%catEngine - same syntax as for full arrays
memliberal=memliberal | ndSparse.oldstyle;
Nargs=length(varargin);
map=cellfun('isclass',varargin,'ndSparse');
first=find(map,1);
if first==1
L=varargin{first};
if Nargs==1, objnew=L; return; end
R=varargin{first+1};
LastProcessed=2;
else
L=cat(dim,varargin{1:first-1});
R=varargin{first};
LastProcessed=first;
end
[szL,ndimsL]=trail1s(size(L),dim);
[szR,ndimsR]=trail1s(size(R),dim);
idx=true(1,length(szL));
idx(dim)=0;
if ndimsL~=ndimsR || any(szL(idx)~=szR(idx))
error 'CAT arguments dimensions are not consistent.'
else
newshape=szL;
newshape(dim)=newshape(dim)+szR(dim);
end
if dim>=ndimsL || (ndimsL<=2) %equivalent to 2D cat
L=mkCompat(L,szL);
R=mkCompat(R,szR);
objnew=ndSparse( cat(min(dim,2),L,R),newshape);
else %Fully Multi-dimensional case
if memliberal
L=mkCompat(L,szL);
R=mkCompat(R,szR);
L=ExtractDim1Reshaping(L,dim,szL);
R=ExtractDim1Reshaping(R,dim,szR);
L=[L.',R.'].';
L=invExtractDim1Reshaping(L,dim,newshape);
objnew=ndSparse(L,newshape);
else
if ~map(1), L=ndSparse(L); end
if ~map(2), R=ndSparse(R); end
[LSubs,Lvals]=getEntryTable(L,ndimsL);
[RSubs,Rvals]=getEntryTable(R,ndimsR);
RSubs{dim}=RSubs{dim}+szL(dim);
objnew=ndSparse.build( [[LSubs{:}];[RSubs{:}]] , [Lvals;Rvals] , newshape);
end
end
%Use recursion to handle further CAT args
if LastProcessed<Nargs
objnew=catEngine(memliberal,dim,objnew,varargin{LastProcessed+1:Nargs});
end
end
function varargout=minEngine(memliberal,L,varargin)
%minEngine
memliberal=memliberal | ndSparse.oldstyle;
if nargin>2,
R=varargin{1};
else
R=[];
end
if nargin>3
dim=varargin{2};
elseif isempty(R);
dim=find(L.ndShape,1,'first');
end
if isempty(R)
sz=size(L);
nn=length(sz);
if dim>nn || sz(dim)<=1,
varargout{1}=L;
if nargout>1,
varargout{2}=ones(size(L));
end
return;
end
if memliberal
data=ExtractDim1Reshaping(L.data,dim,sz);
[varargout{1:nargout}]=min(data,[],1);
sz(dim)=1;
varargout{1}=invExtractDim1Reshaping(varargout{1},dim,sz);
varargout{1}=ndSparse(varargout{1},sz);
if nargout>1
varargout{2}=invExtractDim1Reshaping(varargout{2},dim,sz);
varargout{2}=reshape(varargout{2},sz);
end
else
data=ExtractDimNReshaping(L.data,dim,sz);
[varargout{1:nargout}]=min(data,[],2);
sz(dim)=1;
varargout{1}=invExtractDimNReshaping(varargout{1},dim,sz);
if nargout>1
varargout{2}=invExtractDimNReshaping(varargout{2},dim,sz);
end
end
elseif length(varargin)>1 %R and dim given
error('MIN with two matrices to compare and a working dimension is not supported.')
elseif nargout>1%R given, dim not given, but 2nd argout requested
error('MIN with two matrices to compare and two output arguments is not supported.')
else%bi-operand min
s=getshape(L,R);
varargout{1} = finalObject( min(mkCompat(L,s),mkCompat(R,s)),s);
end
end
function varargout=maxEngine(memliberal,L,varargin)
%maxEngine
memliberal=memliberal | ndSparse.oldstyle;
if nargin>2,
R=varargin{1};
else
R=[];
end
if nargin>3
dim=varargin{2};
elseif isempty(R);
dim=find(L.ndShape,1,'first');
end
if isempty(R)
sz=size(L);
nn=length(sz);
if dim>nn || sz(dim)<=1,
varargout{1}=L;
if nargout>1,
varargout{2}=ones(size(L));
end
return;
end
if memliberal
data=ExtractDim1Reshaping(L.data,dim,sz);
[varargout{1:nargout}]=max(data,[],1);
sz(dim)=1;
varargout{1}=invExtractDim1Reshaping(varargout{1},dim,sz);
varargout{1}=ndSparse(varargout{1},sz);
if nargout>1
varargout{2}=invExtractDim1Reshaping(varargout{2},dim,sz);
varargout{2}=reshape(varargout{2},sz);
end
else
data=ExtractDimNReshaping(L.data,dim,sz);
[varargout{1:nargout}]=max(data,[],2);
sz(dim)=1;
varargout{1}=invExtractDimNReshaping(varargout{1},dim,sz);
if nargout>1
varargout{2}=invExtractDimNReshaping(varargout{2},dim,sz);
end
end
elseif length(varargin)>1 %R and dim given
error('MAX with two matrices to compare and a working dimension is not supported.')
elseif nargout>1%R given, dim not given, but 2nd argout requested
error('MAX with two matrices to compare and two output arguments is not supported.')
else%bi-operand max
s=getshape(L,R);
varargout{1} = finalObject( max(mkCompat(L,s),mkCompat(R,s)),s);
end
end
end%hidden methods
methods (Static)
function obj=loadobj(B)
obj=ndSparse(B.data,B.ndShape);
end
function obj=build(varargin)
%ndSparse.build - This is a generalization of sparse(i,j,s,m,n,nzmax)
%to N dimensions, allowing one to build an ndSparse array from a
%user-supplied table of explicit array entries.
%
%SYNTAX:
%
% S = ndSparse.build(explicitCoords,s,ndShape, nzmax)
%
%
%out:
%
% S: the N-dimensional sparse output array as an ndSparse object.
%
%in:
%
% explicitCoords: a PxN matrix whos rows are the N-dimensionsal
% coordinates of the explicit entries in the desired
% output array S.
%
% s: a P-vector of array values at the locations
% supplied in explicitCoords. Can also be a scalar.
%
% ndShape: The desired dimensions of S. Defaults to max(explicitCoords,[],1).
%
% nzmax: The output S will have memory allocated for nzmax
% non-zeros. Default = nnz(S).
%
%
%To create an ndSparse array of zeros, one can do
%
% S=ndSparse.build(ndShape)
%
%which is a simplifcation of S=ndSparse.build([],[],ndShape,0).
%
%See also ndSparse.spalloc, ndSparse.sprand, ndSparse.sprandn
nzmaxARG={};
switch nargin
case 1 %ndSparse.build(ndShape)
ndShape=varargin{1};
[M,N,ndShape]=ndShape2MN(ndShape);
data=sparse(M,N);
obj=ndSparse(data,ndShape);
return;
case 2 %ndSparse.build(explicitCoords,s)
explicitCoords=varargin{1};
s=varargin{2};
if ~iscell(explicitCoords)
ndShape=max(explicitCoords,[],1);
%nzmax=size(explicitCoords,1);
else
ndShape=cellfun(@max,explicitCoords);
%nzmax=length(explicitCoords{1});
end
case 3 %ndSparse.build(explicitCoords,s,ndShape)
[explicitCoords,s,ndShape]=deal(varargin{:});
% if ~iscell(explicitCoords)
% nzmax=size(explicitCoords,1);
% else
% nzmax=length(explicitCoords{1});
% end
case 4 %ndSparse.build(explicitCoords,s,ndShape,nzmax)
[explicitCoords,s,ndShape,nzmaxARG{1}]=deal(varargin{:});
otherwise
error 'Too many inputs'
end
if iscell(explicitCoords)
emptyCoords=any(cellfun('isempty',explicitCoords));
else
emptyCoords=isempty(explicitCoords);
end
if emptyCoords && isempty(s)
[M,N]=ndShape2MN(ndShape);
data=sparse([],[],[], M, N, nzmaxARG{:});
obj=ndSparse(data,ndShape);
return
end
[I,J,dims]=ndCoords2IJ(explicitCoords,ndShape);
S=s;
M=dims(1);
N=dims(2);
data=sparse(I,J,S,M,N,nzmaxARG{:});
obj=ndSparse(data,ndShape);
end
function obj=spalloc(ndShape,nzmax)
%spalloc method (Static) for ndSparse
%
%SYNTAX:
%
% obj=ndSparse.spalloc(ndShape,nzmax)
%
%in:
%
% ndShape: desired dimensions of the output ndSparse array obj
% nzmax : desired capacity for nonzeros
%
%out: output all-zero ndSparse object of dimensions ndShape
obj=ndSparse.build([],[],ndShape,nzmax);
end
function obj=sprand(ndShape,density)
%sprand method (Static) for ndSparse
%
%SYNTAX:
%
% obj=ndSparse.sprand(ndShape,density)
%
%in:
%
% ndShape: desired dimensions of the output ndSparse array obj
% density : desired density of nonzeros
%
%out: output ndSparse object of dimensions ndShape
[mm,nn,ndShape]=ndShape2MN(ndShape);
obj=ndSparse( sprand(mm,nn,density) , ndShape);
end
function obj=sprandn(ndShape,density)
%sprandn method (Static) for ndSparse
%
%SYNTAX:
%
% obj=ndSparse.sprand(ndShape,density)
%
%in:
%
% ndShape: desired dimensions of the output ndSparse array obj
% density : desired density of nonzeros
%
%out: output ndSparse object of dimensions ndShape
[mm,nn,ndShape]=ndShape2MN(ndShape);
obj=ndSparse( sprandn(mm,nn,density) , ndShape);
end
function obj=accumarray(subs,val,sz,fun)
%ndSparse.accumarray
%
%SYNTAX:
%
% obj = ndSparse.accumarray(subs,val,sz,fun)
%
%This works pretty much like MATLAB's built-in accumarray(subs,val,sz,fun)
%except that the output obj is always ndSparse. Defaults for sz and fun
%are the same as for the native accumarray as well.
%
%Note that the normal built-in accumarray also supports a 5th and 6th argument
%FILLVAL and an ISSPARSE, which are not supported here. They make no sense
%for an ndSparse output obj, since it will always be sparse, and hence also the only
%sensible FILLVAL is zero.
if nargin==4, do_fun=true; else do_fun=false; end
if iscell(subs)
msg='Cells of first input SUBS must contain real, full, numeric vectors of equal length.';
subs=cellfun(@(c) c(:), subs, 'uni', 0);
try
subs=[subs{:}];
catch
error(msg)
end
if ~isreal(subs) || issparse(subs) || ~isnumeric(subs)
error(msg)
end
end
N=size(subs,2);
if ~isreal(subs) || issparse(subs) || ~isnumeric(subs) || N<1
error 'First input SUBS must be a real, full, numeric matrix with at least one column.'
end
minsz=max(subs,[],1);
if isempty(minsz), minsz=zeros(1,N); end
if isscalar(minsz), minsz=[minsz,1]; end
if nargin<3 || isempty(sz),
sz=minsz;
end
if N>1 && length(sz)~=N
error 'Third input SZ must be a full row vector with one element for each column of SUBS.'
elseif N==1 && length(sz)<2
error 'When input SUBS is a column vector, input SZ must be of the form [M,1]'
end
if ~all(sz>=minsz)
error 'First input SUBS and third input SZ must satisfy ALL(MAX(SUBS)<=SZ).'
end
if isempty(subs)
obj=ndSparse.build(sz); return
end
desiredShape=sz;
[ii,jj,dims]=ndCoords2IJ(subs,desiredShape);
if ~isa(val,'double')
val=double(val);
end
if do_fun
data=accumarray([ii,jj],val,dims,fun,0,true);
else
data=accumarray([ii,jj],val,dims,[],0,true);
end
obj=ndSparse(data,desiredShape);
end
end%STATIC METHODS
end
function [E,targetshape,boolGrow,other]=nd2matrixIndex(ndSubs,ndShape,operation)
%nd2matrixIndex
%
% ndShape is expected to be untrailed by ones
N=length(ndSubs);
[ndShape, numdims]=untrail1s(ndShape);
[objIs1D,obj1Ddim]=proc1D(ndShape);
LinearIndexing=(N==1);
boolConsol = 1<N && N<numdims; %consolidating indexing of trailing dims.
if ~LinearIndexing%subscript indexing
ShapeCell=parseSize(ndShape,N);
quasiShape=[ShapeCell{:}];
else%linear indexing
idx1=ndSubs{1};
quasiShape=prod(ndShape);
end
ndSubs0=ndSubs;%initial substruct
ndSubs=cellfun(@mkCompat, ndSubs,'uni',0);
targetshape=zeros(1,N);
accessedshape=zeros(1,N);
subshape=zeros(1,N);
boolGrow=false;
for ii=1:N %loop over indices
idx=ndSubs{ii};
if islogical(idx) %logical index
idx=find(idx);
maxAccessed=max(idx);
if isempty(maxAccessed), maxAccessed=quasiShape(ii); end
subshape(ii)=numel(idx);
elseif ischar(idx) && ~LinearIndexing %colon index
idx=(1:quasiShape(ii)).';
maxAccessed=quasiShape(ii);
subshape(ii)=maxAccessed;
elseif ischar(idx) && LinearIndexing %colon index
idx=':';
maxAccessed=quasiShape(ii);
subshape(ii)=maxAccessed;
else %numeric index
maxAccessed=max(idx);
subshape(ii)=numel(idx);
end
switch operation
case 'subsref'
if maxAccessed>quasiShape(ii)
error 'Index exceeds matrix dimensions.'
end
targetshape(ii)=subshape(ii);
accessedshape(ii)=quasiShape(ii);
case 'subsasgn'
accessedshape(ii)=max(maxAccessed,quasiShape(ii));
targetshape(ii)=accessedshape(ii);
if maxAccessed>quasiShape(ii)
boolGrow=true;
end
otherwise
error 'Unknown operation'
end
ndSubs{ii}=idx;
end
if ~LinearIndexing, %Clean up trailing dimensions
[accessedshape,N]=untrail1s(accessedshape);
targetshape=targetshape(1:N);
subshape=subshape(1:N);
ndSubs=ndSubs(1:N);
end
if boolGrow
if LinearIndexing && ~objIs1D %linear indexing of non-vector
if numdims==2
error 'In an assignment A(I) = B, a matrix A cannot be resized.'
elseif ndims>2
error 'Attempt to grow array along ambiguous dimension.'
end
elseif boolConsol
error 'Attempt to grow array along ambiguous dimension.'
end
elseif isequal(operation,'subsasgn')
targetshape=ndShape;
end
if LinearIndexing && isequal(operation,'subsref') %linear indexing in subsref
idx1=ndSubs{1};
if ~ischar(idx1)
idxShape=size(idx1);
else%single colon
idxShape=[quasiShape(1),1];
end
idxIs1D=proc1D(idxShape);
if idxIs1D && objIs1D
tmp=ndShape;
tmp(obj1Ddim)=targetshape;
targetshape=tmp;
else % shape of idx1 dictates resulting shape
targetshape=idxShape;
end
end
%%Everything up to here was to convert indices to numeric form
%%and to compute targetshape and accessedshape - now prepare
%%subscripts
E.type='()';
E.subs{1}=ndSubs{1};
multicolons=false;
switch N
case 1 %linear indexing, use index as is.
case 2
E.subs{2}=ndSubs{2};
otherwise%N>2
E.subs{2}=ndSubs{N};
%=numdims;
nnn=N;
if cellfun('isclass', ndSubs0(1:nnn-1),'char')
%Colons in the first numdims-1 subscripts
multicolons=true;
E.subs{1}=':';
else
E.subs{1}=sub2allind(accessedshape(1:N-1),ndSubs{1:N-1});
end
end
other.subshapeND=subshape;
other.subshape=cellfun('length',E.subs);
if LinearIndexing,
other.subshape=subshape;
elseif multicolons
other.subshape(1)=prod(ndShape(1:N-1));
end
other.boolConsol=boolConsol;
other.quasiShape=quasiShape;
other.LinearIndexing=LinearIndexing;
end
function outcell=parseSize(sz,numargsout,dim)
%parseSize
idx=find(sz~=1,1,'last');
if isempty(idx)
sz=[1,1];
elseif idx==1
sz=[sz(1), 1];
else
sz=sz(1:idx);
end
dimSpecified=logical(exist('dim','var'));
if dimSpecified,
if ~isscalar(dim), error('DIM argument must be scalar.'); end
if dim>length(sz), sz(end+1:dim)=1; end
sz=sz(dim);
end
if numargsout<=1,
outcell{1}=sz;
return;
elseif dimSpecified,
error('Too many output arguments.')
else
if length(sz)>=numargsout
sz(numargsout)=prod(sz(numargsout:end));
sz(numargsout+1:end)=[];
else
sz(end+1:numargsout)=1;
end
outcell=num2cell(sz);
end
end
function X=mkCompat(X,s)
%mkCompat
if isa(X,'ndSparse')
X=X.data;
elseif ~(isa(X,'double') || islogical(X) || ischar(X))
%char is for handling the case where X is the index ':'
X=double(X);
end
if nargin>1 && ~isscalar(X)
X=reshape(X,[],s(end));
end
end
function out=finalObject(obj, ndShape)
%finalObject
if nargin<2,
out=ndSparse(obj);
elseif ~issparse(obj)%a full object type is produced (not ndSparse or native sparse)
out=reshape(obj,ndShape);
else
out=ndSparse(obj,ndShape);
end
end
function shape=getshape(L,R)
%getshape
szL=size(L);
szR=size(R);
if isscalar(L)
shape=szR;
elseif isscalar(R)
shape=szL;
elseif ~isequal(szL,szR)
error 'Array dimensions must match for binary array op.'
else
shape=szL;
end
end
function data=ExtractDim1Reshaping(data,dim,ndShape)
%ExtractDim1Reshaping
%ndShape=obj.ndShape;
%data=obj.data;
nd=length(ndShape);
if dim==1
data=reshape(data,ndShape(1),[]);
elseif dim==nd
data=data.';
elseif dim<nd %it is an implied condition here that dim>1
data=reshape(data, prod(ndShape(1:dim-1)),[]);
data=reshape(data.',ndShape(dim),[]);
else% dim>nd
%ndShape=trail1s(ndShape,dim);
error 'DIM>nD. This situation should be handled by separate processing'
end
end
function data=invExtractDim1Reshaping(data,dim,newshape)
%invExtractDim1Reshaping - not a precise inverse, only up to a reshaping
nd=length(newshape);
if dim==1
elseif dim==nd
data=data.';
elseif dim<nd %it is an implied condition here that dim>1
data=reshape( data, [], prod( newshape(1:dim-1) ) ).' ;
else% dim>nd
error 'DIM>nD. This situation should be handled by separate processing'
end
%obj=ndSparse(data,newshape);
end
function data=ExtractDimNReshaping(data,dim,ndShape)
%ExtractDimNReshaping
nd=length(ndShape);
if dim==nd
%do nothing
elseif dim<nd
order=[1:dim-1, dim+1:nd, dim];
obj=ndSparse(data,ndShape);
obj=permute(obj,order);
data=obj.data;
if ndShape(dim)==1
data=data(:);
end
else% dim>nd
error 'DIM>nD. This situation should be handled by separate processing'
end
end
function obj=invExtractDimNReshaping(data,dim,newshape)
%invExtractDimNReshaping
nd=length(newshape);
if dim<=nd
order=[1:dim-1, dim+1:nd, dim];
obj=ndSparse(data,newshape(order));
obj=ipermute(obj,order);
else% dim>nd
error 'DIM>nD. This situation should be handled by separate processing'
end
end
function [ndShape,N]=trail1s(ndShape,dim)
%trail1s
N=length(ndShape);
if N==1,
ndShape=[ndShape,1];
N=2;
end
if dim>N, %add trailing ones
ndShape=[ndShape, ones(1,dim-N)];
N=dim;
end
end
function [ndShape,N]=untrail1s(ndShape)
%untrail1s
idx=find(ndShape~=1,1,'last');
if isempty(idx)
ndShape=[1,1];
elseif idx==1
ndShape=[ndShape(1),1];
else
ndShape=ndShape(1:idx);
end
N=length(ndShape);
end
function [MM,NN,ndShape]=ndShape2MN(ndShape)
%ndShape2MN
ndShape=untrail1s(ndShape);
MM=prod(ndShape(1:end-1));
NN=ndShape(end);
end
function [ndSubs,vals,ndShape]=getEntryTable(obj,maxdim)
%getEntryTable
%
%[ndSubs,vals,ndShape]=getEntryTable(obj,maxdim)
%
%Produces a cell array ndSubs{1:maxdim} of N-dimensional
%subscripts of explicit entries of obj and a corresponding array,
%val, of values. The output NDSHAPE will be the dimension of obj with any
%required trailing 1 dimensions.
N=maxdim;
ndShape=trail1s(obj.ndShape,N);
[II,JJ,vals]=find(obj.data);
ndSubs=IJ2ndCoords(II,JJ,ndShape);
vals=vals(:);
end
function [II,JJ,dims]=ndCoords2IJ(ndCoords,ndShape)
%ndCoords2IJ
if ~iscell(ndCoords), ndCoords=num2cell(ndCoords,1); end
if any(cellfun('isempty',ndCoords))
II=[]; JJ=[];
dims=[ prod( ndShape(1:end-1) ) , ndShape(end)];
return
end
[ndShape,N]=trail1s(ndShape,length(ndCoords));
minShape=trail1s( cellfun(@max, ndCoords) , N );
if any(minShape>ndShape)
error 'Given array ndShape is not large enough to hold given n-D coordinates.'
end
[ndShape,numdims]=untrail1s(ndShape);
if ndShape(2)==1 && numdims<=2,
numdims=1;
end
ndCoords=ndCoords(1:numdims);
if numdims==1
II=ndCoords{1}(:);
JJ=ones(length(II),1);
dims=ndShape;
elseif numdims==2
II=ndCoords{1}(:);
JJ=ndCoords{2}(:);
dims=ndShape;
else
II=sub2ind(ndShape(1:end-1),ndCoords{1:end-1});
II=II(:);
JJ=ndCoords{end}(:);
dims=[ prod( ndShape(1:end-1) ) , ndShape(end)];
end
end
function [ndCoords,ndShape]=IJ2ndCoords(II,JJ,ndShape)
%IJ2ndCoords
N0=length(ndShape);
[ndShape,N]=untrail1s(ndShape);
ndCoords=[ cell( 1, N-1 ) , {JJ} ];
[ndCoords{1:N-1}]=ind2sub(ndShape(1:end-1), II);
ndCoords=cellfun(@(c) c(:), ndCoords, 'uni',0);
if N0>N
ndCoords=[ndCoords, num2cell( ones(size(ndCoords{1},1),N0-N) , 1) ];
end
end
function bool=iscolumn(A)
sz=size(A);
bool=length(sz)==2 && sz(2)==1;
end
function bool=isrow(A)
sz=size(A);
bool=length(sz)==2 && sz(1)==1;
end
%%%%%%%%%%%%%%%%%%%%% EXTERNAL CONTRIBUTIONS %%%%%%%%%%%%%%%%%%%%%%%%
function idx = sub2allind(sz, varargin)
% idx = sub2allind(sz, sub1, sub2, sub3, ... )
%
% Like sub2ind, sub2allind computes the equivalent linear indices for
% given subscripts of an array with size SZ.
% Unlike sub2ind, it computes them for all combinations of
% subscripts. So, instead of calling A( 2:3, 1, 4:11) you might
% use
% linIdx = sub2allind( size(A), 2:3, 1, 4:11 );
%
% and then call A(linIdx).
%
% Using the colon operator is allowed:
% linIdx = sub2allind( sz, sub1, :, sub3 );
% Michael Völker, 2011
% This is Matlab FileEx #30096
% ==================================================================
% Make sure, "sz" is fine
%
if nargin < 2
error('sub2allind:NoOfInputs', 'At least two inputs, please.')
end
if ~exist('sz', 'var') || length(sz(:)) < 2 || ~all(isnumeric(sz(:))) || ~all(isreal(sz(:)))
error('sub2allind:szNotSane', 'Size vector must be real and have at least 2 elements.');
end
sz = sz(:).';
Ndims = length(sz);
if Ndims < nargin-1
ResDims = nargin-1 - Ndims;
sz = [sz ones(1,ResDims)]; % Adjust for trailing singleton dimensions
elseif Ndims > nargin-1
sz = [sz(1:nargin-2) prod(sz(nargin-1:end))]; % Adjust for linear indexing on last element
end
% ==================================================================
Ndims = length(sz); % No. of dimensions
% ==================================================================
% Make sure all subscripts are fine, including the use of ":"
%
for d = 1:Ndims
flagNum = isnumeric( varargin{d}(:) );
flagReal = isreal( varargin{d}(:) );
flagInt = all( varargin{d}(:) == floor(varargin{d}(:)) );
if ~flagNum || ~flagReal || ~flagInt
if isequal( varargin{d}, ':' ) % We allow exactly one type of non-numeric input
varargin{d} = 1:sz(d); % interpret ":" as in usual subscript syntax
else
error('sub2allind:SubsNotSane', 'Subscripts must either be valid matrix subscripts or the well known '':''.')
end
end
end
% ==================================================================
% ==================================================================
% Compute linear indices, e.g. in 3D, with edges L1, L2, L3:
%
% idx = x1 + (x2-1)*L1 + (x3-1)*L1*L2,
%
% for every permutation (x1,x2,x3)
k = [1 cumprod(sz(1:end-1))]; % k(d) holds the accumulated No. of array elements
% up to the d'th subscript (or dimension)
idx = 1; % smallest possible index
for d = 1:Ndims
xd = varargin{d}(:); % the d'th subscripts
if any(xd < 1) || any(xd > sz(d)) % Verify subscripts are within range
error('sub2allind:IndexOutOfRange', 'Out of range subscript.');
end
reshrule = ones(1,d+1); % how to reshape the size of xd
reshrule(end) = length(xd); %
xd = reshape( xd, reshrule ); % prepare for bsxfun's needs
idx = bsxfun( @plus, idx, (xd-1)*k(d) ); % iteratively calculate the sum pointed out above
end
% ==================================================================
idx = idx(:); % linearify indices
end
function [bool,dim] = proc1D(sz)
dim=find(sz>1);
ldim=length(dim);
map=double(sz>0);
bool= (ldim<=1) & all(map);
if ~bool
dim=[];
elseif ldim==0 %means the object is a scalar
dim=1;
end
end