decompose the cone of valid polytopes into regions of equal a-type rayMatrix contains a "scaledToClose" extreme ray in each row. THIS FUNCTION IS NO USER FUNCTION.
0001 function partitionData = decomposeCone(dec, A, genericPartition, ... 0002 zerotol, verbose) 0003 % decompose the cone of valid polytopes into regions of equal a-type 0004 % 0005 % rayMatrix contains a "scaledToClose" extreme ray in each row. 0006 % 0007 % THIS FUNCTION IS NO USER FUNCTION. 0008 0009 % The elk-library: convex geometry applied to crystallization modeling. 0010 % Copyright (C) 2012 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 % the input cone is provided in hrep with its facet normals in A; we 0027 % require the closed cone 0028 hrepCone = struct('A', A, 'h', zeros(size(A, 1), 1)); 0029 hrepClosed = closeCone(hrepCone, dec.centerRay); 0030 coneVolume = computeVolumeHrep(hrepClosed); 0031 % this is also the starting remaining domain 0032 remainingDomain = {hrepClosed}; 0033 0034 %% Initializing 0035 % initialize random number generator 0036 randState = rng; 0037 rng(1); 0038 % loop initialization 0039 partitionData = genericPartition(); 0040 gatheredVolume = 0; 0041 0042 if verbose,fprintf('......');end 0043 0044 while 1 == 1 0045 if verbose 0046 fprintf('\b\b\b\b\b\b%5.1f%%', (gatheredVolume/coneVolume*100)); 0047 % fprintf('%5.1f%%, %f ', (gatheredVolume/coneVolume*100), ... 0048 % length(remainingDomain)); 0049 end 0050 0051 %% Decompose a sample 0052 % We randomly generate a hC-vector from the last hrep that is entered. 0053 % This hrep can be assumed to be relatively small according to the 0054 % polytope difference algorithm. For the random vector, an offset of 0055 % 0.1 is applied so that the resulting random vector is safely in the 0056 % interior of the selected cone. 0057 thisVrep = convertHrepToVrep(remainingDomain{end}); 0058 thisVrep = openCone(thisVrep, zerotol); 0059 lam = rand([1, size(thisVrep.V, 1)]) + 0.1; 0060 hC = (lam * thisVrep.V)'; 0061 0062 %% use generic partition generator 0063 thisPartitionData = genericPartition(dec, hC, zerotol); 0064 0065 % repeat, if failed 0066 if thisPartitionData.flagResult == 0, continue, end; 0067 0068 % check if sample is part of partition 0069 if any((thisPartitionData.hrep.A * hC) > zerotol) 0070 error('decomposition:elk:numericalError', ['Sample is not part of ' ... 0071 'the generated partition. This should not happen.']); 0072 end 0073 0074 if thisPartitionData.flagResult == 1 0075 % accept and add partition 0076 partitionData(end+1) = thisPartitionData; 0077 elseif thisPartitionData.flagResult == -1 0078 % partition will be ignored, nothing to do 0079 else 0080 error('elk:decomposition:numericalError', ['flagResult is unknown. ' ... 0081 'This should not happen.']); 0082 end 0083 % fprintf(' part '); 0084 0085 % compute volume 0086 thisHrep = closeCone(thisPartitionData.hrep, dec.centerRay); 0087 thisVolume = computeVolumeHrep(thisHrep); 0088 gatheredVolume = gatheredVolume + thisVolume; 0089 % fprintf(' %f \n', thisVolume); 0090 0091 %% Update remaining domain 0092 % we substract the unified partition, we just found. note that we 0093 % subtract the unclosed cone from the closed remaining domain, which 0094 % is perfectly OK. 0095 remainingDomain = subtractPolytopeHrep(remainingDomain, ... 0096 thisPartitionData.hrep, zerotol); 0097 0098 % if this is empty, we have obviously finished 0099 if isempty(remainingDomain) 0100 break; 0101 end 0102 end 0103 if verbose, fprintf('\b\b\b\b\b\b');end 0104 0105 % recover random number generator 0106 rng(randState); 0107 0108 end