project points based on precalculated projection data
Syntax: xMatrix = projectOrthogonalFromData(xMatrix, projectionData)
The input/output xMatrix contains vectors in each column. The
projectionData contains several possible projections as a cell array
.projectionList. It also contains a set of linear constraints in the
rows of projectionData.A for which .responsibleMatrix indicates which
constraint line (1st column) belongs to which projection (2nd column)
with a possible sign change (3rd column). Despite of the selected
subset, the indicated projection is chosen by:
projectionData.A * xMatrix <= 0
with sign=1 and by:
projectionData.A * xMatrix > 0
with sign=-1.
The additional argument is checkCompletion enable (=1) warnings when not
exactly one projection is selected. For checkCompletion=2, the warnings
are followed by an error.
THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE
OF APPLICATIONS (decomposition case study; technical reports)
See also: generateOrthogonalProjectionDataForCone0001 function [xMatrixOutT projectionIndexVector constraintMatrixList] = ... 0002 projectFromData(xMatrixT, projectionData, zerotol, checkCompletion, ... 0003 xMatrixOutT) 0004 % project points based on precalculated projection data 0005 % 0006 % Syntax: xMatrix = projectOrthogonalFromData(xMatrix, projectionData) 0007 % 0008 % The input/output xMatrix contains vectors in each column. The 0009 % projectionData contains several possible projections as a cell array 0010 % .projectionList. It also contains a set of linear constraints in the 0011 % rows of projectionData.A for which .responsibleMatrix indicates which 0012 % constraint line (1st column) belongs to which projection (2nd column) 0013 % with a possible sign change (3rd column). Despite of the selected 0014 % subset, the indicated projection is chosen by: 0015 % projectionData.A * xMatrix <= 0 0016 % with sign=1 and by: 0017 % projectionData.A * xMatrix > 0 0018 % with sign=-1. 0019 % 0020 % The additional argument is checkCompletion enable (=1) warnings when not 0021 % exactly one projection is selected. For checkCompletion=2, the warnings 0022 % are followed by an error. 0023 % 0024 % THIS FUNCTION IS PUBLISHED BUT TESTED AND MAINTAINED ONLY IN SCOPE 0025 % OF APPLICATIONS (decomposition case study; technical reports) 0026 % 0027 % See also: generateOrthogonalProjectionDataForCone 0028 0029 % The elk-library: convex geometry applied to crystallization modeling. 0030 % Copyright (C) 2013 Alexander Reinhold 0031 % 0032 % This program is free software: you can redistribute it and/or modify it 0033 % under the terms of the GNU General Public License as published by the 0034 % Free Software Foundation, either version 3 of the License, or (at your 0035 % option) any later version. 0036 % 0037 % This program is distributed in the hope that it will be useful, but 0038 % WITHOUT ANY WARRANTY; without even the implied warranty of 0039 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0040 % General Public License for more details. 0041 % 0042 % You should have received a copy of the GNU General Public License along 0043 % with this program. If not, see <http://www.gnu.org/licenses/>. 0044 0045 0046 %% Input handling 0047 if ~exist('xMatrixOutT', 'var') 0048 xMatrixOutT = xMatrixT; 0049 end 0050 if size(xMatrixT, 1) ~= size(projectionData.projectionList{1}, 1) || ... 0051 size(xMatrixOutT, 1) ~= size(projectionData.projectionList{1}, 1) 0052 error('elk:decomposition:wrongInput', ['Please, check the provided ' ... 0053 'xMatrix, the vectors are the columns of this matrix and must ' ... 0054 'equal the dimensions of the projection matrices']); 0055 end 0056 if ~exist('zerotol', 'var') || isempty(zerotol) 0057 zerotol = 0; 0058 end 0059 if ~exist('checkCompletion', 'var') || isempty(checkCompletion) 0060 checkCompletion = 0; 0061 end 0062 0063 %% Initialize 0064 % return vector for used projections 0065 if nargout >= 2 0066 projectionIndexVector = zeros(1, size(xMatrixT, 2)); 0067 end 0068 if nargout >= 3 0069 constraintMatrixList = cell(1, size(xMatrixT, 2)); 0070 end 0071 % vector for checking whether it was projected properly 0072 numberProjectionsVector = zeros(1, size(xMatrixT, 2)); 0073 0074 %% Evaluate constraints 0075 % determine which constraints are fulfilled, and which not - this CAN be 0076 % based on a second matrix, when provided. This behavior is required for 0077 % the grwoth rate limitations. 0078 constraintMatrix = projectionData.A*xMatrixT; 0079 0080 % constraints are evaluated in favour of the negative halfspace; zerotol 0081 % must be applied because the demand of a projection already uses this 0082 % zerotol. 0083 %! DEV: We change this decision such that a double-projection is possible - 0084 % but only in regions where the projection should be almost the same, 0085 % anyway.. ..in the scope of zerotol. 0086 % old code: 0087 % constraintMatrixNegative = constraintMatrix <= zerotol; 0088 % constraintMatrixPositive = constraintMatrix > zerotol; 0089 % new code: 0090 constraintMatrixNegative = bsxfun(@le, constraintMatrix, zerotol); 0091 constraintMatrixPositive = bsxfun(@ge, constraintMatrix, -1*zerotol); 0092 %/ 0093 0094 %% Find projection 0095 % To execute the projection, it is necessary to cycle either through the 0096 % vectors or through the projections. The number of projections should 0097 % typically be less in applications. 0098 for iProjection = 1:length(projectionData.projectionList) 0099 thisConstraintIndexVector = projectionData.responsibleMatrix(... 0100 projectionData.responsibleMatrix(:, 2) == iProjection, 1); 0101 thisConstraintSignVector = projectionData.responsibleMatrix(... 0102 projectionData.responsibleMatrix(:, 2) == iProjection, 3); 0103 % the projection is valid, when all constraints are fulfilled 0104 thisFilterMatrix = ... 0105 (bsxfun(@or, thisConstraintSignVector < 0, ... 0106 constraintMatrixNegative(thisConstraintIndexVector, :)) & ... 0107 bsxfun(@or, thisConstraintSignVector > 0, ... 0108 constraintMatrixPositive(thisConstraintIndexVector, :))); 0109 thisFilterVector = all(thisFilterMatrix, 1); 0110 0111 if ~isempty(thisFilterVector) 0112 xMatrixOutT(:, thisFilterVector) = ... 0113 projectionData.projectionList{iProjection} * ... 0114 xMatrixOutT(:, thisFilterVector); 0115 if nargout >= 2 0116 projectionIndexVector(thisFilterVector) = iProjection; 0117 end 0118 if nargout >= 3 0119 [constraintMatrixList{thisFilterVector}] = deal(... 0120 bsxfun(@times, ... 0121 projectionData.A(thisConstraintIndexVector, :), ... 0122 thisConstraintSignVector(:))); 0123 end 0124 numberProjectionsVector(thisFilterVector) = ... 0125 numberProjectionsVector(thisFilterVector) + 1; 0126 end 0127 end 0128 0129 %% error handling 0130 flagError = 0; 0131 if checkCompletion 0132 %! Info: any of the following warnings / errors is only thrown with 0133 % respect to validity mapping. Otherwise, the "should not happen" is 0134 % wrong. 0135 % if max(numberProjectionsVector) > 1 0136 % warning('elk:decomposition:general', ['A vector was ' ... 0137 % 'projected more than once. This should not happen.']); 0138 % flagError = 1; 0139 % end 0140 if min(numberProjectionsVector) < 1 0141 warning('elk:decomposition:general', ['A vector was ' ... 0142 'not projected. This should not happen.']); 0143 flagError = 1; 0144 end 0145 if checkCompletion == 2 && flagError 0146 error('elk:decomposition:numericalError', 'See the warning above.'); 0147 end 0148 end