% Arrow project

clear
clc

recalculateUniqueProfiles = 1;

%{

EXPLANATION OF PREFERENCE PROFILES

A preference profile is a 3-by-5 matrix of integers 1 to 3

Example: 

A B C B B
B A B C A
C C A A C

Each column stands for one voter. 
Each line stands for a preference rank. 

In the example above, voter 5 ranks option C first, option A second, and
option B last. 

Because matlab more easily deals with numbers than letters, the options are
denotedf by numbers, according to the following convention: 

A = 1
B = 2 
C = 3

Accordingly, the above profile is represented as: 

1 2 3 2 2
2 1 2 3 1
3 3 1 1 3

That is, each integer in the matrix represents an option, it does NOT
represent a preference rank. Preference ranks are indicated by the position
of an integer, not by its value. 

(The above example can easily be misinterpreted as showing that voter 4
puts option A in second place, which is not true). 

%}





%% Generate data

nRankings = 126;
nCitizens = 5; 
nOptions = 3;


%% Generate all preference profiles with 5 voters

profile = ones(1, 5);

index = 0;

for i2 = 1 : 6
    for i3 = i2 : 6
        for i4 = i3 : 6
            for i5 = i4 : 6
               index = index + 1; 
               profile(index,:) = [1, i2, i3, i4, i5];
            end
        end
    end
end

% Expand profile indexes into preference profiles 

[nProfiles, nVoters] = size(profile);

profileIndices = profile;

profile1 = profileIndices == 1;
profile2 = profileIndices == 2;
profile3 = profileIndices == 3;
profile4 = profileIndices == 4;
profile5 = profileIndices == 5;
profile6 = profileIndices == 6;

profileCoordinates = zeros(nRankings, 6);
profileCoordinates(:, 1) = sum(profile1, 2);
profileCoordinates(:, 2) = sum(profile2, 2);
profileCoordinates(:, 3) = sum(profile3, 2);
profileCoordinates(:, 4) = sum(profile4, 2);
profileCoordinates(:, 5) = sum(profile5, 2);
profileCoordinates(:, 6) = sum(profile6, 2);




profile = reshape(profile, 1, nProfiles, nVoters);
profile = repmat(profile, 3, 1, 1);

order1 = [1; 2; 3];
order2 = [1; 3; 2];
order3 = [3; 1; 2];
order4 = [3; 2; 1];
order5 = [2; 3; 1];
order6 = [2; 1; 3];

order = [order1, order2, order3, order4, order5, order6];

for profileNr = 1 : nProfiles 
    for voterNr = 1 : nVoters
        profileIndex = profile(1, profileNr, voterNr);
        profile(:, profileNr, voterNr) = order(:, profileIndex);
    end
end

pref = profile;


savename='data/uniqueProfiles.mat';


if recalculateUniqueProfiles == 1
    uniqueProfiles = profilePermutations(pref);
    uniqueProfiles = logical(uniqueProfiles);
    save(savename,'uniqueProfiles')
end

load(savename)

pref = pref(:, uniqueProfiles, :);
profileCoordinates = profileCoordinates(uniqueProfiles, :);

for i = 1:42
    disp(i)
    disp(squeeze(pref(:,i,:)))
end



nRankings = sum(uniqueProfiles);
nProfiles = nRankings;

% Implied binary preference relations

pref12 = pref;
pref13 = pref;
pref23 = pref;

pref12((pref12 ~= 1 & pref12 ~= 2 )) = [];
pref13((pref13 ~= 1 & pref13 ~= 3 )) = [];
pref23((pref23 ~= 2 & pref23 ~= 3 )) = [];

pref12 = reshape(pref12, 2, nRankings, nCitizens);
pref13 = reshape(pref13, 2, nRankings, nCitizens);
pref23 = reshape(pref23, 2, nRankings, nCitizens);


%% Get rank order statistics of the options for each profile


locationA = pref(:,:,:) == 1;
locationB = pref(:,:,:) == 2;
locationC = pref(:,:,:) == 3;

orderStatsA = sum(locationA, 3);
orderStatsB = sum(locationB, 3);
orderStatsC = sum(locationC, 3);

orderStatsAcum = cumsum(orderStatsA, 1);
orderStatsBcum = cumsum(orderStatsB, 1);
orderStatsCcum = cumsum(orderStatsC, 1);

AdomB = orderStatsAcum >= orderStatsBcum;
BdomA = orderStatsBcum >= orderStatsAcum;
BdomC = orderStatsBcum >= orderStatsCcum;
CdomB = orderStatsCcum >= orderStatsBcum;
AdomC = orderStatsAcum >= orderStatsCcum;
CdomA = orderStatsCcum >= orderStatsAcum;

AdomB = min(AdomB, [], 1);
BdomA = min(BdomA, [], 1);
BdomC = min(BdomC, [], 1);
CdomB = min(CdomB, [], 1);
AdomC = min(AdomC, [], 1);
CdomA = min(CdomA, [], 1);

domAB = max(AdomB, BdomA);
domBC = max(BdomC, CdomB);
domAC = max(AdomC, CdomA);

dominanceOrderingFull = min([domAB; domBC; domAC], [], 1);


disp([[1:42]', AdomB', BdomA', BdomC', CdomB', AdomC', CdomA'])
disp([[1:42]', dominanceOrderingFull'])





pointsScore = zeros(nRankings * nCitizens, nOptions);

[pointsScore(:, 1), ~] = find(pref == 1);
[pointsScore(:, 2), ~] = find(pref == 2);
[pointsScore(:, 3), ~] = find(pref == 3);

pointsScore = pointsScore';
pointsScore = reshape(pointsScore, nOptions, nRankings, nCitizens);

AstatewiseDominatesB = pointsScore(1, :, :) < pointsScore(2, :, :);
BstatewiseDominatesA = pointsScore(2, :, :) < pointsScore(1, :, :);
BstatewiseDominatesC = pointsScore(2, :, :) < pointsScore(3, :, :);
CstatewiseDominatesB = pointsScore(3, :, :) < pointsScore(2, :, :);
AstatewiseDominatesC = pointsScore(1, :, :) < pointsScore(3, :, :);
CstatewiseDominatesA = pointsScore(3, :, :) < pointsScore(1, :, :);

ABstatewiseDominateC = min(pointsScore(1, :, :), pointsScore(2, :, :)) < pointsScore(3, :, :);
ACstatewiseDominateB = min(pointsScore(1, :, :), pointsScore(3, :, :)) < pointsScore(2, :, :);
BCstatewiseDominateA = min(pointsScore(2, :, :), pointsScore(3, :, :)) < pointsScore(1, :, :);

AstatewiseDominatesB = min(AstatewiseDominatesB , [], 3);
BstatewiseDominatesA = min(BstatewiseDominatesA , [], 3);
BstatewiseDominatesC = min(BstatewiseDominatesC , [], 3);
CstatewiseDominatesB = min(CstatewiseDominatesB , [], 3);
AstatewiseDominatesC = min(AstatewiseDominatesC , [], 3);
CstatewiseDominatesA = min(CstatewiseDominatesA , [], 3);

ABstatewiseDominateC = min(ABstatewiseDominateC , [], 3);
ACstatewiseDominateB = min(ACstatewiseDominateB , [], 3);
BCstatewiseDominateA = min(BCstatewiseDominateA , [], 3);

swDomAB = max(AstatewiseDominatesB, BstatewiseDominatesA);
swDomBC = max(BstatewiseDominatesC, CstatewiseDominatesB);
swDomAC = max(AstatewiseDominatesC, CstatewiseDominatesA);

swDominance = max([swDomAB', swDomBC', swDomAC'], [], 2);

swPairDominance = max([ABstatewiseDominateC', ACstatewiseDominateB', BCstatewiseDominateA'], [], 2);

disp([[1:42]', swDominance, swPairDominance])




%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%      Get outcomes from specific rules           %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% %%%%%%%%%%%%%%%%%%%%%%%     Scoring rules           %%%%%%%%%%%%%%%%%%%%


%% Get Borda-count outcome

pointsScore = zeros(nRankings * nCitizens, nOptions);

% The variables pointScore are how many points options A, B, C receive
[pointsScore(:, 1), ~] = find(pref == 1);
[pointsScore(:, 2), ~] = find(pref == 2);
[pointsScore(:, 3), ~] = find(pref == 3);

pointsScore = pointsScore';
pointsScore = reshape(pointsScore, nOptions, nRankings, nCitizens);

bordaCount = sum(pointsScore, 3);
[minBordaCount, bordaWinner] = min(bordaCount, [], 1);

winSetBorda = repmat(minBordaCount, nOptions, 1) == bordaCount;
winNumBorda = sum(winSetBorda, 1);

bordaWin1 = bordaCount(1,:) == minBordaCount;
bordaWin2 = bordaCount(2,:) == minBordaCount;
bordaWin3 = bordaCount(3,:) == minBordaCount;


tieBorda = winNumBorda > 1; 

bordaWinner = max(bordaWinner - nOptions * tieBorda, -1);





%% Get plurality outcome

weights = [1, 0, 0];
pluralityScore = weights(1) * (pointsScore == 1) ...
               + weights(2) * (pointsScore == 2) ...
               + weights(3) * (pointsScore == 3) ;
              

pluralityCount = sum(pluralityScore, 3);
[maxPluralityCount, pluralityWinner] = max(pluralityCount, [], 1);



pluralityWin1 = maxPluralityCount == pluralityCount(1, :);
pluralityWin2 = maxPluralityCount == pluralityCount(2, :);
pluralityWin3 = maxPluralityCount == pluralityCount(3, :);

winSetPlurality = repmat(maxPluralityCount, nOptions, 1) == pluralityCount;
winNumPlurality = sum(winSetPlurality, 1);

tiePlurality = winNumPlurality > 1; 
pluralityWinner = max(pluralityWinner - nOptions * tiePlurality, -1);


%% Get general scoring outcome

nScoringRules = 101;
s = linspace(0, 1, nScoringRules);

generalScore = zeros(3, nRankings, nCitizens, nScoringRules);

generalScoreWinner = zeros(nRankings, nScoringRules);
generalScoreWin1 = zeros(nRankings, nScoringRules);
generalScoreWin2 = zeros(nRankings, nScoringRules);
generalScoreWin3 = zeros(nRankings, nScoringRules);

tieGeneralScore = zeros(nRankings, nScoringRules);

for sNr = 1:nScoringRules 
    weights = [1, s(sNr), 0];
    generalScore(:,:,:,sNr) =   weights(1) * (pointsScore == 1) ...
                   + weights(2) * (pointsScore == 2) ...
                   + weights(3) * (pointsScore == 3) ;


    generalScoreCount = sum(generalScore(:,:,:,sNr), 3);
    [maxGeneralScore, generalScoreWinnerS] = max(generalScoreCount, [], 1);

    tol = 1e-10;
    generalScoreWin1(:, sNr) = generalScoreCount(1,:) >= maxGeneralScore - tol;
    generalScoreWin2(:, sNr) = generalScoreCount(2,:) >= maxGeneralScore - tol;
    generalScoreWin3(:, sNr) = generalScoreCount(3,:) >= maxGeneralScore - tol;

    winSetGeneralScore = repmat(maxGeneralScore, nOptions, 1) == generalScoreCount;
    winNumGeneralScore = sum(winSetGeneralScore, 1);

    tieGeneralScore(:,sNr) = [winNumGeneralScore > 1]'; 

    generalScoreWinner(:,sNr) = max(generalScoreWinnerS - nOptions * tieGeneralScore(:,sNr)', -1);
end

s = repmat(s, nRankings, 1);

unavailable1 = generalScoreWin1;
unavailable2 = generalScoreWin2;
unavailable3 = generalScoreWin3;

unavailable1Shift = [unavailable1(:, 2:nScoringRules), unavailable1(:,nScoringRules)];
unavailable2Shift = [unavailable2(:, 2:nScoringRules), unavailable2(:,nScoringRules)];
unavailable3Shift = [unavailable3(:, 2:nScoringRules), unavailable3(:,nScoringRules)];


availabilitySwitch1 = +(unavailable1 ~= unavailable1Shift);
availabilitySwitch2 = +(unavailable2 ~= unavailable2Shift);
availabilitySwitch3 = +(unavailable3 ~= unavailable3Shift);



availabilitySwitch1 = s .* availabilitySwitch1;
availabilitySwitch2 = s .* availabilitySwitch2;
availabilitySwitch3 = s .* availabilitySwitch3;


availabilitySwitchOverall = cat(3, availabilitySwitch1, availabilitySwitch2, availabilitySwitch3);

availabilitySwitchOverall = max(availabilitySwitchOverall, [], 3);

switchingValues = 0 * availabilitySwitchOverall;
numSwitchingValues = 1;
for i = 1:nProfiles
    uniqueSwitches = unique(availabilitySwitchOverall(i,:));
    numSwitchingValues = max(numSwitchingValues,numel(uniqueSwitches));
    switchingValues(i,1:numel(uniqueSwitches)) = uniqueSwitches;
end

switchingValues = switchingValues(:, 1:numSwitchingValues);
uniqueSwitchingValues = unique(switchingValues);
disp(uniqueSwitchingValues);
 
disp([[1:42]', switchingValues(:,2)])


unavailable1(generalScoreWin1 == 0) = NaN;
unavailable2(generalScoreWin2 == 0) = NaN;
unavailable3(generalScoreWin3 == 0) = NaN;


generalScoreWinScale1 = s .* generalScoreWin1 .* unavailable1;
generalScoreWinScale2 = s .* generalScoreWin2 .* unavailable2;
generalScoreWinScale3 = s .* generalScoreWin3 .* unavailable3;

minScore1 = min(generalScoreWinScale1, [], 2);
minScore2 = min(generalScoreWinScale2, [], 2);
minScore3 = min(generalScoreWinScale3, [], 2);

maxScore1 = max(generalScoreWinScale1, [], 2);
maxScore2 = max(generalScoreWinScale2, [], 2);
maxScore3 = max(generalScoreWinScale3, [], 2);

rangeScore1 = [minScore1, maxScore1];
rangeScore2 = [minScore2, maxScore2];
rangeScore3 = [minScore3, maxScore3];



%% %%%%%%%%%%%%%%%%     Condorcet extensions           %%%%%%%%%%%%%%%%%%%%


%% Get Condorcet outcome

% Get binary comparisons


bin1beats2 = find(pref == 1) < find(pref == 2);
bin1beats2 = reshape(bin1beats2, nRankings, nCitizens);
bin2beats1 = 1 - bin1beats2;

bin1beats3 = find(pref == 1) < find(pref == 3);
bin1beats3 = reshape(bin1beats3, nRankings, nCitizens);
bin3beats1 = 1 - bin1beats3;

bin2beats3 = find(pref == 2) < find(pref == 3);
bin2beats3 = reshape(bin2beats3, nRankings, nCitizens);
bin3beats2 = 1 - bin2beats3;

% Get winner in each binary comparison

vote1beats2 = sum(bin1beats2, 2) >= nCitizens/2;
vote2beats1 = sum(bin2beats1, 2) >= nCitizens/2;
vote1beats3 = sum(bin1beats3, 2) >= nCitizens/2;
vote3beats1 = sum(bin3beats1, 2) >= nCitizens/2;
vote2beats3 = sum(bin2beats3, 2) >= nCitizens/2;
vote3beats2 = sum(bin3beats2, 2) >= nCitizens/2;


condorcetWin1 = vote1beats2 & vote1beats3; 
condorcetWin2 = vote2beats1 & vote2beats3; 
condorcetWin3 = vote3beats1 & vote3beats2; 

winSet = [condorcetWin1 condorcetWin2 condorcetWin3 ];
condorcetWinNum = sum(winSet, 2);

condorcetWin1 = (condorcetWinNum == 0) | (condorcetWinNum == 1 & condorcetWin1 == 1);
condorcetWin2 = (condorcetWinNum == 0) | (condorcetWinNum == 1 & condorcetWin2 == 1);
condorcetWin3 = (condorcetWinNum == 0) | (condorcetWinNum == 1 & condorcetWin3 == 1);

condorcetTie = condorcetWinNum > 1; 


%% Get unanimity (Pareto) outcome

unanimity1beats2 = sum(bin1beats2, 2) == 5;
unanimity2beats1 = sum(bin2beats1, 2) == 5;
unanimity1beats3 = sum(bin1beats3, 2) == 5;
unanimity3beats1 = sum(bin3beats1, 2) == 5;
unanimity2beats3 = sum(bin2beats3, 2) == 5;
unanimity3beats2 = sum(bin3beats2, 2) == 5;

unanimityWin1 = ~unanimity3beats1 & ~unanimity2beats1;
unanimityWin2 = ~unanimity3beats2 & ~unanimity1beats2;
unanimityWin3 = ~unanimity1beats3 & ~unanimity2beats3;


%% %%%%%%%%%%%%%%%%     Runoff methods           %%%%%%%%%%%%%%%%%%%%

%% Plurality with runoff

pluralityRanks = zeros(size(pluralityCount));
for i = 1 : nRankings
    [~,pluralityRanks(:,i)]=ismember(pluralityCount(:,i),sort(pluralityCount(:,i),'descend'));
end

% Keep options ranked first, drop all other options
keepOption = pluralityRanks == 1; 
keepOption = repmat(keepOption, 1, 1, 5);
dropOption = pluralityRanks ~= 1; 
dropOption = repmat(dropOption, 1, 1, 5);

prefKept = pref.*keepOption + (1+nOptions) * dropOption; 

prefKeptRank = zeros(size(prefKept));
for i = 1 : nRankings
    for j = 1 : nCitizens
        [~,prefKeptRank(:,i, j)]=ismember(prefKept(:,i, j),sort(prefKept(:,i, j),'ascend'));
    end
end

pointsScoreRunoff = sum(prefKeptRank == 1, 3); 

[maxPointsScoreRunoff, runoffWinner] = max(pointsScoreRunoff, [], 1);


runoffWin1 = maxPointsScoreRunoff == pointsScoreRunoff(1, :);
runoffWin2 = maxPointsScoreRunoff == pointsScoreRunoff(2, :);
runoffWin3 = maxPointsScoreRunoff == pointsScoreRunoff(3, :);

uniquePlurality = ~tiePlurality;
runoffWinner(uniquePlurality) = pluralityWinner(uniquePlurality);

runoffWin1(uniquePlurality) = pluralityWin1(uniquePlurality);
runoffWin2(uniquePlurality) = pluralityWin2(uniquePlurality);
runoffWin3(uniquePlurality) = pluralityWin3(uniquePlurality);


runoffWinnerAlt = 1*runoffWin1 + 2*runoffWin2 + 3*runoffWin3; 
runoffTest = runoffWinner == runoffWinnerAlt';


%% Scoring runoff NEW


% In the variable prefRankInCell, row 1 represents option A, row 2
% represents option B, row 3 represents option C. The entries in the cells
% (integers) represent the preference rank, with lower values meaning
% stronger preference. 

[prefRankInCell(:, 1), ~] = find(pref == 1);
[prefRankInCell(:, 2), ~] = find(pref == 2);
[prefRankInCell(:, 3), ~] = find(pref == 3);

prefRankInCell = prefRankInCell';
prefRankInCell = reshape(prefRankInCell, nOptions, nRankings, nCitizens);


nScoringRules = 101;
s = linspace(0, 1, nScoringRules);

scoringRunoffWinner = zeros(nOptions, nRankings, nScoringRules);

for sIndex = 1 : nScoringRules
    % Get score

    
    
    scores =    1 * (prefRankInCell == 1) ... 
              + s(sIndex) * (prefRankInCell == 2) ...
              + 0 * (prefRankInCell == 0);

    scoreSum = sum(scores, 3);    
    scoreMin = min(scoreSum, [], 1);
    scoreMin = repmat(scoreMin, [3 1]);

    % Replace rank for option with lowest score across citizens with 10

    dropOption = scoreSum == scoreMin;
    tie1stStage = sum(dropOption, 1) > 1;
    dropOption = repmat(dropOption, [1 1 5]);

    pref2ndStage = prefRankInCell .* (1 - dropOption);
    pref2ndStage = pref2ndStage + 10 * (pref2ndStage == 0);

    % Assign numbers 1, 2, 3 to the preference profile created above

    [pref2ndStageSorted, sortIndex] = sort(pref2ndStage, 1, 'ascend');
    [~, reverseSortIndex] = ismember(pref2ndStage, pref2ndStageSorted); 

    pref2ndStage = pref(sortIndex);
    
    clear prefRankInCell2
    [prefRankInCell2(:, 1), ~] = find(pref2ndStage == 1);
    [prefRankInCell2(:, 2), ~] = find(pref2ndStage == 2);
    [prefRankInCell2(:, 3), ~] = find(pref2ndStage == 3);

    prefRankInCell2 = prefRankInCell2';
    prefRankInCell2 = reshape(prefRankInCell2, nOptions, nRankings, nCitizens);

    scores2ndStage  =    1 * (prefRankInCell2 == 1) ... 
                       + s(sIndex) * (prefRankInCell2 == 2) ...
                       + 0 * (prefRankInCell2 == 0);
                   
    %OLD starts here
    %pref2ndStageSortedNew = ones(size(pref2ndStageSorted));

    %pref2ndStageSortedNew(1, :, :) = 1 * pref2ndStageSortedNew(1, :, :);
    %pref2ndStageSortedNew(2, :, :) = 2 * pref2ndStageSortedNew(2, :, :);
    %pref2ndStageSortedNew(3, :, :) = 3 * pref2ndStageSortedNew(3, :, :);

    %pref2ndStageNew = pref2ndStageSortedNew(reverseSortIndex);

    % Run the scoring rule again

    %scores2ndStage  =    1 * (pref2ndStageNew == 1) ... 
    %                   + s(sIndex) * (pref2ndStageNew == 2) ...
    %                   + 0 * (pref2ndStageNew == 0);
   
                   
    scoreSum2ndStage = sum(scores2ndStage, 3);    
    scoreSum2ndStageMax = max(scoreSum2ndStage, [], 1);
    
    tie2ndStage = sum(scoreSum2ndStageMax == scoreSum2ndStage, 1) > 1;

    scoringRunoffWinner(:,:,sIndex) = scoreSum2ndStageMax == scoreSum2ndStage;

    % Profile 20 is the only one in which some scoring rule (those with
    % parameter weakly exceeding 1/2 are completely irresolute. For those 
    % cases, no options can be removed. Design tie across all options as
    % result. 
    
    for i = 51:101
        scoringRunoffWinner(:,20,i) = [1 1 1];
    end
    
    noTie1stStage = 1 - tie1stStage;
    noTie2ndStage = 1 - tie2ndStage;
    noTie = noTie1stStage & noTie2ndStage;
    
    tieNaN = zeros(size(noTie));
    tieNaN(noTie == 0) = 10;
    noTie = tieNaN;
    noTie = repmat(noTie, 3, 1);
    
    % scoringRunoffWinner(:,:,sIndex) = scoringRunoffWinner(:,:,sIndex) + noTie;
    
end


% Number of scoring runoffs that can be distinguished


profilesMainExperiment = [0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0];
profilesPolitical =      zeros(1, nProfiles);
profilesPolitical(3) = 1;
profilesPolitical(4) = 1;
profilesPolitical(9) = 1;
profilesPolitical(10) = 1;
profilesPolitical(15) = 1;
profilesPolitical(30) = 1;
profilesPolitical(31) = 1;

 
profilesGermanPoliticians =      zeros(1, nProfiles);
profilesGermanPoliticians(3) = 1;
profilesGermanPoliticians(4) = 1;
profilesGermanPoliticians(9) = 1;
profilesGermanPoliticians(21) = 1;
profilesGermanPoliticians(30) = 1;
profilesGermanPoliticians(31) = 1;

profilesMainExperiment = logical(profilesMainExperiment);
profilesPolitical = logical(profilesPolitical);
profilesGermanPoliticians = logical(profilesGermanPoliticians);


scoringRunoffWinnerShort = scoringRunoffWinner(:, profilesMainExperiment, :);

scoringRunoffWinnerShortShifted = scoringRunoffWinnerShort(:,:,2:nScoringRules);

scoringRunoffWinnerShortTest = scoringRunoffWinnerShortShifted - scoringRunoffWinnerShort(:,:,1:nScoringRules - 1);


scoringRunoffWinnerShortTest = scoringRunoffWinnerShortTest ~= 0; 
scoringRunoffWinnerShortTest = sum(scoringRunoffWinnerShortTest, 1);
scoringRunoffWinnerShortTest = sum(scoringRunoffWinnerShortTest, 2);
scoringRunoffWinnerShortTest = squeeze(scoringRunoffWinnerShortTest);

numDiff = scoringRunoffWinnerShortTest > 0;
numDiff = sum(numDiff);


indexMapping = 1:42;
indexMapping = indexMapping(profilesMainExperiment);



% Testing for differences other than those created through ties

scoringRunoffWinnerLargeMagnitude = (1 - (abs(scoringRunoffWinnerShortTest) > 5) );
scoringRunoffWinnerShortTest = scoringRunoffWinnerLargeMagnitude.*scoringRunoffWinnerShortTest;





%% %% Other methods %%%%

%% Supermajority


bin1beats2 = find(pref == 1) < find(pref == 2);
bin1beats2 = reshape(bin1beats2, nRankings, nCitizens);
bin2beats1 = 1 - bin1beats2;

bin1beats3 = find(pref == 1) < find(pref == 3);
bin1beats3 = reshape(bin1beats3, nRankings, nCitizens);
bin3beats1 = 1 - bin1beats3;

bin2beats3 = find(pref == 2) < find(pref == 3);
bin2beats3 = reshape(bin2beats3, nRankings, nCitizens);
bin3beats2 = 1 - bin2beats3;

% Get winner in each binary comparison

superm1beats2 = sum(bin1beats2, 2) >= 4;
superm2beats1 = sum(bin2beats1, 2) >= 4;
superm1beats3 = sum(bin1beats3, 2) >= 4;
superm3beats1 = sum(bin3beats1, 2) >= 4;
superm2beats3 = sum(bin2beats3, 2) >= 4;
superm3beats2 = sum(bin3beats2, 2) >= 4;

supermWin1 = superm1beats2 & superm1beats3; 
supermWin2 = superm2beats1 & superm2beats3; 
supermWin3 = superm3beats1 & superm3beats2; 

winSet = [supermWin1 supermWin2 supermWin3 ];
supermWinNum = sum(winSet, 2);

supermWin1 = (supermWinNum == 0) | (supermWinNum == 1 & supermWin1 == 1);
supermWin2 = (supermWinNum == 0) | (supermWinNum == 1 & supermWin2 == 1);
supermWin3 = (supermWinNum == 0) | (supermWinNum == 1 & supermWin3 == 1);



%% %%%%%%%%%%%%%%%%     Non-benevolent rules           %%%%%%%%%%%%%%%%%%%%

%% Malevolent supermajority

bin1beats2 = find(pref == 1) > find(pref == 2);
bin1beats2 = reshape(bin1beats2, nRankings, nCitizens);
bin2beats1 = 1 - bin1beats2;

bin1beats3 = find(pref == 1) > find(pref == 3);
bin1beats3 = reshape(bin1beats3, nRankings, nCitizens);
bin3beats1 = 1 - bin1beats3;

bin2beats3 = find(pref == 2) > find(pref == 3);
bin2beats3 = reshape(bin2beats3, nRankings, nCitizens);
bin3beats2 = 1 - bin2beats3;

% Get winner in each binary comparison

badSuperm1beats2 = sum(bin1beats2, 2) >= 4;
badSuperm2beats1 = sum(bin2beats1, 2) >= 4;
badSuperm1beats3 = sum(bin1beats3, 2) >= 4;
badSuperm3beats1 = sum(bin3beats1, 2) >= 4;
badSuperm2beats3 = sum(bin2beats3, 2) >= 4;
badSuperm3beats2 = sum(bin3beats2, 2) >= 4;

badSupermWin1 = badSuperm1beats2 & badSuperm1beats3; 
badSupermWin2 = badSuperm2beats1 & badSuperm2beats3; 
badSupermWin3 = badSuperm3beats1 & badSuperm3beats2; 

winSet = [badSupermWin1 badSupermWin2 badSupermWin3 ];
badSupermWinNum = sum(winSet, 2);

badSupermWin1 = (badSupermWinNum == 0) | (badSupermWinNum == 1 & badSupermWin1 == 1);
badSupermWin2 = (badSupermWinNum == 0) | (badSupermWinNum == 1 & badSupermWin2 == 1);
badSupermWin3 = (badSupermWinNum == 0) | (badSupermWinNum == 1 & badSupermWin3 == 1);


%% Get malevolent scoring outcome

nScoringRules = 101;
s = linspace(0, 1, nScoringRules);


malevolentScore = zeros(3, nRankings, nCitizens, nScoringRules);

malevolentScoreWinner = zeros(nRankings, nScoringRules);
malevolentScoreWin1 = zeros(nRankings, nScoringRules);
malevolentScoreWin2 = zeros(nRankings, nScoringRules);
malevolentScoreWin3 = zeros(nRankings, nScoringRules);


for sNr = 1:nScoringRules 
    weights = [0, s(sNr), 1];
    malevolentScore(:,:,:,sNr) =   weights(1) * (pointsScore == 1) ...
                   + weights(2) * (pointsScore == 2) ...
                   + weights(3) * (pointsScore == 3) ;


    malevolentScoreCount = sum(malevolentScore(:,:,:,sNr), 3);
    [maxmalevolentScore, malevolentScoreWinnerS] = max(malevolentScoreCount, [], 1);

    malevolentScoreWin1(:, sNr) = malevolentScoreCount(1,:) == maxmalevolentScore;
    malevolentScoreWin2(:, sNr) = malevolentScoreCount(2,:) == maxmalevolentScore;
    malevolentScoreWin3(:, sNr) = malevolentScoreCount(3,:) == maxmalevolentScore;

    winSetmalevolentScore = repmat(maxmalevolentScore, nOptions, 1) == malevolentScoreCount;
    winNummalevolentScore = sum(winSetmalevolentScore, 1);

    tiemalevolentScore = winNummalevolentScore > 1; 

    malevolentScoreWinner(:,sNr) = max(malevolentScoreWinnerS - nOptions * tiemalevolentScore, -1);
end


%% Malevolent scoring runoff

if 1
    generalbadScore = squeeze(sum(malevolentScore, 3)); 
    badScoreRanks = zeros(size(generalbadScore));
    for i = 1 : nRankings
        for s = 1 : nScoringRules
            [~,badScoreRanks(:,i,s)]=ismember(generalbadScore(:,i,s),sort(generalbadScore(:,i,s),'descend'));
        end
    end


    % Drop options ranked third, keep remaining other options
    keepOption = badScoreRanks ~= 3; 
    keepOption = reshape(keepOption, nOptions, nProfiles, 1, nScoringRules);
    keepOption = repmat(keepOption, 1, 1, 5, 1);
    dropOption = badScoreRanks == 3; 
    dropOption = reshape(dropOption, nOptions, nProfiles, 1, nScoringRules);
    dropOption = repmat(dropOption, 1, 1, 5, 1);

    % New preference profile, in which the preference for the dropped option is
    % set to 4: 
    prefKept = pref.*keepOption + (1+nOptions) * dropOption; 

    % Relabel preference ranks to 1 and 2, such that the unavailable option
    % receives preference 3 for all Choosers. 
    for i = 1 : nRankings
        for s = 1 : nScoringRules
            for j = 1 : 5 
                [~,prefKept(:,i,j,s)]=ismember(prefKept(:,i,j,s),sort(prefKept(:,i,j,s),'ascend'));
            end
        end
    end

    % Find winner among remaining two options
    % (ties occur in prefKept only if first stage is resolute)

    badScoreSecondStage = prefKept == 1;
    badScoreSecondStage = squeeze(sum(badScoreSecondStage, 3));
    maxbadScoreSecondStage = max(badScoreSecondStage, [], 1);
    maxbadScoreSecondStage = repmat(maxbadScoreSecondStage, [3 1 1]);
    badScoreWinnerSecondStage = badScoreSecondStage == maxbadScoreSecondStage;


    % Declare winner as first or second stage winner

    badScoringRunoffWin1 = squeeze(badScoreWinnerSecondStage(1,:,:));
    badScoringRunoffWin2 = squeeze(badScoreWinnerSecondStage(2,:,:));
    badScoringRunoffWin3 = squeeze(badScoreWinnerSecondStage(3,:,:));
    
end



%% Get malevolent Condorcet outcome

% Get binary comparisons

bin1beats2 = find(pref == 1) < find(pref == 2);
bin1beats2 = reshape(bin1beats2, nRankings, nCitizens);
bin2beats1 = 1 - bin1beats2;

bin1beats3 = find(pref == 1) < find(pref == 3);
bin1beats3 = reshape(bin1beats3, nRankings, nCitizens);
bin3beats1 = 1 - bin1beats3;

bin2beats3 = find(pref == 2) < find(pref == 3);
bin2beats3 = reshape(bin2beats3, nRankings, nCitizens);
bin3beats2 = 1 - bin2beats3;

% Get loser in each binary comparison

vote1beats2 = sum(bin1beats2, 2) <= nCitizens/2;
vote2beats1 = sum(bin2beats1, 2) <= nCitizens/2;
vote1beats3 = sum(bin1beats3, 2) <= nCitizens/2;
vote3beats1 = sum(bin3beats1, 2) <= nCitizens/2;
vote2beats3 = sum(bin2beats3, 2) <= nCitizens/2;
vote3beats2 = sum(bin3beats2, 2) <= nCitizens/2;

badCondorcetWins1 = vote1beats2 + vote1beats3;
badCondorcetWins2 = vote2beats1 + vote2beats3;
badCondorcetWins3 = vote3beats1 + vote3beats2;

badCondorcetWins = [badCondorcetWins1 badCondorcetWins2 badCondorcetWins3];

[maxbadCondorcetWins, badCondorcetWinner] = max(badCondorcetWins, [], 2); 


badCondorcetWin1 = maxbadCondorcetWins == badCondorcetWins1;
badCondorcetWin2 = maxbadCondorcetWins == badCondorcetWins2;
badCondorcetWin3 = maxbadCondorcetWins == badCondorcetWins3;


%% Malevolent plurality outcome


badPointsScore = zeros(nRankings * nCitizens, nOptions);

[badPointsScore(:, 1), ~] = find(pref == 1);
[badPointsScore(:, 2), ~] = find(pref == 2);
[badPointsScore(:, 3), ~] = find(pref == 3);

badPointsScore = badPointsScore';
badPointsScore = reshape(badPointsScore, nOptions, nRankings, nCitizens);


weights = [0, 0, 1];
badPluralityScore = weights(1) * (badPointsScore == 1) ...
              	  + weights(2) * (badPointsScore == 2) ...
              	  + weights(3) * (badPointsScore == 3) ;
             
badPluralityCount = sum(badPluralityScore, 3);

[maxbadPluralityCount, badPluralityWinner] = max(badPluralityCount, [], 1);



badPluralityWin1 = maxbadPluralityCount == badPluralityCount(1, :);
badPluralityWin2 = maxbadPluralityCount == badPluralityCount(2, :);
badPluralityWin3 = maxbadPluralityCount == badPluralityCount(3, :);

winSetbadPlurality = repmat(maxbadPluralityCount, nOptions, 1) == badPluralityCount;
winNumbadPlurality = sum(winSetbadPlurality, 1);

tiebadPlurality = winNumbadPlurality > 1; 
badPluralityWinner = max(badPluralityWinner - nOptions * tiebadPlurality, -1);

uniqueBadPlurality = ~tiebadPlurality;

%% Malevolent unanimity (Pareto) outcome

badUnanimity1beats2 = sum(bin1beats2, 2) == 0;
badUnanimity2beats1 = sum(bin2beats1, 2) == 0;
badUnanimity1beats3 = sum(bin1beats3, 2) == 0;
badUnanimity3beats1 = sum(bin3beats1, 2) == 0;
badUnanimity2beats3 = sum(bin2beats3, 2) == 0;
badUnanimity3beats2 = sum(bin3beats2, 2) == 0;

badUnanimityWin1 = ~badUnanimity3beats1 & ~badUnanimity2beats1;
badUnanimityWin2 = ~badUnanimity3beats2 & ~badUnanimity1beats2;
badUnanimityWin3 = ~badUnanimity1beats3 & ~badUnanimity2beats3;


%% Malevolent Runoff
 
badPluralityRanks = zeros(size(badPluralityCount));
for i = 1 : nRankings
    [~,badPluralityRanks(:,i)]=ismember(badPluralityCount(:,i),sort(badPluralityCount(:,i),'ascend'));
end
 
% Keep options ranked first, drop all other options
keepOption = badPluralityRanks == 2; 
keepOption = repmat(keepOption, 1, 1, 5);
dropOption = badPluralityRanks ~= 2; 
dropOption = repmat(dropOption, 1, 1, 5);
 
prefKept = pref.*keepOption + (1+nOptions) * dropOption; 
 
prefKeptRank = zeros(size(prefKept));
for i = 1 : nRankings
    for j = 1 : nCitizens
        [~,prefKeptRank(:,i, j)]=ismember(prefKept(:,i, j),sort(prefKept(:,i, j),'descend'));
    end
end
 
pointsScorebadRunoff = sum(prefKeptRank == 2, 3); 
 
[maxPointsScorebadRunoff, badRunoffWinner] = max(pointsScorebadRunoff, [], 1);
 
badRunoffWin1 = maxPointsScorebadRunoff == pointsScorebadRunoff(1, :);
badRunoffWin2 = maxPointsScorebadRunoff == pointsScorebadRunoff(2, :);
badRunoffWin3 = maxPointsScorebadRunoff == pointsScorebadRunoff(3, :);
 
uniquebadPlurality = ~tiebadPlurality;
badRunoffWinner(uniquebadPlurality) = badPluralityWinner(uniquebadPlurality);
 
badRunoffWin1(uniquebadPlurality) = badPluralityWin1(uniquebadPlurality);
badRunoffWin2(uniquebadPlurality) = badPluralityWin2(uniquebadPlurality);
badRunoffWin3(uniquebadPlurality) = badPluralityWin3(uniquebadPlurality);
 
 
badRunoffWinnerAlt = 1*badRunoffWin1 + 2*badRunoffWin2 + 3*badRunoffWin3; 
badRunoffTest = badRunoffWinner == badRunoffWinnerAlt';
 
 
result = [badRunoffWin1', badRunoffWin2', badRunoffWin3', runoffWin1', runoffWin2', runoffWin3'];


%% Get full randomization rule

randomizationWin1 = ones(size(pluralityWin1));
randomizationWin2 = ones(size(pluralityWin2));
randomizationWin3 = ones(size(pluralityWin3));




%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%   Organize data for export and plotting    %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


bordaPlot = [+bordaWin1; ...
             +bordaWin2; ...
             +bordaWin3];
         

condorcetPlot = [+condorcetWin1'; ...
                 +condorcetWin2'; ...
                 +condorcetWin3'];

unanimityPlot = [+unanimityWin1'; ...
                 +unanimityWin2'; ...
                 +unanimityWin3'];
             
             
runoffPlot = [+runoffWin1; ...
              +runoffWin2; ...
              +runoffWin3];

supermPlot = [+supermWin1'; ...
              +supermWin2'; ...
              +supermWin3'];
 
badSupermPlot = [+badSupermWin1'; ...
                 +badSupermWin2'; ...
                 +badSupermWin3'];
  
badRunoffPlot = [+badRunoffWin1; ...
                 +badRunoffWin2; ...
                 +badRunoffWin3];

          
% Boundary scoring cases         
          
scoring0Plot = [+generalScoreWin1(:, 1 + 100*0)'; ...
                +generalScoreWin2(:, 1 + 100*0)'; ...
                +generalScoreWin3(:, 1 + 100*0)'];          
          
scoring100Plot = [+generalScoreWin1(:, 1 + 100)'; ...
                  +generalScoreWin2(:, 1 + 100)'; ...
                  +generalScoreWin3(:, 1 + 100)'];   
              

malevolentScore0Plot = [+malevolentScoreWin1(:, 1 + 100*0)'; ...
                        +malevolentScoreWin2(:, 1 + 100*0)'; ...
                        +malevolentScoreWin3(:, 1 + 100*0)'];          
 
malevolentScore100Plot = [+malevolentScoreWin1(:, 1 + 100)'; ...
                          +malevolentScoreWin2(:, 1 + 100)'; ...
                          +malevolentScoreWin3(:, 1 + 100)'];             
              
           
                      
                      
% Interior scoring cases     
              

            
scoring0to33Plot = [+generalScoreWin1(:, 1 + 30)'; ...
                    +generalScoreWin2(:, 1 + 30)'; ...
                    +generalScoreWin3(:, 1 + 30)'];  
          
scoring33to50Plot = [+generalScoreWin1(:, 1 + 40)'; ...
                     +generalScoreWin2(:, 1 + 40)'; ...
                     +generalScoreWin3(:, 1 + 40)'];          

scoring50to60Plot = [+generalScoreWin1(:, 1 + 55)'; ...
                     +generalScoreWin2(:, 1 + 55)'; ...
                     +generalScoreWin3(:, 1 + 55)'];          
            
scoring60to67Plot = [+generalScoreWin1(:, 1 + 63)'; ...
                     +generalScoreWin2(:, 1 + 63)'; ...
                     +generalScoreWin3(:, 1 + 63)'];  
                      
scoring67to75Plot = [+generalScoreWin1(:, 1 + 70)'; ...
                     +generalScoreWin2(:, 1 + 70)'; ...
                     +generalScoreWin3(:, 1 + 70)'];          

scoring75to80Plot = [+generalScoreWin1(:, 1 + 78)'; ...
                     +generalScoreWin2(:, 1 + 78)'; ...
                     +generalScoreWin3(:, 1 + 78)'];          

scoring80to100Plot = [+generalScoreWin1(:, 1 + 85)'; ...
                      +generalScoreWin2(:, 1 + 85)'; ...
                      +generalScoreWin3(:, 1 + 85)'];     
                  

scoring33Plot = max(scoring0to33Plot,  scoring33to50Plot);                  
scoring50Plot = max(scoring33to50Plot, scoring50to60Plot);
scoring60Plot = max(scoring50to60Plot, scoring60to67Plot);
scoring67Plot = max(scoring60to67Plot, scoring67to75Plot);
scoring75Plot = max(scoring67to75Plot, scoring75to80Plot);
scoring80Plot = max(scoring75to80Plot, scoring80to100Plot);

                       
% Interior scoring runoff cases     
              

scoringRunoff0to33Plot = scoringRunoffWinner(:,:,1); 
scoringRunoff33to50Plot = scoringRunoffWinner(:,:,40); 
scoringRunoff50to100Plot = scoringRunoffWinner(:,:,60); 
scoringRunoff100Plot = scoringRunoffWinner(:,:,101); 
             
                  
% Malevolent scoring rules                  
                   
malevolentScore0to33Plot = [+malevolentScoreWin1(:, 1 + 30)'; ...
                            +malevolentScoreWin2(:, 1 + 30)'; ...
                            +malevolentScoreWin3(:, 1 + 30)'];  
          
malevolentScore33to50Plot = [+malevolentScoreWin1(:, 1 + 40)'; ...
                             +malevolentScoreWin2(:, 1 + 40)'; ...
                             +malevolentScoreWin3(:, 1 + 40)'];          

malevolentScore50to60Plot = [+malevolentScoreWin1(:, 1 + 55)'; ...
                             +malevolentScoreWin2(:, 1 + 55)'; ...
                             +malevolentScoreWin3(:, 1 + 55)'];          
            
malevolentScore60to67Plot = [+malevolentScoreWin1(:, 1 + 63)'; ...
                             +malevolentScoreWin2(:, 1 + 63)'; ...
                             +malevolentScoreWin3(:, 1 + 63)'];  
                      
malevolentScore67to75Plot = [+malevolentScoreWin1(:, 1 + 70)'; ...
                             +malevolentScoreWin2(:, 1 + 70)'; ...
                             +malevolentScoreWin3(:, 1 + 70)'];          

malevolentScore75to80Plot = [+malevolentScoreWin1(:, 1 + 78)'; ...
                             +malevolentScoreWin2(:, 1 + 78)'; ...
                             +malevolentScoreWin3(:, 1 + 78)'];          
        
malevolentScore80to100Plot = [+malevolentScoreWin1(:, 1 + 85)'; ...
                              +malevolentScoreWin2(:, 1 + 85)'; ...
                              +malevolentScoreWin3(:, 1 + 85)'];          


malevolentScore33Plot = max(malevolentScore0to33Plot,  malevolentScore33to50Plot);                  
malevolentScore50Plot = max(malevolentScore33to50Plot, malevolentScore50to60Plot);
malevolentScore60Plot = max(malevolentScore50to60Plot, malevolentScore60to67Plot);
malevolentScore67Plot = max(malevolentScore60to67Plot, malevolentScore67to75Plot);
malevolentScore75Plot = max(malevolentScore67to75Plot, malevolentScore75to80Plot);
malevolentScore80Plot = max(malevolentScore75to80Plot, malevolentScore80to100Plot);                  
         


% Malevolent scoring runoffs

 
badScoringRunoff0Plot = [+badScoringRunoffWin1(:, 1)'; ...
                      +badScoringRunoffWin2(:, 1)'; ...
                      +badScoringRunoffWin3(:, 1)'];  
                                
badScoringRunoff0to33Plot = [+badScoringRunoffWin1(:, 1 + 30)'; ...
                    +badScoringRunoffWin2(:, 1 + 30)'; ...
                    +badScoringRunoffWin3(:, 1 + 30)'];  
          
badScoringRunoff33to50Plot = [+badScoringRunoffWin1(:, 1 + 40)'; ...
                     +badScoringRunoffWin2(:, 1 + 40)'; ...
                     +badScoringRunoffWin3(:, 1 + 40)'];          
 
badScoringRunoff50to60Plot = [+badScoringRunoffWin1(:, 1 + 55)'; ...
                     +badScoringRunoffWin2(:, 1 + 55)'; ...
                     +badScoringRunoffWin3(:, 1 + 55)'];          
            
badScoringRunoff60to67Plot = [+badScoringRunoffWin1(:, 1 + 63)'; ...
                     +badScoringRunoffWin2(:, 1 + 63)'; ...
                     +badScoringRunoffWin3(:, 1 + 63)'];  
                      
badScoringRunoff67to75Plot = [+badScoringRunoffWin1(:, 1 + 70)'; ...
                     +badScoringRunoffWin2(:, 1 + 70)'; ...
                     +badScoringRunoffWin3(:, 1 + 70)'];          
 
badScoringRunoff75to80Plot = [+badScoringRunoffWin1(:, 1 + 78)'; ...
                     +badScoringRunoffWin2(:, 1 + 78)'; ...
                     +badScoringRunoffWin3(:, 1 + 78)'];          
 
badScoringRunoff80to100Plot = [+badScoringRunoffWin1(:, 1 + 85)'; ...
                      +badScoringRunoffWin2(:, 1 + 85)'; ...
                      +badScoringRunoffWin3(:, 1 + 85)'];     
               
badScoringRunoff100Plot = [+badScoringRunoffWin1(:, 101)'; ...
                      +badScoringRunoffWin2(:, 101)'; ...
                      +badScoringRunoffWin3(:, 101)'];  

% Other malevolent rules                 
                  

badCondorcetPlot = [+badCondorcetWin1'; ...
                    +badCondorcetWin2'; ...
                    +badCondorcetWin3'];
             
randomizationPlot = [+randomizationWin1; ...
                     +randomizationWin2; ...
                     +randomizationWin3];          
              
badUnanimityPlot = [+badUnanimityWin1'; ...
                    +badUnanimityWin2'; ...
                    +badUnanimityWin3'];





%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%    Export data                  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

profileIndex = 1:42;
profileIndex = profileIndex';


ruleResults = cat(3, ...
			scoring0Plot, ...
			scoring0to33Plot, ...
            scoring33Plot, ...
			scoring33to50Plot, ...
            scoring50Plot, ...
			scoring50to60Plot, ...
            scoring60Plot, ...
			scoring60to67Plot, ...
            scoring67Plot, ...
			scoring67to75Plot, ...
            scoring75Plot, ...
			scoring75to80Plot, ...
            scoring80Plot, ...
			scoring80to100Plot, ...
			scoring100Plot, ...
            scoringRunoff0to33Plot, ...    
            scoringRunoff33to50Plot, ...
            scoringRunoff50to100Plot, ...
            scoringRunoff100Plot, ...            
            condorcetPlot, ...
            supermPlot, ...
            unanimityPlot, ...
			malevolentScore0Plot, ...
			malevolentScore0to33Plot, ...
            malevolentScore33Plot, ...
			malevolentScore33to50Plot, ...
            malevolentScore50Plot, ...
			malevolentScore50to60Plot, ...
            malevolentScore60Plot, ...
			malevolentScore60to67Plot, ...
            malevolentScore67Plot, ...
			malevolentScore67to75Plot, ...
            malevolentScore75Plot, ...
			malevolentScore75to80Plot, ...
            malevolentScore80Plot, ...
			malevolentScore80to100Plot, ...
			malevolentScore100Plot, ...
            badScoringRunoff0to33Plot, ...
            badScoringRunoff33to50Plot, ...
            badScoringRunoff50to60Plot, ...
            badScoringRunoff100Plot, ...            
			badCondorcetPlot, ...
            badSupermPlot, ...
            badUnanimityPlot);      
        
        
profilesMainExperiment = [0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0];
profilesMainExperiment = logical(profilesMainExperiment);

% Testing number of rules without Pareto violations
prod(2.^sum(unanimityPlot(:,profilesMainExperiment))-1)
prod(2.^sum(unanimityPlot)-1)





ruleNames= [ ...
"scoring"; ...      
"scoring"; ...      
"scoring"; ...      
"scoring"; ...      
"scoring"; ...      
"scoring"; ...      
"scoring"; ...      
"scoring"; ...      
"scoring"; ...      
"scoring"; ...      
"scoring"; ...      
"scoring"; ...      
"scoring"; ...      
"scoring"; ...      
"scoring"; ...      
"scoringRunoff"; ...      
"scoringRunoff"; ...      
"scoringRunoff"; ...      
"scoringRunoff"; ...     
"Condorcet"; ...    
"Superm"; ...
"Unanimity"; ...
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoring"; ...  
"bad scoringRunoff"; ...      
"bad scoringRunoff"; ...      
"bad scoringRunoff"; ...      
"bad scoringRunoff"; ...     
"bad Condorcet"; ...
"bad Superm"; ...   
"bad unanimity"]; 


   
rulePara = ...
           [01, -1, 0; ...
            02, 0, 33; ...
	        03, 33, 33; ...
            04, 33, 50; ...
	        05, 50, 50; ...
            06, 50, 60; ...
	        07, 60, 60; ...
            08, 60, 67; ...
	        09, 67, 67; ...
            10, 67, 75; ...
	        11, 75, 75; ...
            12, 75, 80; ...
	        13, 80, 80; ...
            14, 80, 100; ...
            15, 100, 101; ...
            16, -1, 33; ...
            17, 33, 50; ...
            18, 50, 100; ...
            19, 100, 101; ...
            20, 0, 0; ...
            21, 0, 0; ...
            22, 0, 0; ...
            23, -1, 0; ...
            24, 0, 33; ...
	        25, 33, 33; ...
            26, 33, 50; ...
	        27, 50, 50; ...
            28, 50, 60; ...
	        29, 60, 60; ...
            30, 60, 67; ...
	        31, 67, 67; ...
            32, 67, 75; ...
	        33, 75, 75; ...
            34, 75, 80; ...
	        35, 80, 80; ...
            36, 80, 100; ...
            37, 100, 101; ...
            38, -1, 33; ...
            39, 33, 50; ...
            40, 50, 100; ...
            41, 100, 101; ...
            42, 0, 0; ...
            43, 0, 0; ...
            44, 0, 0];        
     

inclusionExtendedRules = [1 1 0 1 0 1 0 1 0 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 0 1 0 1 0 1 0 1 1 1 1 1 1 1 1];    
inclusionExtendedRules = logical(inclusionExtendedRules);

ruleResultsExtended = ruleResults(:,:,inclusionExtendedRules);
ruleNamesExtended = ruleNames(inclusionExtendedRules);
ruleParaExtended = rulePara(inclusionExtendedRules,:);

rulesInteriorScoring = [1 1 0 1 0 1 0 1 0 1 0 1 0 1 1 0 0 0 0 0 0 0 1 1 0 1 0 1 0 1 0 1 0 1 0 1 1 0 0 0 0 0 0];
rulesInteriorScoring = logical(rulesInteriorScoring);

ruleResultsIntSc = ruleResults(:,:,rulesInteriorScoring);
ruleNamesIntSc = ruleNames(rulesInteriorScoring);
ruleParaIntSc = rulePara(rulesInteriorScoring,:);

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


%% Store the rule predictions data

        
rulePara = repmat(rulePara, 1, 1, nProfiles);
rulePara = permute(rulePara, [2 3 1]);
        
ruleNames = repmat(ruleNames, 1, 1, nProfiles);
ruleNames = permute(ruleNames, [2 3 1]);
        
ruleParaIntSc = repmat(ruleParaIntSc, 1, 1, nProfiles);
ruleParaIntSc = permute(ruleParaIntSc, [2 3 1]);
        
ruleNamesIntSc = repmat(ruleNamesIntSc, 1, 1, nProfiles);
ruleNamesIntSc = permute(ruleNamesIntSc, [2 3 1]);

  
        
ruleData = cat(1, ruleResults, ruleNames, rulePara); 

prefExport = [pref(:,:,1)', pref(:,:,2)', pref(:,:,3)', pref(:,:,4)', pref(:,:,5)'];

ruleDataLong = [profileIndex, ones(nRankings, 1), ruleData(:,:,1)', prefExport];
for r = 2 : nRules 
    ruleDataLong = cat(1, ruleDataLong, [profileIndex, r*ones(nRankings, 1), ruleData(:,:,r)', prefExport]);
end


rowNames = 1:nRankings;
colNames = {'profileIndex', 'ruleIndex', 'optionA', 'optionB', 'optionC', 'ruleName', 'ruleIndex2',  'scLow', 'scHigh', '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(ruleDataLong,'VariableNames',colNames);
writetable(table, 'data/rulePredictions3optNew.csv')



%% Plot rule differences matrix

% Work
ruleResultsExp = ruleResults(:, profilesMainExperiment, :);
ruleResultsExpIndex = 100*ruleResultsExp(1,:,:) + 10*ruleResultsExp(2,:,:) + ruleResultsExp(3,:,:);
ruleResultsExpIndex = squeeze(ruleResultsExpIndex);

[~, numRules] = size(ruleResultsExpIndex);

diffMatrixSize = numRules/2;
diffMatrix = zeros(diffMatrixSize);
for i = 1:diffMatrixSize
    for j = 1:diffMatrixSize
        diffMatrix(i,j) = sum(ruleResultsExpIndex(:, i) ~= ruleResultsExpIndex(:,j));
    end
end

heatmap(diffMatrix)
xlabel('Rule number')
ylabel('Rule number')
title('Number of profiles on which rules differ')
saveas(gcf,'../writeup/graphs/rulesDistanceMatrix.png')

% Politics
ruleResultsExp = ruleResults(:, profilesPolitical, :);
ruleResultsExpIndex = 100*ruleResultsExp(1,:,:) + 10*ruleResultsExp(2,:,:) + ruleResultsExp(3,:,:);
ruleResultsExpIndex = squeeze(ruleResultsExpIndex);

[~, numRules] = size(ruleResultsExpIndex);

diffMatrixSize = numRules/2;
diffMatrix = zeros(diffMatrixSize);
for i = 1:diffMatrixSize
    for j = 1:diffMatrixSize
        diffMatrix(i,j) = sum(ruleResultsExpIndex(:, i) ~= ruleResultsExpIndex(:,j));
    end
end

heatmap(diffMatrix)
xlabel('Rule number')
ylabel('Rule number')
title('Number of profiles on which rules differ')

saveas(gcf,'../writeup/graphs/rulesDistanceMatrixPolitics.png')



% Omitted profiles
ruleResultsExp = ruleResults(:, ~profilesMainExperiment, :);
ruleResultsExpIndex = 100*ruleResultsExp(1,:,:) + 10*ruleResultsExp(2,:,:) + ruleResultsExp(3,:,:);
ruleResultsExpIndex = squeeze(ruleResultsExpIndex);

[~, numRules] = size(ruleResultsExpIndex);

diffMatrixSize = numRules/2;
diffMatrix = zeros(diffMatrixSize);
for i = 1:diffMatrixSize
    for j = 1:diffMatrixSize
        diffMatrix(i,j) = sum(ruleResultsExpIndex(:, i) ~= ruleResultsExpIndex(:,j));
    end
end

heatmap(diffMatrix)
xlabel('Rule number')
ylabel('Rule number')
title('Number of profiles on which rules differ')
saveas(gcf,'../writeup/graphs/rulesDistanceMatrixOmittedProfiles.png')







%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%    Decompose rules into resolute parts             %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Do not use knife-edge scoring rules, as they cannot be distinguished

rule1 = scoring0Plot(:,profilesMainExperiment);
rule2 = scoring100Plot(:,profilesMainExperiment);

ruleSet = cat(3,rule1, rule2);




% Write function that splits a rule at a given index j

% Assume the rule is already resolute before the splitIndex

[X, ruleLength, numExpandRules] = size(ruleSet);
r = 1;
while r < numExpandRules 
    
    rule = ruleSet(:,:,r);
    
    multiplicity = sum(rule, 1);
    
    % Advance to the next rule if the current rule is exhausted
    if max(multiplicity) == 1
        r = r+1;
    else
        
        % Set i to 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)];
            end
            if numel(subrule) == 0 && numel(subrule_j) ~= 0
                subrule = subrule_j;
            elseif numel(subrule_j) ~= 0
                subrule = cat(3, subrule, subrule_j);
            end
        end

        % Replace rule under consideration by expanded rule
        ruleSet = cat(3, ruleSet(:,:,1:r-1), subrule, ruleSet(:,:,r+1:numExpandRules));
    end
    
end






%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%    Plot rules                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% Define sets etc. of profiles to display

all = true(nRankings, 1);

mainExperiment = [0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0];
mainExperiment = logical(mainExperiment);

% Choose set

subset = mainExperiment;


prefPlot = profileCoordinates';          
       
          
                    
xAxis = 1:1:sum(subset);


sortingMatrix = permute(ruleResults, [2 1 3]);             

%[~, permutation] = sortrows(sortingMatrix(:,:), 'descend'); 
permutation = 1:42;

subset              = subset(:, permutation);

scoring0Plot		= scoring0Plot(:, permutation);
scoring0to33Plot		= scoring0to33Plot(:, permutation);
scoring33Plot		= scoring33Plot(:, permutation);
scoring33to50Plot		= scoring33to50Plot(:, permutation);
scoring50Plot		= scoring50Plot(:, permutation);
scoring50to60Plot		= scoring50to60Plot(:, permutation);
scoring60Plot		= scoring60Plot(:, permutation);
scoring60to67Plot		= scoring60to67Plot(:, permutation);
scoring67Plot		= scoring67Plot(:, permutation);
scoring67to75Plot		= scoring67to75Plot(:, permutation);
scoring75Plot		= scoring75Plot(:, permutation);
scoring75to80Plot		= scoring75to80Plot(:, permutation);
scoring80Plot		= scoring80Plot(:, permutation);
scoring80to100Plot		= scoring80to100Plot(:, permutation);
scoring100Plot		= scoring100Plot(:, permutation);
condorcetPlot		= condorcetPlot(:, permutation);
runoffPlot          = runoffPlot(:, permutation);
supermPlot          = supermPlot(:, permutation);
malevolentScore0Plot		= malevolentScore0Plot(:, permutation);
malevolentScore0to33Plot		= malevolentScore0to33Plot(:, permutation);
malevolentScore33Plot		= malevolentScore33Plot(:, permutation);
malevolentScore33to50Plot		= malevolentScore33to50Plot(:, permutation);
malevolentScore50Plot		= malevolentScore50Plot(:, permutation);		
malevolentScore50to60Plot		= malevolentScore50to60Plot(:, permutation);
malevolentScore60Plot		= malevolentScore60Plot(:, permutation);
malevolentScore60to67Plot		= malevolentScore60to67Plot(:, permutation);
malevolentScore67Plot		= malevolentScore67Plot(:, permutation);
malevolentScore67to75Plot		= malevolentScore67to75Plot(:, permutation);
malevolentScore75Plot		= malevolentScore75Plot(:, permutation);
malevolentScore75to80Plot		= malevolentScore75to80Plot(:, permutation);
malevolentScore80Plot		= malevolentScore80Plot(:, permutation);
malevolentScore80to100Plot		= malevolentScore80to100Plot(:, permutation);
malevolentScore100Plot		= malevolentScore100Plot(:, permutation);
badCondorcetPlot		= badCondorcetPlot(:, permutation);
badRunoffPlot          = badRunoffPlot(:, permutation);
badSupermPlot          = badSupermPlot(:, permutation);
randomizationPlot			= randomizationPlot	(:, permutation);		
	 

% Selected rules
if 1 
    
    permutationPaper = [3 21 9 15 32 4 23 31 10 16 30 19 20 24 26 25 35];
    
    clf 
    nRules = 6;
    sPlot = 0;

    sPlot = sPlot + 1;
    subplot(nRules,1,sPlot)
    heatmap(xAxis,{'A>B>C','A>C>B','C>A>B', 'C>B>A', 'B>C>A', 'B>A>C'}, prefPlot(:,permutationPaper) ...
        , 'ColorbarVisible','off' ...
        ); 

    title("Preference profile")

    sPlot = sPlot + 1;
    subplot(nRules,1,sPlot)
    heatmap(xAxis,{'A','B','C'},condorcetPlot(:,permutationPaper));  
    title("Condorcet")

    sPlot = sPlot + 1;
    subplot(nRules,1,sPlot)
    heatmap(xAxis,{'A','B','C'},scoring0Plot(:,permutationPaper));  
    title("Plurality")

    sPlot = sPlot + 1;
    subplot(nRules,1,sPlot)
    heatmap(xAxis,{'A','B','C'},scoring50Plot(:,permutationPaper));  
    title("Borda")

    sPlot = sPlot + 1;
    subplot(nRules,1,sPlot)
    heatmap(xAxis,{'A','B','C'},scoring100Plot(:,permutationPaper));  
    title("Antiplurality")

    sPlot = sPlot + 1;
    subplot(nRules,1,sPlot)
    heatmap(xAxis,{'A','B','C'}, runoffPlot(:,permutationPaper));             
    title("Plurality runoff")

    saveas(gcf,'../writeup/graphs/selectedRules3opt.png')

end



