% %% Generating Structured Representatations and Measures

% data is hopefuly self-explanatory

%% File
SRfile = 'data/data.mat';

%% Options
toleranceDecomposition = 1e-10;
toleranceMeasure = 1e-6;
verboseDecomposition = 2;

%% Switches
% normal option does generate the data, if not already existent (it
% enforces 1 for all preceeding options, as well)
% redo options eliminate any calculation afterwards, which are
% automatically generated again, if normal option is set to one
flagCleanAll      = 0;

flagRmLoadCrystal = 0; 
flagLoadCrystal   = 1; % load and generate new data

flagRmDecomposition = 0;
flagDecomposition   = 1; % generate partitioning

flagRmValidateDecomposition    = 0;
flagValidateDecomposition      = 1; % verify partitioning

flagRmMeasure    = 0;
flagRmAllMeasure = 0; % reset all measure stuff
flagMeasure      = 1; % generate measures

flagRmVerifyMeasure = 0;
flagVerifyMeasure   = 1; % verify that all simplex partitions are identic in hC
% 
flagRmMeasTest = 0;
flagMeasTest   = 1;

% % enforce 1 for preceding options:
flagMeasure = flagMeasure || flagVerifyMeasure || flagMeasTest;
% empty: measure generation does not require verification
flagDecomposition = flagDecomposition || flagValidateDecomposition || flagMeasure;
flagLoadCrystal = flagLoadCrystal || flagDecomposition;


% 
% % enforce rm 1 for subsequent options
flagRmDecomposition = flagRmDecomposition || flagRmLoadCrystal;
flagRmValidateDecomposition = flagRmValidateDecomposition || flagRmDecomposition;
flagRmMeasure = flagRmMeasure || flagRmDecomposition; 
% ^- don't remake measures if verifying again
flagRmVerifyMeasure = flagRmVerifyMeasure || flagRmMeasure || flagRmAllMeasure;
flagRmMeasTest = flagRmMeasTest || flagRmMeasure;

%% Processing affects the following set of examples
% code to generate index combinations:
%   n=7; k=2; sort(sum(2.^(obtainOrderedSetMatrix(n,k)-1), 2)')
cdefList = cell(0);
cdefCaseList = cell(0);

% STATUS: 5D complete
% ERRORS: 
% PROBLEMS: 
% cdefList{end+1} = 'generic_potash_alum'; % 1 + max 4; par: max 15 (=2^2-1)
% cdefCaseList{end+1}  = 0; % 1D (1 case)
% cdefCaseList{end}  = [cdefCaseList{end} 1 2 4 8]; % 2D (4 cases)
% cdefCaseList{end}  = [cdefCaseList{end} 3 5 9 6 10 12]; % 3D (6 cases)
% cdefCaseList{end}  = [cdefCaseList{end} 7 11 13 14]; % 4D (4 cases)
% cdefCaseList{end}  = [cdefCaseList{end} 15]; % 5D (1 case)
% testing/debugging
% cdefCaseList{end+1} =  3;

% STATUS: 7D complete
% ERRORS:
% PROBLEMS: 15 23 27 currently while DVD viewing
% cdefList{end+1} = 'generic_paracetamol'; % 3 + max 7; par: max 127 (=2^7-1)
% cdefCaseList{end+1} = 0; % 3D (1 case)
% cdefCaseList{end}   = [cdefCaseList{end} 1 2 4 8 16 32 64]; % 4D (7 cases)
% cdefCaseList{end}   = [cdefCaseList{end} 3 5 6 9 10 12]; % 5D (7 cases)
% cdefCaseList{end}   = [cdefCaseList{end} 7 11 13 14 19 21]; % 6D (5 cases)
% cdefCaseList{end}   = [cdefCaseList{end} 15 23 27]; % 7D (3)
% cdefCaseList{end}   = [cdefCaseList{end} 31]; % 8D (1)
%cdefCaseList{end}   = []; % 9D
%cdefCaseList{end}   = []; % 10D
%cdefCaseList{end}   = []; % 11D
% testing/debugging
% cdefCaseList{end} = 19;

% STATUS: 7D complete
% ERRORS: 61 (verification conversion Hrep to Vrep)
% PROBLEMS: 
cdefList{end+1} = 'generic_MgSO4_11H2O'; % 3 + max 10; par: max 1023 (=2^2-1)
cdefCaseList{end+1} = 0; % 3D (1 case)
cdefCaseList{end}   = [cdefCaseList{end} 1 2 4 8 16 32 64 128 256 512]; % 4D (8 cases)
% cdefCaseList{end}   = [cdefCaseList{end} 3 5 6 9 10 12 17 18 20 24]; % 5 (10)
% cdefCaseList{end}   = [cdefCaseList{end} 7 11 13 14 19 21 25 26 28 35 37 38 41 42]; % 6 (14)
% cdefCaseList{end}   = [cdefCaseList{end} 15 23 27 29 30 39 43 45 46 51 53 54 57 58 60]; % 7 (15) %PROBLEM: 60 15
% cdefCaseList{end}   = [cdefCaseList{end} 31 47 55 59 61 62 79 87]; % 8 (8) 
%                      47 (61, 62, 55, 87, 79, not verified) omitted (31, 61, 62, 79, 87 delaunayn problem volume)
% cdefCaseList{end}   = [cdefCaseList{end} ]; % 9 (2) % 95 omitted 119 (ja) 123 111 63(nein)
%cdefCaseList{end}   = [cdefCaseList{end} 127]; % 10
% testing/debugging
% cdefCaseList{end} = 61;

%% Some setup
clc

%% Initialize data structure & clean_all-option

% delete data file
if flagCleanAll && exist(SRfile, 'file') 
    delete(SRfile);
end

% load data / prepare as struct
if ~exist(SRfile, 'file');
    data = struct;
else
    data = load(SRfile);
    data = data.data;
end

caseNameList = cell(0);
for i = 1:length(cdefList)
    if isempty(cdefCaseList{i})
        if ~isfield(data, cdefList{i})
            data.(cdefList{i}).filename = cdefList{i};
            data.(cdefList{i}).par = [];
        end
        caseNameList{end+1} = cdefList{i};
    else
        for j = 1:length(cdefCaseList{i})
            if ~isfield(data, [cdefList{i} '_' num2str(cdefCaseList{i}(j))])
                data.([cdefList{i} '_' num2str(cdefCaseList{i}(j))]).filename =...
                    cdefList{i};
                data.([cdefList{i} '_' num2str(cdefCaseList{i}(j))]).par = ...
                    cdefCaseList{i}(j);
            end
            caseNameList{end+1} = [cdefList{i} '_' num2str(cdefCaseList{i}(j))];
        end
    end
end

%% Crystal Data
for i = 1:length(caseNameList)
    if flagRmLoadCrystal && isfield(data.(caseNameList{i}), 'cdef')
        data.(caseNameList{i}) = rmfield(data.(caseNameList{i}), 'cdef');
    end
    if flagLoadCrystal && ~isfield(data.(caseNameList{i}), 'cdef')
        tic
        disp(['Loading ' caseNameList{i} ' ..']);
        data.(caseNameList{i}).cdef = obtainCrystal(...
            data.(caseNameList{i}).filename, ...
            data.(caseNameList{i}).par);
        data.(caseNameList{i}).time.loadCdef = toc;
        fprintf('\b done.\n');
    end
end

%% Decomposition
for i = 1:length(caseNameList)
    if flagRmDecomposition && isfield(data.(caseNameList{i}), 'dec')
        data.(caseNameList{i}) = rmfield(data.(caseNameList{i}), 'dec');
    end
    if flagDecomposition && ~isfield(data.(caseNameList{i}), 'dec')
        tic
        fprintf(['Generating Domain for ' caseNameList{i} ' ..\n']);
        data.(caseNameList{i}).dec = obtainDec(...
            data.(caseNameList{i}).cdef.A, 'mapping', ...
            data.(caseNameList{i}).cdef.groupMappingMatrix, ...
            'verbose', verboseDecomposition);
        data.(caseNameList{i}).time.obtainDec = toc;
        fprintf('\b done.\n');
    end
end

%% Verify Decomposition
for i = 1:length(caseNameList)
    if flagRmValidateDecomposition && isfield(data.(caseNameList{i}), 'decVerified')
        data.(caseNameList{i}) = rmfield(data.(caseNameList{i}), 'decVerified');
        data.(caseNameList{i}) = rmfield(data.(caseNameList{i}), 'decTolerance');
    end
    if flagValidateDecomposition && ~isfield(data.(caseNameList{i}), 'decVerified')
        fprintf(['Validating decomposition for ' caseNameList{i} ' .. \n']);
        tic
        data.(caseNameList{i}).decTolerance = toleranceDecomposition;
        [data.(caseNameList{i}).dec data.(caseNameList{i}).decVerified] = ...
            validateDec(data.(caseNameList{i}).dec, ...
              'verbose', 1, ...
              'zerotol', toleranceDecomposition, ...
              'output', 'both');
        data.(caseNameList{i}).time.validateDec = toc;
        if data.(caseNameList{i}).decVerified
            fprintf('  valid.\n');
        else
            fprintf('  NOT VALID.\n');
        end
    end
end

%% Add Mixed Volumes & Conversion & Verifying (Simplex <-> Unified) &
%% Measure Testing

% This handles FOUR seperate steps!
for i = 1:length(caseNameList)
if isfield(data.(caseNameList{i}), 'dec')
    
    % load empty or already generated decomposition for measures
    if flagRmAllMeasure || ~isfield(data.(caseNameList{i}), 'decMeas')
        decMeas = data.(caseNameList{i}).dec;
%         decMeas.meas = struct;
%         data.(caseNameList{i}).meas_name = cell(0);
%         data.(caseNameList{i}).time_meas1 = [];
%         data.(caseNameList{i}).time_meas2 = [];
%         data.(caseNameList{i}).time_meas3 = [];
%         data.(caseNameList{i}).time_meas4 = [];
%         data.(caseNameList{i}).vermeas    = [];
%         data.(caseNameList{i}).meastest   = testStructRepMeas_hC(decMeas, 'measures', {}, 'summarized', 1);
    else
        decMeas = data.(caseNameList{i}).decMeas;
    end
    
    % make sure the field measureData exists in decMeas
    if ~isfield(decMeas, 'measureData')
        decMeas.measureData = struct();
    end
    
    for j = 1:length(measureList)
%         % search, if measure already exists, when (where) it was added
%         if isfield(decMeas.meas, measureList(j).name)
%             % increase search index so that it throws error, if not found
%             for k = 1:(length(data.(caseNameList{i}).meas_name)+1)
%                 if strcmp(data.(caseNameList{i}).meas_name{k}, measureList(j).name)
%                     break;
%                 end
%             end
%             % k is now the appropriate index
%         else
%             k = -1;
%         end

        % rm has effect on each measure seperately
        if flagRmMeasure && isfield(decMeas.measureData, measureList(j).name)
            decMeas.measureData = rmfield(decMeas.measureData, measureList(j).name);
        end
        
        % generating measure
        if flagMeasure && ~isfield(decMeas.measureData, measureList(j).name)
            fprintf(['Calculating Measure ' measureList(j).name ' for ' caseNameList{i} ' ..\n']);
            decMeas = addMeasureDataDec(decMeas, measureList(j).name, ...
                measureList(j).par);
            fprintf('  done.\n');
        end
        
        % rm verify measure
        if flagRmVerifyMeasure && isfield(decMeas.measureData.(measureList(j).name), ...
                'verifyMatrix')
            decMeas.measureData.(measureList(j).name) = ...
                rmfield(decMeas.measureData.(measureList(j).name), ...
                'verifyMatrix');
            decMeas.measureData.(measureList(j).name) = ...
                rmfield(decMeas.measureData.(measureList(j).name), ...
                'verifyTolerance');
        end
        
        % verifying measure
        if flagVerifyMeasure && ...
                ~isfield(decMeas.measureData.(measureList(j).name), 'verifyMatrix')
            fprintf(['Validating measure ' measureList(j).name ' for '...
                     caseNameList{i} ' ..\n']);
            decMeas.measureData.(measureList(j).name).verifyTolerance = ...
                toleranceMeasure;
            tic
            decMeas.measureData.(measureList(j).name).verifyMatrix = ...
                validateMeasureDataDec(decMeas, 'detailed', 1, ...
                  'zerotol', toleranceMeasure, 'measure', measureList(j).name);
            decMeas.timing.validateMeasureData = toc;
            if all(decMeas.measureData.(measureList(j).name).verifyMatrix(:))
                fprintf('  valid.\n');
            else
                fprintf('  NOT VALID.\n');
            end
        end
        % rm measure testing
        if flagRmMeasTest && isfield(decMeas.measureData.(measureList(j).name), 'timing')
            decMeas.measureData.(measureList(j).name) = ...
                rmfield(decMeas.measureData.(measureList(j).name), 'timing');
        end
        % testing measure
        % The goal here is to take at least minumumTestNumber(1) hC-vectors
        %   equally distributed over the simplex partitions for the
        %   computation of the partition. Then, only minimumTestNumber(2)
        %   preferably equally distributed vectors shall be used to get the
        %   timing for measure calculation. minimumTestNumber(1) is usually
        %   much larger than minimumTestNumber(2).
        if flagMeasTest && ~isfield(decMeas.measureData.(measureList(j).name), 'timing')
            % screen message
            fprintf(['Testing Measure ' measureList(j).name ' for ' caseNameList{i} ' ..\n']);
            
            % ensure that at least minimumTestNumber(1) are generated
            sampleNumberPerSimplexPartition = ...
                ceil(measureList(j).minimumTestNumber(1) / decMeas.nPS);
            % generate hC-vectors per simplex partition
            sampleHcMatrix = zeros([...
                decMeas.nPS * sampleNumberPerSimplexPartition, ...
                decMeas.nC]);
            for iSimplex = 1:decMeas.nPS
                sampleHcMatrix(...
                    (iSimplex-1)*sampleNumberPerSimplexPartition + ...
                    (1:sampleNumberPerSimplexPartition), :) = ...
                    rand([sampleNumberPerSimplexPartition, decMeas.nC]) * ...
                      decMeas.se(...
                        decMeas.simplexPartition{iSimplex}.seIndexVector, ...
                        :);
            end
            % obtain unified partitions
            tic
            thisSimplexVector = computePartitionIndexDec(decMeas, ...
                sampleHcMatrix, [], 1);
            decMeas.measureData.(measureList(j).name).timing.obtainUnifiedTime = toc;
            decMeas.measureData.(measureList(j).name).timing.obtainUnifiedSamples = ...
                size(sampleHcMatrix, 1);
            % obtain simplex partitions
            tic
            thisSimplexVector = computePartitionIndexDec(decMeas, ...
                sampleHcMatrix, [], 0);
            decMeas.measureData.(measureList(j).name).timing.obtainSimplexTime = toc;
            decMeas.measureData.(measureList(j).name).timing.obtainSimplexSamples = ...
                size(sampleHcMatrix, 1);
            
            downsampleRate = max(1, ...
                                 fix(measureList(j).minimumTestNumber(1) / ...
                                     measureList(j).minimumTestNumber(2))...
                                );
            sampleHcMatrix = sampleHcMatrix(1:downsampleRate:end, :);
            thisSimplexVector = thisSimplexVector(1:downsampleRate:end);
            
            % obtain measure
            if isempty(measureList(j).par)
                nPar = 1;
                tic
                tmp = computeMeasureDec(decMeas, 'hC', ...
                    {thisSimplexVector, sampleHcMatrix}, ...
                    measureList(j).name, ...
                    [], 1);
            else
                nPar = size(measureList(j).par, 1);
                tic
                tmp = computeMeasureDec(decMeas, 'hC', ...
                    {thisSimplexVector, sampleHcMatrix}, ...
                    {measureList(j).name, measureList(j).par}, ...
                    [], 1);
            end
            decMeas.measureData.(measureList(j).name).timing.obtainMeasureTime = toc;
            decMeas.measureData.(measureList(j).name).timing.obtainMeasureSamples = ...
                size(sampleHcMatrix, 1);
            decMeas.measureData.(measureList(j).name).timing.obtainMeasureNPar = ...
                nPar;
            
            % obtain real measure
            downsampleRate = max(1, ...
                                 fix(measureList(j).minimumTestNumber(2) / ...
                                     measureList(j).minimumTestNumber(3))...
                                );
            sampleHcMatrix = sampleHcMatrix(1:downsampleRate:end, :);
            thisSimplexVector = thisSimplexVector(1:downsampleRate:end);
            sampleHMatrix = sampleHcMatrix * decMeas.mappingReducedToFull';
            
            if isempty(measureList(j).par)
                nPar = 1;
                tic
                for iSample = 1:size(sampleHcMatrix, 1)
                    tmp = decMeas.measureData.(measureList(j).name).handleHrep(...
                            struct('A', decMeas.A, 'h', sampleHMatrix(iSample, :)')...
                            );
                end
            else
                nPar = size(measureList(j).par, 1);
                tic
                for iSample = 1:size(sampleHcMatrix, 1)
                    tmp = decMeas.measureData.(measureList(j).name).handleHrep(...
                            struct('A', decMeas.A, 'h', sampleHMatrix(iSample, :)'), ...
                            measureList(j).par);
                end
            end
            decMeas.measureData.(measureList(j).name).timing.realMeasureTime = toc;
            decMeas.measureData.(measureList(j).name).timing.realMeasureSamples = ...
                size(sampleHcMatrix, 1);
            decMeas.measureData.(measureList(j).name).timing.realMeasureNPar = ...
                nPar;
            
            
            fprintf('  done.\n');
%             tic
%             hCtests = max(1, round(30/...
%                 size(data.(caseNameList{i}).SR_part.simplex.se_partitions, 1)));
%             data.(caseNameList{i}).meastest(k) = testStructRepMeas_hC(decMeas, ...
%                 'measures', measureList(j).name, 'verbose', 0, 'summarized', 1, ...
%                 'tests', hCtests, 'extra_tests', measureList(j).extra_tests);
%             %data.(example{i}).meastest(k).measure = meas(j).name;
%             data.(caseNameList{i}).time_meas4(k) = toc;
%             relmin = data.(caseNameList{i}).meastest(k).relmin;
%             relmax = data.(caseNameList{i}).meastest(k).relmax;
%             if max(abs(relmin),abs(relmax)) < 1e-10, tag='[+]';else tag='[-]';end
%             fprintf([' done. \n' tag ' (min %.2e, max %.2e)'], relmin, relmax);
        end
    end

    % write shortcut back into data structure
    data.(caseNameList{i}).decMeas = decMeas;
end
end
fprintf('\n');


%% Save data
save(SRfile, 'data');
