classdef (Hidden) SamplingGrid
   % Manages multi-dimensional sampled data.
   
   %   Author(s): P. Gahinet
   %   Copyright 2010-2013 The MathWorks, Inc.
   properties (Access = protected)
      % Grid data (struct)
      Data_
      % Version
      Version_ = ltipack.ver();
   end
   
   
   methods
      
      % Constructor
      function G = SamplingGrid(D)
         % Constructs sampling grid from structure or table
         if isstruct(D)
            % Struct input
            f = fieldnames(D);
            if ~isscalar(D) || isempty(f)
               error(message('Control:general:SamplingGrid1'))
            end
            % Check data type and consistency
            fs = size(D.(f{1}));
            for ct=1:numel(f)
               fV = D.(f{ct});
               if ~(isnumeric(fV) || islogical(fV))
                  error(message('Control:general:SamplingGrid2'))
               elseif ~isequal(size(fV),fs)
                  error(message('Control:general:SamplingGrid3',f{1},f{ct}))
               end
            end
            G.Data_ = D;
         elseif isa(D,'table') || isa(D,'dataset')
            % Table input
         else
            error(message('Control:general:SamplingGrid5'))
         end
      end
      
      function disp(G)
         % Display method
         f = fieldnames(G.Data_);
         Dims = sprintf('%dx',getSize(G));
         Vars = sprintf('%s, ',f{:});
         M = message('Control:general:SamplingGrid4',Dims(1:end-1),Vars(1:end-2));
         fprintf('%s\n',getString(M))
      end
      
      %% Data access
      function s = getSize(G)
         % Returns grid size
         D = G.Data_;
         f = fieldnames(D);
         s = size(D.(f{1}));
      end
      
      function D = getData(G)
         % Returns sampling grid as a struct
         D = G.Data_;
      end
      
      function C = getCellData(G)
         % Returns sampling grid as a cell array
         C = struct2cell(G.Data_);
      end
      
      function V = getVariable(G)
         % Returns list of sampling grid variables
         V = fieldnames(G.Data_);
      end
      
      %% Subselection/assignment
      function G = subselect(G,indices)
         % Selects portion of the grid
         D = G.Data_;
         f = fieldnames(D);
         for ct=1:numel(f)
            D.(f{ct}) = D.(f{ct})(indices{:});
         end
         G.Data_ = D;
      end
      
      function G = subdelete(G,indices)
         % Deletes portion of the grid
         D = G.Data_;
         f = fieldnames(D);
         for ct=1:numel(f)
            D.(f{ct})(indices{:}) = [];
         end
         G.Data_ = D;
      end
      
      function Glhs = reassign(Glhs,indices,Grhs,FinalArraySize)
         % Overwrites a portion of the grid
         % Determine if new grid is fully specified
         MaskLHS = false(FinalArraySize);
         if ~isequal(Glhs,[])
            SizeLHS = getSize(Glhs);
            nd = numel(SizeLHS);
            idx = cell(1,nd);
            for ct=1:nd
               idx{ct} = 1:SizeLHS(ct);
            end
            MaskLHS(idx{:}) = true;
         end
         MaskRHS = false(FinalArraySize);
         if ~isequal(Grhs,[])
            MaskRHS(indices{:}) = true;
         end
         Mask = MaskLHS | MaskRHS;
         if ~all(Mask(:))
            % Reassignment leads to partially specified grid
            warning(message('Control:general:SamplingGrid7'))
            Glhs = [];
         elseif ~isequal(Grhs,[])
            % RHS grid specified
            if isequal(Glhs,[])
               % full reassignment -> use RHS grid as initial template
               Glhs = Grhs;
            end
            % Check field compatibility
            Dlhs = Glhs.Data_;
            Drhs = Grhs.Data_;
            f = fieldnames(Dlhs);
            if isequal(sort(f),sort(fieldnames(Drhs)))
               % Proceed with overwrite
               for ct=1:numel(f)
                  Dlhs.(f{ct})(indices{:}) = Drhs.(f{ct});
               end
               Glhs.Data_ = Dlhs;
            else
               % Incompatible grids
               warning(message('Control:general:SamplingGrid6'))
               Glhs = [];
            end
         end
      end
      
      %% Grid manipulations
      function G = reshape(G,newSize)
         % Reshape grid
         D = G.Data_;
         f = fieldnames(D);
         try
            for ct=1:numel(f)
               D.(f{ct}) = reshape(D.(f{ct}),newSize);
            end
         catch ME
            throw(ME)
         end
         G.Data_ = D;
      end
      
      function G = repmat(G,Pattern)
         % Replicate and tile grid
         D = G.Data_;
         f = fieldnames(D);
         for ct=1:numel(f)
            D.(f{ct}) = repmat(D.(f{ct}),Pattern);
         end
         G.Data_ = D;
      end
      
      function G = permute(G,Order)
         % Permute grid dimensions
         D = G.Data_;
         f = fieldnames(D);
         for ct=1:numel(f)
            D.(f{ct}) = permute(D.(f{ct}),Order);
         end
         G.Data_ = D;
      end
      
      function G = squeeze(G)
         % Eliminate singleton dimensions
         D = G.Data_;
         f = fieldnames(D);
         for ct=1:numel(f)
            D.(f{ct}) = squeeze(D.(f{ct}));
         end
         G.Data_ = D;
      end
      
      function G1 = cat(dim,G1,G2)
         % Concatenates grid data
         D1 = G1.Data_;
         D2 = G2.Data_;
         f = fieldnames(D1);
         if isequal(sort(f),sort(fieldnames(D2)))
            for ct=1:numel(f)
               D1.(f{ct}) = cat(dim,D1.(f{ct}),D2.(f{ct}));
            end
            G1.Data_ = D1;
         else
            % Incompatible fields: clear grid
            warning(message('Control:general:SamplingGrid8'))
            G1 = [];
         end
      end
      
   end
   
   methods (Static)
      
      function G1 = merge(G1,varargin)
         % Grid consolidation in BINOP(M1,M2) where M1 and M2 are two
         % model arrays of identical size
         for j=1:numel(varargin)
            Gj = varargin{j};
            if isempty(G1)
               G1 = Gj;
            elseif ~isempty(Gj)
               % Both grids are specified. Use binop(f(X),g(Y)) model to merge grids
               D1 = G1.Data_;
               Dj = Gj.Data_;
               f1 = fieldnames(D1);
               fj = fieldnames(Dj);
               % Check consistency of common fields
               f = intersect(f1,fj);
               for ct=1:numel(f)
                  if ~isequal(D1.(f{ct}),Dj.(f{ct}))
                     warning(message('Control:general:SamplingGrid9',f{ct}))
                     G1 = [];
                     return
                  end
               end
               % Combine grids
               C = [struct2cell(D1) ; struct2cell(Dj)];
               [f,iu] = unique([f1 ; fj],'stable');
               G1.Data_ = cell2struct(C(iu,:),f,1);
            end
         end
      end
      
      function [GridInfo,RG] = getGridStructure(SG)
         % Takes an arbitrary sampling grid SG (struct) and tries to
         % decompose it into a multi-dimensional rectangular grid RG.
         % The struct GRIDINFO has the following fields:
         %    GridSize        Rectangular grid size
         %    GridVectors     Variable(s) and data associated with each grid
         %                    dimension. Cell array with as many entries as
         %                    grid dimensions, each entry being a two-column
         %                    cell array listing the variable names and grid
         %                    vectors associated with this dimension.
         %    SamplePerm      sample permutation mapping SG to RG
         %    Duplicate       true if SG has duplicated entries
         % The sampling grid RG is the sorted/reshaped version of SG that
         % would be obtained by applying NDGRID to the grid vectors.
         %
         % This function handles the following situations:
         %   * The dimensionality of SG need not match the grid dimension
         %   * The samples in SG are in arbitrary order and may be repeated
         %   * There may be more than one variable associated with a given grid
         %     dimension, and these variables are not necessarily contiguous
         %     fields of SG
         
         % Note: Algorithm leverages periodicity of lexicographical sorting for
         % regular grids
         f = fieldnames(SG);
         c = struct2cell(SG);
         nv = numel(f);
         for ct=1:nv
            c{ct} = c{ct}(:);
         end
         vals = cat(2,c{:});
         ns = size(vals,1);
         
         % Handle degenerate cases
         if ns<2
            % Treat each variable as a dimension
            GridVectors = cell(nv,1);
            if ns==0
               for ct=1:nv
                  GridVectors{ct} = {f{ct} []};
               end
               GridSize = zeros(1,max(nv,2));
            else
               for ct=1:nv
                  GridVectors{ct} = {f{ct} vals(ct)};
               end
               GridSize = ones(1,max(nv,2));
            end
            GridInfo = struct('GridSize',GridSize,'GridVectors',{GridVectors},...
               'SamplePerm',zeros(0,1),'Duplicate',false);
            RG = SG; 
            return
         end
         
         % Sort grid points
         [~,SamplePerm] = sortrows(fliplr(vals));
         vals = vals(SamplePerm,:);
         
         % Compute periodicity of value changes in a given column. In rectangular
         % grids, the period increases with the grid dimension. Handle case
         % when several variables are associated with same grid dimension.
         dv = [diff(vals,[],1) ; ones(1,nv)];
         for ct=nv:-1:1
            period(1,ct) = find(dv(:,ct)~=0,1);
         end
         Duplicate = any(all(dv==0,2));
         if Duplicate
            ndim = 1;  % Can't be rectangular grid, treat as 1D
         else
            [period,~,ju] = unique(period);   % sorted by increasing period
            nextperiod = [period(2:end),ns];  % p(k+1)
            ndim = numel(period);
         end
         
         % For rectangular grids, the periods p(k) form a sequence such that
         %     p(1) = 1,  p(k) divides p(k+1),  p(end) = ns.
         % The ratios p(k+1)/p(k) are the grid sizes.
         if ndim>1 && all(rem(nextperiod,period)==0)
            gridSize = nextperiod./period;
            gridVec = cell(1,nv);
            for ct=1:nv
               % Variable f{ct}
               k = ju(ct);       % associated grid dimension
               v = vals(:,ct);   % all data
               gvec = v(1:period(k):nextperiod(k));  % grid vector
               gap = reshape(v,gridSize) - reshape(gvec,[ones(1,k-1) numel(gvec) ones(1,ndim-k)]);
               if any(gap(:))
                  % Data is not consistent with rectangular grid
                  ndim = 1;   break
               else
                  gridVec{ct} = gvec;
               end
            end
         else
            ndim = 1;  % treat as 1D
         end
         
         % Build outputs
         if ndim>1
            % Multi-dimensional grid
            GV = cell(1,ndim);
            for ct=1:nv
               k = ju(ct);  % grid dimension
               GV{k} = cat(1,GV{k},{f{ct},gridVec{ct}});
               c{ct} = reshape(vals(:,ct),gridSize);
            end
         else
            % Scattered data or single-dimensional grid
            gridSize = [ns 1];
            for ct=1:nv
               c{ct} = vals(:,ct);
            end
            GV = {[f c]};
         end
         GridInfo = struct('GridSize',gridSize,'GridVectors',{GV},...
            'SamplePerm',SamplePerm,'Duplicate',Duplicate);
         RG = cell2struct(c,f,1);
      end
      
      function varargout = getGridVectors(SG)
         % Given a multi-dimensional sampling grid SG (struct), computes the 
         % grid vectors associated with each dimension. Variables whose value
         % is not constant along all but one grid dimensions are ignored,
         % and there can be more than one variable associated with a given
         % grid dimension. For example, for the 2D grid created by
         %    [alpha,beta] = ndgrid(alphav,betav)
         %    SG = struct('alpha',alpha,'beta',beta)
         % the call
         %   [GX,GY] = ltipack.SamplingGrid.getGridVectors(SG)
         % returns two structures
         %   GX = struct('alpha',alphav)
         %   GY = struct('beta',betav)
         % This can be seen as the inverse of the NDGRID operation.
         no = nargout;
         f = fieldnames(SG);
         nf = numel(f);
         if nf==0
            varargout = repmat({struct},[1 no]);
         else
            s = size(SG.(f{1}));
            ndim = numel(s);
            VarMap = zeros(nf,max(ndim,no));
            for dim=1:ndim
               VarMap(:,dim) = structfun(@(x) localVarIndicator(x,dim),SG);
            end
            % Consider only variables that vary only along one dimension
            VarMap(sum(VarMap,2)~=1,:) = 0;
            % Populate output arguments
            varargout = repmat({struct},[1 no]);
            ix0 = repmat({1},[1 ndim]);
            for dim=1:ndim
               ixv = find(VarMap(:,dim));
               nv = numel(ixv);
               if nv>0
                  vals = cell(nv,1);
                  ix = ix0;   ix{dim} = ':';
                  for ct=1:nv
                     vals{ct} = reshape(SG.(f{ixv(ct)})(ix{:}),[s(dim) 1]);
                  end
                  varargout{dim} = cell2struct(vals,f(ixv),1);
               end
            end
         end
      end
      
      function G = loadobj(s)
         G = s;
         G.Version_ = ltipack.ver();
      end
      
   end
   
end

%---------- local function

function boo = localVarIndicator(A,dim)
D = diff(A,1,dim);
boo = any(D(:)~=0);
end