compute projection matrices and validity regions for orthogonal mapping
Syntax: projectionData = generateOrthogonalProjectionDataForCone(A)
This function generates a set of projection matrices P and validity
regions {B*x <= 0} so that y=P*x is the closest point with (A*y <= 0)
for x. The output herefore is a data structure with fields:
.constraintMatrix - all possible rows of B
.projectionList - cell array with all projection matrices P
.responsibleMatrix - a matrix that assigns the rows of the constraint
matrix to the matrices B for each P. Each rows contains in the
first column a row-index of the constraintMatrix, in the second
column an index to a projection and in the third column +1 or -1
that might change the sign of the row in B (see example).
.targetDimensionVector - the target dimension of each projection
Alternatively, you can use the inputs ..(A, extremeRayMatrix, zerotol) to
provide the a prior calculated extreme rays of the cone A and set the
tolerance.
See also: obtainDec, approximateDec, projectFromData0001 function projectionData = generateOrthogonalProjectionDataForCone(A, ... 0002 extremeRayMatrix, zerotol) 0003 % compute projection matrices and validity regions for orthogonal mapping 0004 % 0005 % Syntax: projectionData = generateOrthogonalProjectionDataForCone(A) 0006 % 0007 % This function generates a set of projection matrices P and validity 0008 % regions {B*x <= 0} so that y=P*x is the closest point with (A*y <= 0) 0009 % for x. The output herefore is a data structure with fields: 0010 % .constraintMatrix - all possible rows of B 0011 % .projectionList - cell array with all projection matrices P 0012 % .responsibleMatrix - a matrix that assigns the rows of the constraint 0013 % matrix to the matrices B for each P. Each rows contains in the 0014 % first column a row-index of the constraintMatrix, in the second 0015 % column an index to a projection and in the third column +1 or -1 0016 % that might change the sign of the row in B (see example). 0017 % .targetDimensionVector - the target dimension of each projection 0018 % 0019 % Alternatively, you can use the inputs ..(A, extremeRayMatrix, zerotol) to 0020 % provide the a prior calculated extreme rays of the cone A and set the 0021 % tolerance. 0022 % 0023 % See also: obtainDec, approximateDec, projectFromData 0024 0025 % The elk-library: convex geometry applied to crystallization modeling. 0026 % Copyright (C) 2013 Alexander Reinhold 0027 % 0028 % This program is free software: you can redistribute it and/or modify it 0029 % under the terms of the GNU General Public License as published by the 0030 % Free Software Foundation, either version 3 of the License, or (at your 0031 % option) any later version. 0032 % 0033 % This program is distributed in the hope that it will be useful, but 0034 % WITHOUT ANY WARRANTY; without even the implied warranty of 0035 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0036 % General Public License for more details. 0037 % 0038 % You should have received a copy of the GNU General Public License along 0039 % with this program. If not, see <http://www.gnu.org/licenses/>. 0040 0041 %% input handling 0042 0043 % zerotol 0044 if ~exist('zerotol', 'var') || isempty(zerotol) 0045 zerotol = elkZerotol; 0046 end 0047 0048 nDim = size(A,2); 0049 0050 % extreme ray matrix 0051 if ~exist('extremeRayMatrix', 'var') || isempty(extremeRayMatrix) 0052 hrep.A = [A; mean(A, 1)]; 0053 hrep.h = [zeros(size(A, 1), 1); 1]; 0054 vrep = convertHrepToVrep(hrep); 0055 extremeRayMatrix = vrep.V; 0056 % eliminate origin 0057 extremeRayMatrix = eliminateRows(extremeRayMatrix, ... 0058 zeros(1, size(A, 2)), zerotol); 0059 end 0060 0061 %% Catch empty cones 0062 % The provided cone might be trivial so that only the origin itself 0063 % fulfills A*x <= 0. We notify this, if we try to calculate the 0064 % V-representation, which results in a single vertex (the origin) in that 0065 % case. 0066 vrep = convertHrepToVrep(struct('A', A, 'h', zeros(size(A, 1), 1))); 0067 if size(vrep.V, 1) == 1 && all(vrep.V < zerotol) 0068 projectionData.A = [1 zeros(nDim-1, 0)]; 0069 projectionData.projectionList = {zeros(nDim, nDim), zeros(nDim, nDim)}; 0070 projectionData.responsibleMatrix = [1 1 1; 1 2 -1]; 0071 projectionData.targetDimensionVector = [0 0]; 0072 return 0073 end 0074 0075 %% Trivial cases 0076 % especially simple are the cases where P is either the zero-matrix and 0077 % where P is the idendity matrix. 0078 0079 % In case of P being the identity matrix, the point x is already valid. 0080 projectionData.A = A; 0081 projectionData.responsibleMatrix = [(1:size(A,1))', ones(size(A, 1), 1), ... 0082 ones(size(A, 1), 1)]; 0083 projectionData.projectionList = {eye(nDim)}; 0084 projectionData.targetDimensionVector = nDim; 0085 0086 % In case of P being the zero-matrix, the closest point is the origin. 0087 % This projection is valid for any point in the positive cone of the 0088 % normal vectors in A. 0089 vrep.V = [A; zeros(1, nDim)]; 0090 hrep = convertVrepToHrep(vrep); 0091 % filter facets that do not hit 0 0092 hrep.A(abs(hrep.h) > zerotol, :) = []; 0093 % assign output (we accept double entries, they are unlikely) 0094 projectionData = appendProjectionData(projectionData, hrep.A, ... 0095 zeros(nDim, nDim), zerotol); 0096 projectionData.targetDimensionVector(end+1) = 0; 0097 0098 %% adjacency information 0099 % which facet (row) is adjacent to which extreme ray (column) 0100 adjacencyMatrix = abs(A * extremeRayMatrix') < zerotol; 0101 0102 % obtain all combinations of 0 up to nC selected elements from the 0103 % boundaries in A 0104 permutationMatrix = concatenateOrderetSetMatrices(size(A,1), nDim, 1); 0105 %! DEV: I've put concatenateOrderetSetMatrices as a seperate function under 0106 % private and adopted it to also return selecting nDim elements. This 0107 % should change things for here, even though it matches the expectation 0108 % of what concatenateOrderetSetMatrices should do. 0109 0110 % initialize 0111 % constraintMatrix = zeros(0, nDim); 0112 % responsibleMatrix = zeros(0, 3); 0113 % projectionList = {}; 0114 % targetDimensionVector = []; 0115 0116 for iPerm = 1:size(permutationMatrix, 1) 0117 %% general case information 0118 % we use this permutation, now 0119 thisPerm = permutationMatrix(iPerm, :); 0120 % delete 'no facet selected' entries 0121 thisPerm(thisPerm == 0) = []; 0122 0123 % basic information 0124 thisNullSpaceDim = length(thisPerm); 0125 thisTargetDim = nDim - length(thisPerm); 0126 thisNormalMatrix = A(thisPerm, :)'; 0127 0128 % skip invalid combinations: each set of facet normals must have 0129 % targetDim extreme rays in common. The intersection of nF facets 0130 % leaves an (nDim - nF) dimensional subspace that is defined by (nDim 0131 % - nF) points, the extreme points. 0132 filterMatrix = adjacencyMatrix(thisPerm, :); 0133 filterVector = sum(filterMatrix, 1) == thisNullSpaceDim; 0134 if ~isempty(thisPerm) && (sum(filterVector) < thisTargetDim) 0135 % the nF facets are not adjacent to (nDim-nF) extreme rays AND they 0136 % are not empty: skip this case 0137 continue 0138 end 0139 0140 %% projection and condition 0141 % projection 0142 thisNullProjection = thisNormalMatrix * ... 0143 inv(thisNormalMatrix'*thisNormalMatrix)*thisNormalMatrix'; 0144 thisProjection = eye(nDim) - thisNullProjection; 0145 0146 % condition 0147 thisConstraintMatrix = [A*thisProjection; ... 0148 -1*pinv(thisNormalMatrix)]; 0149 % filter zero rows (happens for the facets that is projected onto) 0150 thisConstraintMatrix = eliminateRows(thisConstraintMatrix, ... 0151 zeros(1, nDim), zerotol); 0152 % normalize again 0153 thisConstraintMatrix = bsxfun(@times, thisConstraintMatrix, ... 0154 1./sqrt(sum(thisConstraintMatrix.^2, 2))); 0155 0156 %% add information 0157 % add projection 0158 projectionData = appendProjectionData(projectionData, ... 0159 thisConstraintMatrix, thisProjection, zerotol); 0160 projectionData.targetDimensionVector(end+1) = thisTargetDim; 0161 0162 end