%% 

%{

To run this script, first run axiomsCheckingScriptUZHfinalNEW.m to generate
the required matrices. 

%}


%%


% Only use knife-edge scoring rules, as they are just unions of the other
% ones. Hence, they cover the interval-scoring rules, too.


ruleSet = cat(3, ...
            scoring0Plot(:,profilesMainExperiment), ...
            scoring0to33Plot(:,profilesMainExperiment), ...
            scoring33Plot(:,profilesMainExperiment), ...
            scoring33to50Plot(:,profilesMainExperiment), ...
            scoring50Plot(:,profilesMainExperiment), ...
            scoring50to60Plot(:,profilesMainExperiment), ...
            scoring60Plot(:,profilesMainExperiment), ...
            scoring60to67Plot(:,profilesMainExperiment), ...
            scoring67Plot(:,profilesMainExperiment), ...
            scoring67to75Plot(:,profilesMainExperiment), ...
            scoring75Plot(:,profilesMainExperiment), ...
            scoring75to80Plot(:,profilesMainExperiment), ...
            scoring80Plot(:,profilesMainExperiment), ...
            scoring80to100Plot(:,profilesMainExperiment), ...
            scoring100Plot(:,profilesMainExperiment), ...
            condorcetPlot(:,profilesMainExperiment), ...
            malevolentScore0Plot(:,profilesMainExperiment), ...
            malevolentScore0to33Plot(:,profilesMainExperiment), ...
            malevolentScore33Plot(:,profilesMainExperiment), ...
            malevolentScore33to50Plot(:,profilesMainExperiment), ...
            malevolentScore50Plot(:,profilesMainExperiment), ...
            malevolentScore50to60Plot(:,profilesMainExperiment), ...
            malevolentScore60Plot(:,profilesMainExperiment), ...
            malevolentScore60to67Plot(:,profilesMainExperiment), ...
            malevolentScore67Plot(:,profilesMainExperiment), ...
            malevolentScore67to75Plot(:,profilesMainExperiment), ...
            malevolentScore75Plot(:,profilesMainExperiment), ...
            malevolentScore75to80Plot(:,profilesMainExperiment), ...
            malevolentScore80Plot(:,profilesMainExperiment), ...
            malevolentScore80to100Plot(:,profilesMainExperiment), ...
            malevolentScore100Plot(:,profilesMainExperiment), ...
            badCondorcetPlot(:,profilesMainExperiment)); 
            
            %scoringRunoff0to33Plot(:,profilesMainExperiment), ...    
            %scoringRunoff33to50Plot(:,profilesMainExperiment), ...
            %scoringRunoff50to100Plot(:,profilesMainExperiment), ...
            %badScoringRunoff0to33Plot(:,profilesMainExperiment), ...
            %badScoringRunoff33to50Plot(:,profilesMainExperiment), ...
            %badScoringRunoff50to60Plot(:,profilesMainExperiment)
            
mult = sum(ruleSet, 1);
powerOne = squeeze(sum(mult==1));       
powerTwo = squeeze(sum(mult==2));       
powerThree = squeeze(sum(mult==3));   

numExpansions = 2.^powerTwo.*3.^powerThree;
totalExpansions = sum(numExpansions); 

[X, ruleLength, numExpandRules] = size(ruleSet);

ruleLabel = [1:numExpandRules]';
knifeEdgeScoring = [[1 0 1 0 1 0 1 0 1 0 1 0 1 0 1], zeros(1, numExpandRules-15)]';





r = 1;
while r <= numExpandRules 
    
    [X, ruleLength, numExpandRules] = size(ruleSet);

    rule = ruleSet(:,:,r);
    ruleLabelCurrent = ruleLabel(r);
    knifeEdgeCurrent = knifeEdgeScoring(r);
    
    multiplicity = sum(rule, 1);
    
    % Advance to the next rule if the current rule is exhausted
    if max(multiplicity) == 1
        r = r+1;
    else
        
        % Select the index of the first multiplicity
        [m, splitIndex] = max(multiplicity); 

        subrule = zeros(0, 0);
        subrule_j = zeros(0, 0);
        for j = 1 : 3
            if rule(j, splitIndex) == 1
                subrule_j = [rule(:, 1:splitIndex-1), [+(j==1);+(j==2);+(j==3)], rule(:,splitIndex+1:ruleLength)];
                if numel(subrule_j) ~= 0
                    if numel(subrule) == 0 
                        subrule = subrule_j;
                    else
                        subrule = cat(3, subrule, subrule_j);
                    end
                end
            end
        end
        [X, Y, expansionSize] = size(subrule);
        ruleLabelExpansion = ruleLabelCurrent * ones(expansionSize,1);
        knifeEdgeExpansion = knifeEdgeCurrent * ones(expansionSize,1);
        % Replace rule under consideration by expanded rule
        ruleSet = cat(3, ruleSet(:,:,1:r-1), subrule, ruleSet(:,:,r+1:numExpandRules));

        ruleLabel = [ruleLabel(1:r-1); ruleLabelExpansion; ruleLabel(r+1:numExpandRules)];
        knifeEdgeScoring = [knifeEdgeScoring(1:r-1); knifeEdgeExpansion; knifeEdgeScoring(r+1:numExpandRules)];
        
    end
    
end

%% Find and remove duplicate rules

[~,~,nRules] = size(ruleSet);

identicalRules = zeros(nRules);
dropRule = zeros(nRules, 1);


for i = 1 : nRules
    for j = i + 1 : nRules 
        identicalRules(i,j) = min(ruleSet(:,:,i) == ruleSet(:,:,j), [], 'all');
        if identicalRules(i,j) == 1 
            if knifeEdgeScoring(i) == 0 && knifeEdgeScoring(j) == 1
                dropRule(j) = 1;
            elseif knifeEdgeScoring(i) == 1 && knifeEdgeScoring(j) == 0
                dropRule(i) = 1;
            else
                dropRule(j) = 1;
            end
        end
    end
end


keepRule = logical(1-dropRule);
ruleSet = ruleSet(:,:,keepRule);
ruleLabel = ruleLabel(keepRule);

[numExpandRules, ~] = size(ruleLabel);

disp('Number of resolute components')
disp(numExpandRules)

%% Export expanded rule data

prefsRuleExpansion = pref(:, profilesMainExperiment, :);
profileIndexExpansion = profileIndex(profilesMainExperiment);
nRankingsExpanded = sum(profilesMainExperiment);
prefExportExpansion = [prefsRuleExpansion(:,:,1)', prefsRuleExpansion(:,:,2)', prefsRuleExpansion(:,:,3)', prefsRuleExpansion(:,:,4)', prefsRuleExpansion(:,:,5)'];

expansionNumber = 1;
ruleResultsLong = [profileIndexExpansion, ones(nRankingsExpanded, 1), ones(nRankingsExpanded, 1), ruleSet(:,:,1)', prefExportExpansion];
for r = 2 : numExpandRules 
    expansionNumber = expansionNumber + 1;
    if ruleLabel(r) ~= ruleLabel(r-1)
        expansionNumber = 1;
    end
    ruleResultsLong = cat(1, ruleResultsLong, [profileIndexExpansion, ruleLabel(r)*ones(nRankingsExpanded, 1), expansionNumber*ones(nRankingsExpanded, 1), ruleSet(:,:,r)', prefExportExpansion]);
end


rowNames = 1:nRankings;
colNames = {'profileIndex', 'ruleIndex', 'expansionIndex', 'optionA', 'optionB', 'optionC', 'pref1st_1', 'pref2nd_1', 'pref3rd_1', 'pref1st_2', 'pref2nd_2', 'pref3rd_2', 'pref1st_3', 'pref2nd_3', 'pref3rd_3', 'pref1st_4', 'pref2nd_4', 'pref3rd_4', 'pref1st_5', 'pref2nd_5', 'pref3rd_5'};
table = array2table(ruleResultsLong,'VariableNames',colNames);
writetable(table, 'data/rulePredictions3optExpanded.csv')

%% Export rule labels data


save('data/ruleLabelsExpanded', 'ruleLabel');

load('data/ruleLabelsExpanded');





