addMonomialMatrixDec

% generate monomial matrices

0001 function dec = addMonomialMatrixDec(dec, nMax)
0002 
0003 %% generate monomial matrices
0004 mat = cell(nMax, 1);
0005 % mat{1} = [];%zeros(1, dec.nC);
0006 nMonomial = 0;
0007 
0008 for iExpo = 1:nMax
0009     % obtain indices for exponent
0010     mat{iExpo} = obtainOrderedSetMatrix(dec.nC, iExpo, 1);
0011     nRow = size(mat{iExpo}, 1);
0012     nMonomial = nMonomial + nRow;
0013     
0014 %     % exponent index matrix
0015 %     mat{iExpo+1} = zeros(nRow, dec.nC);
0016 %     for iRow = 1:nRow
0017 %         mat{iExpo+1}(iRow
0018 %
0019 %     mat = [mat;
0020 end
0021 
0022 nMonomial = nMonomial + 1;
0023 dec.monomialMatrix = zeros(nMonomial, nMax);
0024 iRow = 0;
0025 for iExpo = nMax:-1:1
0026     dec.monomialMatrix(iRow + (1:size(mat{iExpo}, 1)), 1:iExpo) = mat{iExpo};
0027     iRow = iRow + size(mat{iExpo}, 1);
0028 end
0029 
0030 
0031 %% construct derivative matrix
0032 
0033 % will be: monomialDerivativeMatrix
0034 matIndex = zeros(nMonomial, nMonomial);
0035 matFactor = zeros(nMonomial, nMonomial);
0036 
0037 for iExpo = 1:nMax
0038     % we get the partial derivative monomial matrix
0039     derivativeMonomialMatrix = dec.monomialMatrix(:, ...
0040         [(1:(iExpo-1)) ((iExpo+1):end)]);
0041     % we add a zero column to obtain the same configuration as in
0042     % monomialMatrix
0043     derivativeMonomialMatrix = [derivativeMonomialMatrix zeros(nMonomial, 1)];
0044     
0045     % now we need to find the index to monomialMatrix that matches each
0046     % line in derivativeMonomialMatrix
0047     colVector = zeros(nMonomial, 1);
0048     for iMono = 1:nMonomial
0049         filterMatrix = bsxfun(@eq, dec.monomialMatrix, ...
0050             derivativeMonomialMatrix(iMono, :));
0051         filterVector = sum(filterMatrix, 2) >= size(dec.monomialMatrix, 2);
0052         if sum(filterVector) > 1 || ~any(filterVector)
0053             error('elk:decMonomial:internalError', ['The construction of the ' ...
0054                 'monomial derivatives failed. The partial derivative should ' ...
0055                 'be unique in the monomialMatrix but is not. This should not ' ...
0056                 'happen.']);
0057         end
0058         colVector(iMono) = find(filterVector);
0059     end
0060     % row vector is obvious
0061     rowVector = (1:nMonomial)';
0062     
0063     % but a partial derivative where monomialMatrix has a zero element
0064     %  should not be considered
0065     thisFilter = dec.monomialMatrix(:, iExpo) == 0;
0066     colVector(thisFilter) = [];
0067     rowVector(thisFilter) = [];
0068     
0069     % convert to indices into derivative matrix data
0070     matIndexVector = sub2ind(size(matIndex), rowVector, colVector);
0071     
0072     % add required indices - note that indices might be written multiple
0073     %   times but only identical indices can be written for each row
0074     matIndex(matIndexVector) = dec.monomialMatrix(~thisFilter, iExpo);
0075     
0076     % accumulate factor
0077     matFactor(matIndexVector) = matFactor(matIndexVector) + 1;
0078 end
0079 
0080 dec.monomialDerivativeIndexMatrix = matIndex;
0081 dec.monomialDerivativeFactorMatrix = matFactor;