approximize hE-vectors to the domain covered by an improper decomposition hcMatrix = projectHcVector(dec, heMatrix) see also: obtainDec, projectDec, mapToValidityCone
0001 function [heMatrixOut targetUnifiedPartition] = projectVectorDec(heMatrix, dec, zerotol, ... 0002 useTrueMapping) 0003 % approximize hE-vectors to the domain covered by an improper decomposition 0004 % 0005 % hcMatrix = projectHcVector(dec, heMatrix) 0006 % 0007 % see also: obtainDec, projectDec, mapToValidityCone 0008 0009 % The elk-library: convex geometry applied to crystallization modeling. 0010 % Copyright (C) 2014 Alexander Reinhold 0011 % 0012 % This program is free software: you can redistribute it and/or modify it 0013 % under the terms of the GNU General Public License as published by the 0014 % Free Software Foundation, either version 3 of the License, or (at your 0015 % option) any later version. 0016 % 0017 % This program is distributed in the hope that it will be useful, but 0018 % WITHOUT ANY WARRANTY; without even the implied warranty of 0019 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0020 % General Public License for more details. 0021 % 0022 % You should have received a copy of the GNU General Public License along 0023 % with this program. If not, see <http://www.gnu.org/licenses/> 0024 0025 %% input handling 0026 if ~exist('zerotol', 'var') 0027 zerotol = elkZerotol; 0028 end 0029 if ~exist('useTrueMapping', 'var') 0030 useTrueMapping = 0; 0031 end 0032 0033 if useTrueMapping 0034 error('elk:development', 'not yet supported'); 0035 %% true mapping 0036 heMatrixOut = decProj.properData.projectionMatrix * heMatrix; 0037 %! ToDo: The points in the heMatrix still need to be projected to the 0038 % validity cone of the improper decomposition 0039 %/ 0040 elseif ~dec.isComplete 0041 %% incomplete decomposition data (only true possible) 0042 heMatrixOut = dec.properData.projectionMatrix * heMatrix; 0043 else 0044 %% wrapped 0045 % ensure input is valid 0046 heMatrix = mapToValidityCone(heMatrix, dec.properData, 0); 0047 0048 % preset output matrix 0049 heMatrixOut = nan*heMatrix; 0050 0051 % find closest point in hE-space that matches a unified partition 0052 minDistanceVector = inf(1, size(heMatrix, 2)); 0053 thisDistanceVector = minDistanceVector; 0054 targetUnifiedPartition = nan(1, size(heMatrix, 2)); 0055 for iPart = 1:length(dec.unifiedPartition) 0056 0057 % orthogonal projection to selected subdomain 0058 orthogonalProjectionOnDomain = dec.unifiedPartition{iPart... 0059 }.projection.orthogonalOntoSubspace; 0060 heMatrixProjected = orthogonalProjectionOnDomain * heMatrix; 0061 % the resulting point must reside in the target domain to be 0062 % accepted 0063 inDomainA = dec.unifiedPartition{iPart}.embeddingDomain.hrep.A; 0064 filterInDomain = all(inDomainA * heMatrixProjected < zerotol, 1); 0065 % distance to domain 0066 thisDistanceVector(filterInDomain) = sum((... 0067 heMatrixProjected(:,filterInDomain) - ... 0068 heMatrix(:, filterInDomain)).^2, 1); 0069 thisDistanceVector(~filterInDomain) = Inf; 0070 0071 % evaluate current distances 0072 thisFilter = thisDistanceVector < minDistanceVector; 0073 minDistanceVector(thisFilter) = thisDistanceVector(thisFilter); 0074 targetUnifiedPartition(thisFilter) = iPart; 0075 heMatrixOut(:, thisFilter) = heMatrixProjected(:,thisFilter); 0076 end 0077 0078 %% second chance computations 0079 % From the above, we have projected the heVectors directly onto the 0080 % subspace in reach. We did considered the domain in reach in so far 0081 % that the projected heVector must reside in that domain. 0082 % The following loop handles cases in which the projection onto the 0083 % subspace in reach is not in the domain so that an additional 0084 % projection onto the boundary of the domain is required. As an 0085 % example think of an two-dimensional domain in 3D which has such a 0086 % boundary. 0087 for iPart = 1:length(dec.unifiedPartition) 0088 % orthogonal projection to selected subdomain 0089 orthogonalProjectionOnDomain = dec.unifiedPartition{iPart... 0090 }.projection.orthogonalOntoSubspace; 0091 heMatrixProjected = orthogonalProjectionOnDomain * heMatrix; 0092 % distance to domain 0093 thisDistanceVector = sum((heMatrixProjected - heMatrix).^2, 1); 0094 0095 % try to find closest vector in domain when there is a chance to be 0096 % better than the current values 0097 filterSecondTry = thisDistanceVector < minDistanceVector; 0098 if ~any(filterSecondTry) 0099 continue 0100 end 0101 % find suitable x0 0102 seIndexVector = dec.unifiedPartition{iPart}.seIndexVector; 0103 x0 = sum(dec.se(seIndexVector, :), 1)'; 0104 x0 = dec.unifiedPartition{iPart}.projection.validity*... 0105 dec.properData.mappingNewToProper*x0; 0106 inDomainA = dec.unifiedPartition{iPart}.embeddingDomain.hrep.A; 0107 % perform the projection 0108 heMatrixProjected(:, filterSecondTry) = ... 0109 closestProjectionOntoConfinementPartition(... 0110 heMatrix(:, filterSecondTry), ... 0111 x0, ... 0112 inDomainA, ... 0113 orthogonalProjectionOnDomain); 0114 if any(any(inDomainA * heMatrixProjected(:, filterSecondTry) > zerotol, 1)) 0115 error('elk:decomposition:numericalProblem', ['The he-Vectors ' ... 0116 'after this step must reside in the given domain. This ' .... 0117 'should not happen.']); 0118 end 0119 0120 % recompute distance vector.. ..this time only for the ones that 0121 % ARE in the domain.. ..where only vectors can be excluded that 0122 % already have a closest point assigned. 0123 thisDistanceVector(filterSecondTry) = sum((... 0124 heMatrixProjected(:,filterSecondTry) - ... 0125 heMatrix(:, filterSecondTry)).^2, 1); 0126 thisDistanceVector(~filterSecondTry) = Inf; 0127 0128 % evaluate current distances 0129 thisFilter = thisDistanceVector < minDistanceVector; 0130 minDistanceVector(thisFilter) = thisDistanceVector(thisFilter); 0131 targetUnifiedPartition(thisFilter) = iPart; 0132 heMatrixOut(:, thisFilter) = heMatrixProjected(:,thisFilter); 0133 end 0134 0135 % the above steps found the closest points: 0136 % * by orthogonal projection onto a unified partition presented in 0137 % the proper embedding (target domain) 0138 % * the points actually REALLY reside in that target domain 0139 % BUT: it can happen that no such "closest" point exists at all because 0140 % they must be projected ONTO some lowe dimensional subspace at once. 0141 % (finding the closest point on a polytope surface) 0142 0143 if any(isnan(targetUnifiedPartition)) 0144 % If this fails, this could mean that the decomposition is not 0145 % complete. 0146 error('elk:decomposition:numericalError', ['Could not find an ' ... 0147 'appropriate target partition for all input hE-vectors. ' ... 0148 'This should not happen.']); 0149 end 0150 0151 % handle all found target unified partitions subsequently 0152 % targetUnifiedPartition 0153 % hcMatrix = nan(dec.nC, size(heMatrix, 2)); 0154 % for iPart = unique(targetUnifiedPartition) 0155 % thisFilter = targetUnifiedPartition == iPart; 0156 % if ~any(thisFilter),continue,end 0157 % % projection onto an nC-dimensional subspace of hE-vectors that are 0158 % % described with the target hC-vectors 0159 % heMatrix(:, thisFilter) = ... 0160 % dec.unifiedPartition{iPart}.embedding.orthogonalProjection * ... 0161 % heMatrix(:, thisFilter); 0162 % 0163 % % ensure these vectors are valid: if it would not be valid after the 0164 % % projection, it might not be confined at the end 0165 % % heMatrix(:, thisFilter) = mapToValidityCone(heMatrix(:, thisFilter), ... 0166 % % dec.properData, zerotol); 0167 % 0168 % % projection onto the domain of possibly invalid hE-vectors that 0169 % % would be covered by mapNewToProper*hC 0170 % hcMatrix(:, thisFilter) = pinv(... 0171 % dec.unifiedPartition{iPart}.embedding.validityProjection*... 0172 % dec.properData.mappingNewToProper)*heMatrix(:, thisFilter); 0173 % end 0174 % 0175 % final projection into hC-space 0176 % hcMatrix = dec.properData.mappingProperToNew*heMatrix; 0177 end 0178 0179 if nargout < 2 0180 clear targetUnifiedPartition 0181 end 0182