generateConfinementPartition

 generate a unified partition with confinement condition

 THIS IS NO USER FUNCTION

• computePartitionIndexDec determine the partition of a given hC-vector
• obtainProjectionDataForVector extract projection matrices and validity regions
• trimConstraintMatrix normalize and trim rows and (optionally) remove zero-rows

• projectCompleteDec project a complete proper representation

• function [projectionMatrix constraintMatrix] =

0001 function partitionData = generateConfinementPartition(dec, hC, zerotol)
0002 % generate a unified partition with confinement condition
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 %% initialize acuumulating array structure
0023 if exist('hC', 'var')
0024     partitionData.sampleHc = hC;
0025 else
0026     partitionData.sampleHc = [];
0027 end
0028 partitionData.flagResult = [];
0029 partitionData.hrep = [];
0030 partitionData.isValid = [];
0031 partitionData.isConfined = [];
0032 partitionData.unifiedIndex = [];
0033 partitionData.validityProjectionMatrix = [];
0034 partitionData.validityProjectionConstraintMatrix = [];
0035 partitionData.violatedConstraintIndexVector = [];
0036 % partitionData.confinementProjectionMatrix = [];
0037 % partitionData.confinementConstraintMatrix = [];
0038 partitionData.debug = [];
0039     
0040 % quit if no input is given; decomposeCone uses this call to initialize
0041 %   the accumulating array
0042 if ~exist('dec', 'var')
0043     partitionData(:) = [];
0044     return
0045 end
0046 
0047 %% precalculation
0048 % vector in proper embedding space
0049 hE = dec.properData.mappingNewToProper * hC;
0050 
0051 % projection data
0052 projectionData = obtainProjectionDataForVector(dec.properData, hE, 0);
0053 % save the validity projection for later measure mapping
0054 partitionData.validityProjectionMatrix = projectionData.projectionMatrix;
0055 partitionData.validityProjectionConstraintMatrix = ...
0056     projectionData.constraintMatrix;
0057 % also save the violated constraints
0058 partitionData.violatedConstraintIndexVector = ...
0059     projectionData.violatedConstraintIndexVector;
0060 
0061 % what is the unified partition for hE' ?
0062 unifiedIndex = computePartitionIndexDec(dec.properData, ...
0063     projectionData.projectionMatrix*hE, zerotol, 1);
0064 
0065 %% Cases
0066 % Everything that follows is still in hE-space. The transformation to
0067 %   hC-space is equal for all cases.
0068 if all((dec.properData.validityCone.A * hE) < 0)
0069     %% CASE:  hC is valid
0070 %     disp('DBG: valid');
0071     % this will be a valid parittion
0072     partitionData.flagResult = 1;
0073     % but let's remember that its an originally valid one
0074     partitionData.isValid = 1;
0075     partitionData.isConfined = 1;
0076     
0077     % the hrep is the unified partition itself; projected to hC later
0078     thisA = dec.properData.unifiedPartition{unifiedIndex}.A;
0079     
0080     % we save the unified index
0081     partitionData.unifiedIndex = unifiedIndex;
0082      
0083 elseif rank(projectionData.projectionMatrix * dec.properData.mappingNewToProper, ...
0084         zerotol) == dec.nC
0085     %% CASE: confined hC vector but not valid
0086 %     disp('DBG: confined');
0087     % this will be a valid partition
0088     partitionData.flagResult = 1;
0089     % but not valid
0090     partitionData.isValid = 0;
0091     partitionData.isConfined = 1;
0092     
0093     % add the constraint such that the valid hE always resides in THIS
0094     %   unified partition.
0095     thisA = [projectionData.constraintMatrix; ...
0096              dec.properData.unifiedPartition{unifiedIndex}.A*...
0097              projectionData.projectionMatrix];  
0098          
0099     % we save the unified index
0100     partitionData.unifiedIndex = unifiedIndex;
0101     
0102 else
0103     %% unconfined hC vector
0104 %     disp('DBG: unconfined');
0105     % this partition will also be accepted, but is is neither valid nor
0106     %   confined. But we will store it together with the precalculated
0107     %   projection data.
0108     partitionData.flagResult = 1;
0109     partitionData.isValid = 0;
0110     partitionData.isConfined = 0;
0111     
0112     %! Dev: This was the code previously used. But since we now include
0113     %    unconfined partitions in the decomposition, we need to consider
0114     %    all constraints like for the confined partitions above.
0115     % removed part:
0116 %     % add the constraint that the projection is valid, but do not store the
0117 %     %   target unified partition as it is not relevant.
0118 %     thisA = projectionData.constraintMatrix;
0119     % new part
0120     thisA = [projectionData.constraintMatrix; ...
0121              dec.properData.unifiedPartition{unifiedIndex}.A*...
0122              projectionData.projectionMatrix];
0123     % we save the unified index
0124     partitionData.unifiedIndex = unifiedIndex;
0125     %/
0126     
0127 end
0128 
0129 %% Final conclusion
0130 % save original constraint matrix (in proper embedding)
0131 partitionData.debug.embeddingConstraintMatrix = thisA;
0132 
0133 % project constraints to new space
0134 thisA = thisA*dec.properData.mappingNewToProper;
0135 % while usually a lot of redundant facets will appear, especially
0136 %   considering the transformation of the nonempty-cone, these are no
0137 %   problem. Problematic are close-to zero projections where the facet
0138 %   of the cone is almost parallel to the projection.
0139 %
0140 % remove zero rows and trim
0141 thisA = trimConstraintMatrix(thisA, zerotol, 1);
0142 % remove redundant rows
0143 thisA = unique(thisA, 'rows');
0144 
0145 % accumulate now:
0146 partitionData.hrep.A = thisA;
0147 partitionData.hrep.h = zeros(size(thisA, 1), 1);
0148 partitionData.debug.projectionData = projectionData;
0149 
0150 
0151 end
0152 
0153 function [projectionMatrix constraintMatrix] = ...
0154     concatenateProjectionData(projectionData, filterVector)
0155 % concatenate projection matrices according to a filter vector
0156 %
0157 % Input are the cell arrays of the single projection matrices and
0158 %   constraints for this projection in hE-space.
0159 %
0160 % THIS IS NO USER FUNCTION
0161                           
0162 % input
0163 if ~exist('filterVector', 'var')
0164     filterVector = (1:length(projectionData.projectionList))>0;
0165 end
0166 properNc = size(projectionData.projectionList{1}, 1);
0167 
0168 % init data
0169 projectionMatrix = eye(properNc);
0170 constraintMatrix = zeros(0, properNc);
0171 
0172 % concatenate data
0173 for iProj = find(filterVector)
0174     constraintMatrix = [constraintMatrix; ...
0175                         projectionData.constraintList{iProj}];
0176     projectionMatrix = projectionMatrix * ...
0177                        projectionData.projectionList{iProj}; 
0178 end
0179 
0180 end