addMeasureDataDec

 add mixed volumes for measure calculation

 Syntax: dec = addMixedVolumeData(dec, measure)
   or    dec = addMixedVolumeData(dec, measure, par)

 This function adds the coefficients that are required to calculate the
   given measure based on a lambda-vector in a partition or a hC-vector.

 The parameters in par are usually directions and contain one parameter
   for each row.
 
 Use obtainMeasureInfo for a listing of measures that are supported.

 addMeasureDataDec(.., 'recalculate') enforces the algorithm to
   recalculate the coefficients, if the measure is already present.

 See also: obtainDec, computeMeasureDec, validateMeasureData, validateDec

• computeMixedVolumeTensor calculate the coefficients for measure calculation in S-representation
• convertMeasureTensorSe convert the tensor of mixed volumes for application to hC-vectors
• obtainMeasureInfo returns information of an available measures in a struct

• addUnifyingMappingData add mapping information
• recomputeMeasureData recompute measure information according to original dec

0001 function dec = addMeasureDataDec(dec, measure, par, varargin)
0002 % add mixed volumes for measure calculation
0003 %
0004 % Syntax: dec = addMixedVolumeData(dec, measure)
0005 %   or    dec = addMixedVolumeData(dec, measure, par)
0006 %
0007 % This function adds the coefficients that are required to calculate the
0008 %   given measure based on a lambda-vector in a partition or a hC-vector.
0009 %
0010 % The parameters in par are usually directions and contain one parameter
0011 %   for each row.
0012 %
0013 % Use obtainMeasureInfo for a listing of measures that are supported.
0014 %
0015 % addMeasureDataDec(.., 'recalculate') enforces the algorithm to
0016 %   recalculate the coefficients, if the measure is already present.
0017 %
0018 % See also: obtainDec, computeMeasureDec, validateMeasureData, validateDec
0019 
0020 % The elk-library: convex geometry applied to crystallization modeling.
0021 %   Copyright (C) 2012 Alexander Reinhold
0022 %
0023 % This program is free software: you can redistribute it and/or modify it
0024 %   under the terms of the GNU General Public License as published by the
0025 %   Free Software Foundation, either version 3 of the License, or (at your
0026 %   option) any later version.
0027 %
0028 % This program is distributed in the hope that it will be useful, but
0029 %   WITHOUT ANY WARRANTY; without even the implied warranty of
0030 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0031 %   General Public License for more details.
0032 %
0033 % You should have received a copy of the GNU General Public License along
0034 %   with this program.  If not, see <http://www.gnu.org/licenses/>.
0035 
0036 %% Input Handling
0037 [recalculateFlag verbose zerotol] = mapOption(varargin, ...
0038     'recalculate', 0, 'verbose', 0, 'zerotol', elkZerotol);
0039 
0040 % make parameter (par) existing
0041 if ~exist('par', 'var')
0042     par = [];
0043 end
0044 
0045 % collect measure information (name, handle, dim, ispar)
0046 if isstruct(measure)
0047     measureStruct = measure;
0048 else
0049     measureStruct = obtainMeasureInfo(measure, dec.A);
0050 end
0051 
0052 %% Create measure
0053 if ~isfield(dec, 'measureData')
0054     dec.measureData = struct();
0055 end
0056 
0057 if ~isfield(dec.measureData, measureStruct.name) || recalculateFlag
0058     dec.measureData.(measureStruct.name) = measureStruct;
0059 end
0060 
0061 %% prepare seMatrix as h-vectors
0062 hMatrix  = (dec.mappingReducedToFull*dec.se')';
0063 %! DEV: This is not nice but this correction is necessary in the
0064 %       scope of projectDec. And it is still valid in the scope of
0065 %       crystal decompositions because the crystals are always
0066 %       centered around zero.
0067 % hMatrix(abs(hMatrix) < sqrt(zerotol)) = 0;
0068 % OK - exactly this line causes a problem for decompositions with fixed
0069 %   positions that WILL have facets close to zero. In particular, for
0070 %   generic lactose, case 5, with the feret diameter in direction [1 1 1].
0071 % Hence, this line is removed and measure calculation for projected
0072 %   decompositions should rather rely on the proper embedding than on a new
0073 %   compilation of the measure data.
0074 %/
0075 
0076 %% Handle tensor calculation
0077 % Collect tensor
0078 if ~measureStruct.isPar
0079     
0080     %dec.measureData.(measureStruct.name).par = [];
0081     measureStruct.par = [];
0082     
0083     if ~isfield(dec.measureData.(measureStruct.name), 'tensorSe')
0084         % compute tensor
0085         tic
0086         measureStruct.tensorSe = computeMixedVolumeTensor(...
0087             hMatrix, measureStruct.handleH, measureStruct.dim, ...
0088             dec.seIncidenceMatrix, verbose);
0089         measureStruct.time.mixedVolumeTensor = toc;   
0090         
0091         % convert tensor
0092         tic
0093         if verbose
0094             disp('Cycling through simplex partitions: ')
0095             fprintf('......');
0096         end
0097         for iPart = 1:length(dec.simplexPartition)
0098             if verbose
0099                 fprintf('\b\b\b\b\b\b%5.1f%%', ...
0100                     (iPart/length(dec.simplexPartition))*100);
0101             end
0102             measureStruct.tensorHc{iPart} = convertMeasureTensorSe(...
0103                 measureStruct.tensorSe, ...
0104                 dec.simplexPartition{iPart}.seIndexVector, ...
0105                 dec.simplexPartition{iPart}.mappingHcToLam, ...
0106                 measureStruct.dim); 
0107         end
0108         if verbose, fprintf('\b\b\b\b\b\b');end
0109         measureStruct.time.convertedTensor = toc;
0110     end
0111     
0112     dec.measureData.(measureStruct.name) = measureStruct;
0113     
0114 else
0115     
0116     % initialize tensorSe
0117     if ~isfield(dec.measureData.(measureStruct.name), 'tensorSe')
0118         dec.measureData.(measureStruct.name).tensorSe = cell(0);
0119     end
0120     
0121     % initialize par
0122     if ~isfield(dec.measureData.(measureStruct.name), 'par')
0123         dec.measureData.(measureStruct.name).par = [];
0124     end
0125     
0126     % shortcut
0127     measureStruct = dec.measureData.(measureStruct.name);
0128     
0129     % we handle parameters subsequently
0130     for iPar = 1:size(par, 1)
0131         % only so something, if measure for parameter does not yet exist
0132         if ~isempty(...
0133                 eliminateRows(par(iPar,:), ...
0134                 measureStruct.par, 0)...
0135                 )
0136             
0137             % add parameter
0138             measureStruct.par(end+1, :) = par(iPar, :);
0139             
0140             % add seTensor
0141             tic
0142             measureStruct.tensorSe{end+1} = ...
0143                 computeMixedVolumeTensor(hMatrix, ...
0144                     @(h)(measureStruct.handleH(h, par(iPar, :))), ...
0145                     measureStruct.dim, ...
0146                     dec.seIncidenceMatrix, verbose);
0147            measureStruct.time(iPar).mixedVolumeTensor = toc;   
0148 
0149            % add tensors for hC
0150            tic
0151            if verbose
0152             disp('Cycling through simplex partitions: ')
0153             fprintf('......');
0154            end
0155            for iPart = 1:length(dec.simplexPartition)
0156                if verbose
0157                    fprintf('\b\b\b\b\b\b%5.1f%%', ...
0158                        (iPart/length(dec.simplexPartition))*100);
0159                end
0160                measureStruct.tensorHc{iPar}{iPart} = convertMeasureTensorSe(...
0161                    measureStruct.tensorSe{end}, ...
0162                    dec.simplexPartition{iPart}.seIndexVector, ...
0163                    dec.simplexPartition{iPart}.mappingHcToLam, ...
0164                    measureStruct.dim);
0165            end
0166            if verbose, fprintf('\b\b\b\b\b\b');end
0167            measureStruct.time(iPar).convertedTensor = toc; 
0168         end
0169     end
0170     
0171     % store shortcut
0172     dec.measureData.(measureStruct.name) = measureStruct;
0173 end
0174 
0175 end
0176