Spot: Subset reference and assignment nearly complete

git-svn-id: svn://tuatara.cs.ubc.ca/spot/trunk@15 602325e8-326d-0410-a6db-8800ab266cdc

@opSpot/mtimes.m

@@ -58,9 +58,15 @@
       y = zeros(m,q);
 
       % Perform operator*vector on each column of B
-      for i=1:q
-          y(:,i) = A.apply(B(:,i),1);
+      [m,n,p,q]
+      if isempty(A) % Zero rows or columns
+         % Nothing to be done, result will be y
+      else
+         for i=1:q
+            y(:,i) = A.apply(B(:,i),1);
+         end
       end
+
    end   
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Both args are Spot ops.

@opSpot/opSpot.m

@@ -52,6 +52,9 @@
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     methods (Access = protected)
        function y = apply(op,x,mode)
+          % Update counter
+          op.counter(mode) = op.counter(mode) + 1;
+
           y = op.multiply(x,mode);
        end
     end % methods - private
@@ -60,7 +63,7 @@
     % Abstract methods -- must be implemented by subclass.
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     methods (Abstract, Access = protected)
-        res = multiply(op,x,mode)
+        y = multiply(op,x,mode)
     end % methods - abstract
 
 end % classdef

@opSpot/subsasgn.m

@@ -0,0 +1,129 @@
+function result = subsasgn(op,s,b)
+%   Copyright 2009, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: subsasgn.m 39 2009-06-12 20:59:05Z ewout78 $
+
+switch s.type
+   case {'.'}
+      % Set properties and flags
+      op.(s.subs) = b;
+      result = op;
+
+   case {'{}'}
+      error('Cell-index access is read-only.');
+
+   case {'()'}
+
+      % Get operator size
+      [m,n] = size(op);
+
+      % Do not allow one-dimensional indexing
+      if length(s.subs) == 1
+         error('Vectorized subset assignments are not supported.');
+      end
+
+      % Ensure numeric vectors are full and in vectorized form
+      for i=1:length(s.subs)
+        idx = s.subs{i}; idx = idx(:);
+        if issparse(idx), idx = full(idx); end;
+        s.subs{i} = idx;
+      end
+
+      % Check dimensions (required when more than two dimensions are
+      % given)
+      dims = ones(1,length(s.subs));
+      dims(1) = Inf; % No bound (can grow)
+      dims(2) = Inf; % No bound (can grow)
+      for i=1:length(s.subs)
+         idx = s.subs{i};
+
+         if isempty(idx)
+            % Fine
+         elseif strcmp(idx,':')
+            % Fine
+         elseif islogical(idx)
+            if (length(idx) > dims(i))
+               error('Index exceeds operator dimensions.');
+            end
+         elseif isposintmat(idx)
+            if (max(idx) > dims(i))
+                error('Index exceeds operator dimensions.');
+            end
+         else
+            error(['Subscript indices must either be real positive ' ...
+                   'integers or logicals.']);
+         end
+      end
+
+      % Reset first and second dimensions
+      dims(1) = m; dims(2) = n;
+
+      % Check if all indices in given dimension have been specified
+      allIndex   = zeros(2,1);
+      emptyIndex = zeros(2,1);
+      sizeIndex  = zeros(2,1);
+      for i=1:2
+         idx = s.subs{i};
+
+         if strcmp(idx,':')
+            allIndex(i) = 1;
+            sizeIndex(i)= dims(i);
+         elseif islogical(idx)
+            if (length(idx) == dims(i) && all(idx))
+               allIndex(i) = 1;
+            elseif all(idx == false)
+               emptyIndex(i) = 1;
+            end
+            sizeIndex(i) = length(idx);
+         elseif isnumeric(idx)
+            if (max([idx;0]) == dims(i) && length(unique(idx)) == dims(i))
+               allIndex(i) = 1;
+            elseif isempty(idx)
+               emptyIndex(i) = 1;
+            end
+            sizeIndex(i) = length(idx);
+         end
+      end
+
+
+      % ---------------------------------------------------------------
+      % Handle special case where B = []; this will delete rows or
+      % columns from the operator
+      % ---------------------------------------------------------------
+      if ((size(b,1) == 0) && (size(b,2) == 0))
+         % Both dimensions fully specified
+         if allIndex(1) && allIndex(2)
+            % Empty Sparco operator
+            result = opEmpty(0,0);
+         elseif allIndex(1)
+            % Excise columns
+            result = opExcise(s.subs{2},'Cols');
+         elseif allIndex(2)
+            % Excise rows
+            result = opExcise(s.subs{1},'Rows');
+         elseif emptyIndex(1) && emptyIndex(2)
+            % Assign empty matrix to emtpy subset
+            result= op;
+         else
+            error('Subscripted assignment dimension mismatch.');
+         end
+
+         return;
+      end % B = []
+
+      % ---------------------------------------------------------------
+      % Check if size of indices and b match
+      % ---------------------------------------------------------------
+      if ~(size(b,1) == sizeIndex(1)) || ~(size(b,2) == sizeIndex(2))
+         error('Subscripted assignment dimension mismatch.');
+      end
+
+      % ---------------------------------------------------------------
+      % Handle the general case
+      % ---------------------------------------------------------------
+      if isempty(b)
+         result = op;
+      else
+         result = opSubsAsgn(op,s.subs{1},s.subs{2},b);
+      end
+end

@opSpot/subsref.m

@@ -0,0 +1,195 @@
+function varargout = subsref(op,s)
+%   Copyright 2009, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: subsref.m 39 2009-06-12 20:59:05Z ewout78 $
+
+if length(s) > 1
+   result = op;
+   for i=1:length(s)
+      if iscell(result)
+         if strcmp(s(i).type,'{}')
+            result = builtin('subsref',result,s(i));
+         else
+            % Apply the subsref to each element
+            newresult = cell(1,length(result));
+            for j=1:length(result)
+               newresult{j} = subsref(result{j},s(i));
+            end
+            result = newresult;
+         end
+      else
+         result = subsref(result,s(i));
+      end
+   end
+
+   if nargout > 1
+      for i=2:nargout
+         varargout{i} = [];
+      end
+   end
+   varargout{1} = result;
+
+   return;
+end
+
+
+switch s.type
+   case {'.'}
+      % Set properties and flags
+      varargout{1} = op.(s.subs);
+
+   case {'{}'}
+      error('Cell-indexing is not yet supported.');
+
+   case {'()'}
+      % Get operator size
+      [m,n] = size(op);
+
+      % Ensure numeric vectors are full
+      for i=1:length(s.subs)
+        if issparse(s.subs{i}), s.subs{i} = full(s.subs{i}); end;
+      end
+
+
+      % Process query
+      if (length(s.subs) == 1)
+         % --------------------------------------------------
+         % One-dimensional indexing; return explicit matrix
+         % --------------------------------------------------
+         idx   = s.subs{1};
+         [p,q] = size(idx);
+         idx   = idx(:);
+
+         if isempty(idx)
+            result = [];
+
+         elseif strcmp(idx,':')
+            % Get all entries
+            result = double(op);
+            result = result(:);
+
+         elseif isposintmat(idx)
+            % Check index range
+            if any(idx > m*n)
+               error('Index out of bounds.');
+            end
+
+            % Create matrix version of operator, restricted to the
+            % rows or columns accessed (this reduces the number of
+            % operator applications needed to generate the explicit
+            % matrix)
+            colIdx  = floor((idx+m-0.5)/m);
+            rowIdx  = mod(idx-1,m) + 1;
+            columns = unique(colIdx);
+            rows    = unique(rowIdx);
+
+            if length(rows) < length(columns)
+               rowMap =  zeros(max(rows),1);
+               rowMap(rows) = 0:length(rows)-1;
+
+               % Prepare ssub structure (cannot use (rows,:) in call)
+               ssub = struct();
+               ssub.subs = {rows,':'};
+               ssub.type = '()';
+
+               matrix  = double(subsref(op,ssub)');
+               matrix  = matrix(:);
+               result  = reshape(matrix(colIdx + rowMap(rowIdx)*n), p,q);
+            else
+               colMap = zeros(max(columns),1);
+               colMap(columns) = 0:length(columns)-1;
+
+               % Prepare ssub structure (cannot use (:,columns))
+               ssub = struct();
+               ssub.subs = {':',columns};
+               ssub.type = '()';
+
+               matrix  = double(subsref(op,ssub));
+               matrix  = matrix(:);
+               result  = reshape(matrix(rowIdx + colMap(colIdx)*m), p,q);
+            end
+
+         elseif islogical(idx)
+            % Check index range
+            if (numel(idx) > m*n) ...
+               || ((min(size(idx)) > 1) && (size(idx,1) > m || ...
+                                            size(idx,2) > n))
+               error('Index out of bounds.');
+            end
+
+            % Add padding to index and reshape
+            idx = idx(:);
+            if mod(length(idx),m) ~= 0
+               idx = [idx; false(m-mod(length(idx),m),1)];
+            end            
+            idx = reshape(idx,m,round(length(idx)/m));
+
+            % Create matrix version of operator, restricted to the
+            % rows or columns accessed (this reduces the number of
+            % operator applications needed to generate the explicit
+            % matrix)
+            colIdx = find(any(idx,1));
+            rowIdx = find(any(idx,2));
+
+            if length(rowIdx) < length(colIdx)
+               idx = idx(rowIdx,:);
+
+               % Prepare ssub structure (cannot use (:,colIdx))
+               ssub = struct();
+               ssub.subs = {rowIdx,':'};
+               ssub.type = '()';
+
+               matrix  = double(subsref(op,ssub)')';
+               matrix  = matrix(:);
+               result  = matrix(idx);
+            else
+               idx = idx(:,colIdx);
+
+               % Prepare ssub structure (cannot use (:,colIdx))
+               ssub = struct();
+               ssub.subs = {':',colIdx};
+               ssub.type = '()';
+
+               matrix  = double(subsref(op,ssub));
+               matrix  = matrix(:);
+               result  = matrix(idx);
+            end
+
+            if p == 1, result = result'; end;            
+         else
+            error('Invalid data type used for indexing.');
+         end
+
+         varargout{1} = result; % Numerical result
+      else
+         % --------------------------------------------------
+         % Higher-dimensional indexing -- create sub-operator
+         % --------------------------------------------------
+
+         % Check dimensions
+         dims = ones(length(s.subs),1);
+         dims(1) = m;
+         dims(2) = n;
+         for i=1:length(s.subs)
+             idx = s.subs{i}; idx = idx(:); s.subs{i} = idx;
+             if isempty(idx)
+                % Fine
+             elseif strcmp(idx,':')
+                % Fine
+             elseif islogical(idx)
+                if (length(idx) > dims(i))
+                   error('Index exceeds operator dimensions.');
+                end
+             elseif isposintmat(idx)
+                if (max(idx) > dims(i))
+                   error('Index exceeds operator dimensions.');
+                end
+             else
+                error(['Subscript indices must either be real positive ' ...
+                       'integers or logicals.']);
+             end
+         end
+
+         varargout{1} = opSubsRef(op,s.subs{1},s.subs{2});
+      end
+end