collect information required for computeProjectionFromGrass THIS IS NO USER FUNCTION
0001 function grassStruct = obtainGrassStruct(hcMatrix, par) 0002 % collect information required for computeProjectionFromGrass 0003 % 0004 % THIS IS NO USER FUNCTION 0005 0006 % The elk-library: convex geometry applied to crystallization modeling. 0007 % Copyright (C) 2013 Alexander Reinhold 0008 % 0009 % This program is free software: you can redistribute it and/or modify it 0010 % under the terms of the GNU General Public License as published by the 0011 % Free Software Foundation, either version 3 of the License, or (at your 0012 % option) any later version. 0013 % 0014 % This program is distributed in the hope that it will be useful, but 0015 % WITHOUT ANY WARRANTY; without even the implied warranty of 0016 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0017 % General Public License for more details. 0018 % 0019 % You should have received a copy of the GNU General Public License along 0020 % with this program. If not, see <http://www.gnu.org/licenses/>. 0021 0022 %% Coordinate transformation 0023 % We use PCA to align the main variances to the coordinate axes. The 0024 % diagonal elements of S represent the variance along the 0025 % PCA-coordinates. While U contains the PCA coordinates as column 0026 % vectors. 0027 if numel(par.manualScaling) == 1 0028 % singular value decomposition 0029 [U,S] = svd(hcMatrix, 0); 0030 % We want to represent the original points in positionMatrix such that 0031 % they fill a well defined large area. Therefore, we map them to the 0032 % scaled PCA-space, with: 0033 % y = S^(-1) * U' * x 0034 % The variance in each PCA component is now 1. 0035 % But beforehand S^(-1) is problematic for close-to-zero eigenvalues. 0036 % So that we fix that manually: 0037 0038 %! DEV: try to relax the scaling: 0039 % S = diag(diag(S)); 0040 % S = 0.9*S + 0.1*eye(size(hcMatrix, 1)); 0041 S = S.^par.manualScaling; 0042 %/ 0043 0044 S(abs(S(:)) < par.zerotol) = 0; 0045 Sinf = diag(1./diag(S)); 0046 Sinf(isinf(Sinf(:))) = 0; 0047 % mapping to PCA space: 0048 grassStruct.mapToPcaSpace = Sinf * U'; 0049 % To project coordinates in the scaled PCA-space back to original 0050 % coordinates, we use: 0051 grassStruct.mapToRealSpace = U * diag(diag(S)); 0052 %! DEV: is mapToPcaSpace required at all?? 0053 %/ 0054 % The initial range is now simply represented by the first targetDim 0055 % coordinate axes in scaled space. 0056 grassStruct.baseRange = eye(size(hcMatrix, 1)); 0057 grassStruct.baseRange(:, (par.targetDim+1):end) = []; 0058 % and the null space are the remaining coordinates 0059 grassStruct.baseNull = eye(size(hcMatrix, 1)); 0060 grassStruct.baseNull(:, 1:par.targetDim) = []; 0061 0062 else 0063 % if the scaling matrix is already provided, we have: 0064 grassStruct.mapToPcaSpace = par.manualScaling; 0065 % the scaling might not have full rank (see above PCA scaling approach 0066 % when the date itself is already singular); hence, we schould use 0067 % pinv to be safe, here: 0068 grassStruct.mapToRealSpace = pinv(par.manualScaling, par.zerotol); 0069 0070 % the initial subspace must now be selected from PCA projection and 0071 % manually mapped to the scaled space 0072 pcaData = approxPca(hcMatrix, par); 0073 0074 % Let y be a point in scaled-PCA space and x the corresponding point in 0075 % the original real space. Let further denote xP and yP the corresponding 0076 % projected points. We also have the original projection matrix P: 0077 % xP = P * x 0078 % and an unknown projection matrix Ps in scaled-PCA space: 0079 % yP = Ps y 0080 % This unknown projection matrix can be determined, by: 0081 % xP = P *x 0082 % mapToRealSpace yP = P * mapToRealSpace * y 0083 % yP = mapToPcaSpace * P * mapToRealSpace * y 0084 % Ps = mapToPcaSpace * P * mapToRealSpace 0085 scaledProjectionMatrix = grassStruct.mapToPcaSpace * ... 0086 pcaData.projectionMatrix * grassStruct.mapToRealSpace; 0087 0088 % from these some base vectors for the range and null space 0089 %! DEV: this code was originally used, here but turned out to be 0090 % numerically instable. 0091 grassStruct.baseRange = orth(scaledProjectionMatrix); 0092 grassStruct.baseNull = null(scaledProjectionMatrix); 0093 % The following is adopted from orth/null, but enforced to use the right 0094 % dimensions 0095 % [U, S, V] = svd(initProjectionScaledPca, 0) 0096 % grassStruct.baseRange = U(:, 1:par.targetDim); 0097 % grassStruct.baseNull = V(:, (par.targetDim+1):end); 0098 end 0099 0100 0101 %! DEV: new approach using manual initial condition (deprecated): did only 0102 % shift the Grass space in directions but did not yield any 0103 % distortion. 0104 %/ DEV: the below is very unfortunate for the PCA approach since the base 0105 % and null matrices are actually party of the unity matrix. 0106 % grassStruct.baseRange = orth(grassStruct.mapToPcaSpace*par.manualInitialCondition); 0107 % grassStruct.baseNull = null(grassStruct.baseRange'); 0108 %/ 0109 0110 % grassmann dimensions 0111 grassStruct.nDim = size(grassStruct.baseRange, 1); 0112 grassStruct.targetDim = size(grassStruct.baseRange, 2); 0113 grassStruct.grassDim = (grassStruct.nDim - grassStruct.targetDim) * ... 0114 grassStruct.targetDim; 0115 grassStruct.zerotol = par.zerotol; 0116 0117 0118 %% Data bounding cone 0119 %! DEV: is this still required?? In high dimensions this part will take 0120 % extremely long. And so far, I think the nonlinear constraint wa 0121 % swtiched off, anyway. 0122 % >> need a switch for turning on/off this feature - at least the 0123 % documentation requires the nonlincon to be built. But it is absolutely 0124 % awquard 0125 % if isempty(par.manualBoundingCone) 0126 % 0127 % % these lines are seperated: 0128 % [grassStruct.extremePointMatrix, ... 0129 % grassStruct.coneNormalVectorMatrix] = ... 0130 % computeDataBoundingCone(grassStruct, dataOriginal.positionMatrix); 0131 % 0132 % else 0133 % if isstruct(par.manualBoundingCone) 0134 % % hrep provided 0135 % grassStruct.coneNormalVectorMatrix = par.manualBoundingCone.A; 0136 % vrep = convertHrepToVrep(par.manualBoundingCone); 0137 % grassStruct.extremePointMatrix = vrep.V'; 0138 % else 0139 % % generating rays provided 0140 % grassStruct.extremePointMatrix = par.manualBoundingCone; 0141 % hrep = convertVrepToHrep(struct('V', ... 0142 % [grassStruct.extremePointMatrix'; ... 0143 % zeros(grassStruct.nDim, 1)] ... 0144 % )); 0145 % grassStruct.coneNormalVectorMatrix = hrep.A; 0146 % end 0147 end 0148 %/ 0149 0150 %! Debug 0151 % figure(5) 0152 % viewPolytope(struct('A', grassStruct.coneNormalVectorMatrix', ... 0153 % 'h', zeros(size(grassStruct.coneNormalVectorMatrix, 2), 1)), ... 0154 % 'infinityBoxSize', 1) 0155 % plot3(grassStruct.extremePointMatrix(1,:), ... 0156 % grassStruct.extremePointMatrix(2,:), ... 0157 % grassStruct.extremePointMatrix(3,:), 'rx', 'markerSize', 15); 0158 % figure(12),hold on 0159 % viewPolytope(struct('A', [coneNormalVectorMatrix'; 1 0 0], ... 0160 % 'h', ... 0161 % [zeros(size(grassStruct.coneNormalVectorMatrix, 2), 1); 1]), ... 0162 % 'infinityBoxSize', 1) 0163 % % 0164 % set(gca, 'PlotBoxAspectRatioMode', 'auto', ... 0165 % 'DataAspectRatioMode', 'auto'); 0166 % xlim([-1.5 0.5]);ylim([-1 1]);zlim([-1 1]); 0167 %/