%% This routine compares a number of models of Krumhansl's probe tone data.
% It calls the function Probe_tone_spectral_sim_model.m

% Author: Andrew J. Milne
% Version: 1.1

clear all
% tic
global spec_sim weight_non_roots weight_m3 K_data_train valid_ind ...
    num_harmonics virtual len_shifts exp_or_sal

%% Options - for binary choices, no = 0 and yes = 1
lambda = 1; % Box-Cox transformation (1 gives no transformation)
cross_validate = 0; % To cross-validate or not? That is the question
num_partitions = 20; % Number of different partitionings (i.e., runs) of the
% 12-fold cross validation

% Choose the model
% model = 'Null'
% model = 'DumbTriad'
% model = 'KrumPrev'
% model = 'KrumCons'
% model = 'Smith'
% model = 'Ler'
% model = 'Butler'
% model = 'P89'
% model = 'P11a'
% model = 'P11b'
% model = 'P94'
 model = 'Milne'
if isequal(model,'Milne')
    weight_non_roots = 1 % Weight non-roots lower than other roots?
    weight_m3 = 1 % Weight minor third of minor triad separately
    
    num_harmonics = 12 % a value of 9 seems to generalize slightly better than the chosen value of 12 (I'll use 12 because it corresponds to the other paper)
    virtual = 0 % Use virtual rather than spectral pitches?
    exp_or_sal = 'exp' % use expectation or salience vectors
end

%% Data to be modelled
% data from Krumhansl (1990)
y_Maj_all = [6.35 2.23 3.48 2.33 4.38 4.09 2.52 5.19 2.39 3.66 2.29 2.88]'; % p.30 - these are the major hierarchies
y_min_all = [6.33 2.68 3.52 5.38 2.60 3.53 2.54 4.75 3.98 2.69 3.34 3.17]'; % p.30 - these are the minor hierarchies

y_Maj_I = [6.20 2.55 3.45 2.85 4.22 4.57 2.67 5.25 2.45 3.35 2.70 2.70]'; % from Parncutt 94
y_min_I = [6.03 3.35 3.67 5.28 2.58 3.55 2.87 4.80 4.25 2.67 3.50 3.52]'; % from Parncutt 94

y_Maj_IV_V_I = [6.28 1.92 3.03 2.47 4.90 5.38 2.55 5.43 2.40 3.25 2.15 2.28]';
y_min_IV_V_I = [6.50 2.50 3.10 5.43 2.55 4.25 2.50 5.20 4.05 2.55 3.22 2.80]';

y_Maj_II_V_I = [6.57 2.30 4.03 2.10 4.12 3.47 2.20 5.28 2.15 3.72 2.28 3.22]';
y_min_II_V_I = [6.70 2.33 4.30 5.43 2.70 3.00 2.33 4.05 3.08 2.70 3.40 3.65]';

y_Maj_VI_V_I = [6.35 2.13 3.40 1.88 4.28 2.95 2.65 4.80 2.58 4.32 2.03 3.33]';
y_min_VI_V_I = [6.10 2.53 3.00 5.40 2.58 3.30 2.47 4.95 4.55 2.83 3.25 2.83]';

% Choose which set of data to model
y_Maj = y_Maj_all;
y_min = y_min_all;

%% More data
x_conKK = [190 285 275 255 250 245 265 215 260 230 250 255]';           % K&K p.57 - these are the Kameoka consonance values (the best measure in the set)
x_conall = [1.00 11.29 9.50 6.60 4.65 3.10 8.28 2.05 6.45 4.85 8.83 10.50; ...
    -11.00 -0.00 -1.50 -4.35 -6.85 -7.00 -3.85 -9.50 -6.15 -8.00 -3.30 -1.50; ...
    0 76 25 24 18 3 18 1 22 22 24 48; ...
    0 70 32 20 8 2 20 0 20 3 23 42; ...
    .0019 .4886 .2690 .1109 .0551 .0451 .0930 .0221 .0843 .0477 .0998 .2312; ...
    190 285 275 255 250 245 265 215 260 230 250 255];           % p.67 - this is the entire set of consonance values
x_stM = [3213 194 3001 352 3111 1947 556 3615 348 1840 361 1504]';      % p.67 - this is the major mode frequencies
x_stm = [906 103 550 564 124 430 117 1042 343 100 259 272]';            % p.67 - this is the minor mode frequencies
x_SPD_M = [0.2373 0.9034 0.6072 0.8414 0.3664 0.7129 0.8501 0.2911 0.7876 0.6605 0.7977 0.7124]'; % these are my spectral pitch distance values calculated with a probe tone and a triad - 10 partials equal amplitude
x_SPD_m = [0.2911 0.8501 0.7129 0.3664 0.8414 0.6072 0.9034 0.2373 0.7124 0.7977 0.6605 0.7876]';

Huron_diss = x_conall([2 5 5],:); % This is the subset of consonance values Huron used in his aggregate dyadic consonance paper
Huron_conmeans = mean(Huron_diss,2);
Huron_constd = std(Huron_diss,1,2);
Huron_conctrd = bsxfun(@minus, Huron_diss, Huron_conmeans);
Huron_constd = bsxfun(@times, Huron_conctrd, 1./Huron_constd);
Huron_conmean = -mean(Huron_constd);

% standardizing data (so mean=0, sdev=1)
% y_Maj = (y_Maj - mean(y_Maj))/std(y_Maj,1);
% y_min = (y_min - mean(y_min))/std(y_min,1);
x_conKK = -(x_conKK - mean(x_conKK))/std(x_conKK,1);
x_stM = (x_stM - mean(x_stM))/std(x_stM,1);
x_stm = (x_stm - mean(x_stm))/std(x_stm,1);
x_SPD_M = -(x_SPD_M - mean(x_SPD_M))/std(x_SPD_M,1);
x_SPD_m = -(x_SPD_m - mean(x_SPD_m))/std(x_SPD_m,1);

% this gets the mean value of the differing standardized consonance ratings
% from p.57, which is what Krumhansl does to get the consonance model for her
% linear regression at the end of chap. 3.
x_conmeans = mean(x_conall,2);
x_constd = std(x_conall,1,2);
x_conctrd = bsxfun(@minus, x_conall, x_conmeans);
x_constd = bsxfun(@times, x_conctrd, 1./x_constd);
x_conmean = mean(x_constd);
x_conall = -x_conmean';

%% Null model
Null_M_m = zeros(24,1);

%% Dumb triad
Triad_M_m = [1 0 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0]';

%% Krumhansl prev
% Top-down
Krum_prev_M = x_stM;
Krum_prev_m = x_stm;
Krum_prev_M_m = [Krum_prev_M; Krum_prev_m]; % Krumhansl's model using the frequencies of scale degrees in corpus

%% Krumhansl con
% Bottom-up
Krum_cons_M = x_conall;
Krum_cons_m = x_conall;
Krum_cons_M_m = [Krum_cons_M; Krum_cons_m]; % Krumhansl's model using the average of the different standardized consonance values

%% Smith
% Bottom-up (building scales stepwise so as to minimize aggregate dissonance)
Smith_M = 12 - [1 11.5 5 11.5 3 6.5 9.5 2 8 4 9.5 6.5]';
Smith_m = 12 - [1 9.5 6.5 3 11.5 5 11.5 2 6.5 9.5 4 8]';
Smith_M_m = [Smith_M; Smith_m];

%% Lerdahl
% Top-down (in that the levels are an intuition of prevalence in music practice)
Lerdahl_M = [5 1 2 1 3 2 1 4 1 2 1 2]';
%Lerdahl_m = [5 1 2 3 1 2 1 4 2 1 1 2]'; % These vals are copied from Parncutt 2011, but it's not clear why the diatonic level is based on the harmonic minor rather than Aeolian mode (the latter gives a higher minor correlation)
Lerdahl_m = [5 1 2 3 1 2 1 4 2 1 2 1]'; % Here Aeolian is the scale level
Lerdahl_M_m = [Lerdahl_M; Lerdahl_m];

% 1  1  1  1  1  1  1  1  1  1  1  1
% 1     1  1     1     1  1     1
% 1        1           1
% 1                    1
% 1
% 0  1  2  3  4  5  6  7  8  9 10 11

%% Butler
Butler_M = [6 0 4 0 5 2 0 7 0 3 0 3]'; % Butler's original model is a simple count of the number of times each scale degree occurs in the Krumhansl context
Butler_m = [6 0 4 5 0 2 0 7 3 0 0 3]';
Butler_M_m = [Butler_M; Butler_m];

% Butler_M = [8 0 4 0 7 2 0 9 0 3 0 3]'; % This is Parncutt's alteration of Butler's data which triple weights the 1 chord p.341 Parncutt 2011
% Butler_m = [8 0 4 7 0 2 0 9 3 0 0 3]'; % this is "acceptable" because Krumhansl's final data averages the four contexts (IV-V-I, II-V-I, VI-V-I, I)
% Butler_M_m = [Butler_M; Butler_m];     % however, she also says that each different context gave ratings "all of which were very similar" Krumhansl (p. 29 1997)

% This is is if the context is just IV-V-I
% Butler_M = [2 0 1 0 1 1 0 2 0 1 0 1]';
% Butler_m = [2 0 1 1 0 1 0 2 1 0 0 1]';
% Butler_M_m = [Butler_M; Butler_m];

% This is is if the context is just II-V-I
% Butler_M = [1 0 2 0 1 1 0 2 0 1 0 1]';
% Butler_m = [1 0 2 1 0 1 0 2 1 0 0 1]';
% Butler_M_m = [Butler_M; Butler_m];

% This is is if the context is just VI-V-I
% Butler_M = [2 0 1 0 2 0 0 2 0 1 0 1]';
% Butler_m = [2 0 1 2 0 0 0 2 0 1 0 1]';
% Butler_M_m = [Butler_M; Butler_m];

% This is is if the context is just I
% Butler_M = [1 0 0 0 1 0 0 1 0 0 0 0]';
% Butler_m = [1 0 0 1 0 0 0 1 0 0 0 0]';
% Butler_M_m = [Butler_M; Butler_m];

%% Parncutt 89
% Saliences of all pitch classes resulting from the pure Butler counts
% Parncutt89_M = [8.37 0.15 3.59 0.4 3.43 3.62 0.06 7.24 0.35 3.37 0.68 1.7]';
% Parncutt89_m = [6.98 1.03 2.8 4.99 0.82 3.69 0.04 7.96 4.94 0.23 1.26 1.94]';
% Parncutt89_M_m = [Parncutt89_M; Parncutt89_m];

% Saliences of all pitch classes resulting from the triple tonic weighted Butler counts
Parncutt89_M = [11.51 0.15 3.69 0.5 4.53 3.94 0.08 8.32 0.47 3.69 0.7 1.7]';
Parncutt89_m = [9.58 1.05 2.84 7.11 0.82 4.31 0.04 9.42 5.58 0.27 1.28 2.06]';
Parncutt89_M_m = [Parncutt89_M; Parncutt89_m];

% Saliences of all pitch classes resulting from I counts
% Parncutt89_M_I = [1.57 0 0.05 0.05 0.55 0.16 0.01 0.54 0.06 0.16 0.01 0]';
% Parncutt89_m_I = [1.3 0.01 0.02 1.06 0 0.31 0 0.73 0.32 0.02 0.01 0.06]';
% Parncutt89_M_m_I = [Parncutt89_M_I; Parncutt89_m_I];
% Parncutt89_M_corr_I = corr(Parncutt89_M_I,y_Maj_I)
% Parncutt89_m_corr_I = corr(Parncutt89_m_I,y_min_I)

% Saliences of all pitch classes resulting from IV-V-I counts
% Parncutt89_M_IV_V_I = [2.27 0.07 0.75 0.12 0.71 1.74 0.01 2.16 0.11 0.76 0.22 0.56]';
% Parncutt89_m_IV_V_I = [2.19 0.34 0.58 1.13 0.22 1.62 0.01 2.32 1.38 0.07 0.37 0.61]';
% Parncutt89_M_m_IV_V_I = [Parncutt89_M_IV_V_I; Parncutt89_m_IV_V_I];
% Parncutt89_M_corr_IV_V_I = corr(Parncutt89_M_IV_V_I,y_Maj_IV_V_I)
% Parncutt89_m_corr_IV_V_I = corr(Parncutt89_m_IV_V_I,y_min_IV_V_I)

% Saliences of all pitch classes resulting from II-V-I counts
% Parncutt89_M_II_V_I = [1.74 0.07 1.89 0.12 0.73 1.23 0.01 2.42 0.06 0.94 0.38 0.57]';
% Parncutt89_m_II_V_I = [1.48 0.5 1.63 1.14 0.38 1.28 0.02 2.61 1.35 0.07 0.77 0.61]';
% Parncutt89_M_m_II_V_I = [Parncutt89_M_II_V_I; Parncutt89_m_II_V_I];
% Parncutt89_M_corr_II_V_I = corr(Parncutt89_M_II_V_I,y_Maj_II_V_I)
% Parncutt89_m_corr_II_V_I = corr(Parncutt89_m_II_V_I,y_min_II_V_I)

% Saliences of all pitch classes resulting from VI-V-I counts
% Parncutt89_M_VI_V_I = [2.79 0.01 0.9 0.11 1.44 0.49 0.03 2.12 0.12 1.51 0.07 0.57]';
% Parncutt89_m_VI_V_I = [2.01 0.18 0.57 1.66 0.22 0.48 0.01 2.3 1.89 0.07 0.11 0.66]';
% Parncutt89_M_m_VI_V_I = [Parncutt89_M_VI_V_I; Parncutt89_m_VI_V_I];
% Parncutt89_M_corr_VI_V_I = corr(Parncutt89_M_VI_V_I,y_Maj_VI_V_I)
% Parncutt89_m_corr_VI_V_I = corr(Parncutt89_m_VI_V_I,y_min_VI_V_I)


% Parncutt89_mean_corr = mean([Parncutt89_M_corr_I Parncutt89_m_corr_I ...
%                             Parncutt89_M_corr_IV_V_I Parncutt89_m_corr_IV_V_I ...
%                             Parncutt89_M_corr_II_V_I Parncutt89_m_corr_II_V_I ...
%                             Parncutt89_M_corr_VI_V_I Parncutt89_m_corr_VI_V_I],2)


% Saliences of all pitch classes resulting from IV-V counts
% Parncutt89_M = [0.7 0.07 0.7 0.07 0.16 1.58 0 1.62 0.05 0.6 0.21 0.56]';
% Parncutt89_m = [0.89 0.33 0.56 0.07 0.22 1.31 0.01 1.59 1.06 0.05 0.36 0.55]';
% Parncutt89_M_m = [Parncutt89_M; Parncutt89_m];


%% Parncutt 11
Parncutt11_a_M = [18 0 3 3 10 6 2 10 3 7 1 0]';
Parncutt11_a_m = [15 1 2 13 0 8 0 10 8 2 1 3]';
Parncutt11_a_M_m = [Parncutt11_a_M; Parncutt11_a_m];

Parncutt11_b_M = [226 37 121 57 196 88 45 219 64 113 70 87]'; % here I add virtual pitches to the probe as well
Parncutt11_b_m = [219 45 88 196 57 121 37 226 87 70 113 64]'; % here I add virtual pitches to the probe as well
Parncutt11_b_M_m = [Parncutt11_b_M; Parncutt11_b_m];

%% Parncutt 94
Parncutt94_M_cont = [629 246 310 310 459 374 289 459 310 395 268 246; ...
                     678 206 309 285 474 361 249 510 276 389 236 231; ...
                     645 212 310 284 461 354 251 495 275 382 240 235; ...
                     625 213 304 281 450 345 252 477 274 375 238 235];
Parncutt94_M = sum(Parncutt94_M_cont,1)';

Parncutt94_m_cont = [565 304 322 528 285 434 285 472 434 322 304 341; ...
                     614 260 305 553 248 431 228 522 428 285 263 328; ...
                     557 273 312 511 264 411 246 487 408 292 276 328; ...
                     582 264 304 531 253 415 235 499 416 286 265 327];
Parncutt94_m = sum(Parncutt94_m_cont,1)';
Parncutt94_M_m = [Parncutt94_M; Parncutt94_m];

% mean([corr(Parncutt94_M_cont(1,:)',y_Maj_I) ...
%       corr(Parncutt94_M_cont(2,:)',y_Maj_IV_V_I) ...
%       corr(Parncutt94_M_cont(3,:)',y_Maj_II_V_I) ...
%       corr(Parncutt94_M_cont(4,:)',y_Maj_VI_V_I) ...
%       corr(Parncutt94_m_cont(1,:)',y_min_I) ...
%       corr(Parncutt94_m_cont(2,:)',y_min_IV_V_I) ...
%       corr(Parncutt94_m_cont(3,:)',y_min_II_V_I) ...
%       corr(Parncutt94_m_cont(4,:)',y_min_VI_V_I)],2)

%% Empirical data transforms
K_data = [y_Maj; y_min];

% Box-Cox ('Box-Cox")
if lambda~=1
    if lambda~=0
        K_data = ((K_data.^lambda)-1)./lambda;
    else
        K_data = log(K_data);
    end
end

%% Models
num_folds = 12;
if cross_validate==0
    num_folds = 1;
    num_partitions = 1;
end

K_data_CV_fat = zeros(2,num_folds,num_partitions);
sim_data_CV_fat = zeros(2,num_folds,num_partitions);
if isequal(model,'Milne')
    betahat_CV = zeros(num_folds,5,num_partitions);
else
    betahat_CV = zeros(num_folds,2,num_partitions);
end
maj_min_corr_CV = zeros(num_partitions,1);
RLO_maj_min_CV = zeros(num_partitions,1);
RUP_maj_min_CV = zeros(num_partitions,1);
SSTCV = zeros(num_partitions,1);
SSRCV = zeros(num_partitions,1);
SSECV = zeros(num_partitions,1);
for partition = 1:num_partitions
    for fold=1:num_folds
        if cross_validate==0
            valid_ind = [];
        end
        if cross_validate==1
            % These create the partitions for a 10-fold cross-validation. 5
            % different partitions are created for 5 different 10-fold
            % cross-validations. This is to reduce the variance of the CV.
            
            % This generates random 12 X 2 partitions
            index_rand = reshape(randperm(24),[12 2]);
            
            % 5 previously generated random partitions
            if partition==1
                index_rand = [11,19;4,10;16,21;8,18;20,22;24,17;5,2;1,9;14,12;13,3;6,7;23,15;];
            elseif partition==2
                index_rand = [12,6;14,9;10,13;2,16;11,15;24,3;8,23;1,5;7,4;18,20;21,22;17,19;];
            elseif partition==3
                index_rand = [5,21;16,2;1,14;13,19;15,18;17,9;23,3;6,8;11,24;12,7;4,20;10,22;];
            elseif partition==4
                index_rand = [1,15;16,14;13,4;7,11;9,18;5,17;12,3;22,2;21,20;23,19;6,8;10,24;];
            elseif partition==5
                index_rand = [21,23;10,3;5,24;6,16;4,2;22,14;19,18;17,11;7,20;15,8;12,13;9,1;];
            elseif partition==6
                index_rand = [19,12;21,10;3,23;11,6;13,18;8,2;24,7;1,22;15,5;4,17;9,16;14,20;];
            elseif partition==7
                index_rand = [12,20;15,1;21,10;13,3;23,19;7,11;5,14;9,18;8,22;4,2;16,6;17,24;];
            elseif partition==8
                index_rand = [2,11;8,1;21,14;5,4;24,13;16,10;17,20;9,7;23,22;19,15;18,6;3,12;];
            elseif partition==9
                index_rand = [22,2;7,16;15,24;12,20;4,5;10,11;18,8;23,19;13,1;6,14;9,17;3,21;];
            elseif partition==10
                index_rand = [16,4;8,18;11,3;19,6;14,9;20,21;2,10;17,24;12,13;7,23;22,15;1,5;];
            elseif partition==11
                index_rand = [2,17;21,9;19,12;1,18;16,10;23,11;7,22;14,13;3,4;6,24;5,15;8,20;];
            elseif partition==12
                index_rand = [6,21;7,24;16,3;18,20;8,17;9,14;11,4;19,10;13,5;22,2;15,1;23,12;];
            elseif partition==13
                index_rand = [18,9;6,7;15,11;3,24;12,23;2,19;22,16;8,20;17,5;21,4;10,1;14,13;];
            elseif partition==14
                index_rand = [10,24;11,14;5,23;8,6;7,20;13,2;17,1;16,19;9,4;21,18;12,15;3,22;];
            elseif partition==15
                index_rand = [7,6;23,24;16,13;21,11;1,22;4,19;2,10;3,20;15,9;12,8;18,5;17,14;];
            elseif partition==16
                index_rand = [13,3;1,22;2,14;24,9;21,23;10,20;16,5;4,17;19,15;6,11;12,18;7,8;];
            elseif partition==17
                index_rand = [19,24;3,4;11,1;22,21;6,23;18,7;20,10;15,14;12,5;9,2;17,8;16,13;];
            elseif partition==18
                index_rand = [14,4;15,18;8,22;1,12;2,16;24,11;23,3;9,20;17,10;13,7;6,5;21,19;];
            elseif partition==19
                index_rand = [8,15;17,6;14,3;4,23;11,7;10,19;1,2;9,24;16,18;20,22;5,12;21,13;];
            elseif partition==20
                index_rand = [5,24;17,13;10,1;6,7;14,22;12,15;9,8;11,20;3,2;23,4;21,18;19,16;];
            end
            index_Maj = find(index_rand'<13);
            index_min = find(index_rand'>12);
            valid_ind = index_rand(fold,:); % These are the indices of the validation set
        end
        K_data_valid = K_data(valid_ind);
        K_data_train = K_data;
        K_data_train(valid_ind) = [];
        
        if ~isequal(model,'Milne')
            if isequal(model,'Null')
                model = Null_M_m;
            elseif isequal(model,'DumbTriad')
                model = Triad_M_m;
            elseif isequal(model,'KrumPrev')
                model = Krum_prev_M_m;
            elseif isequal(model,'KrumCons')
                model = Krum_cons_M_m;
            elseif isequal(model,'Smith')
                model = Smith_M_m;
            elseif isequal(model,'Ler')
                model = Lerdahl_M_m;
            elseif isequal(model,'Butler')
                model = Butler_M_m;
            elseif isequal(model,'P89')
                model = Parncutt89_M_m;
            elseif isequal(model,'P11a')
                model = Parncutt11_a_M_m;
            elseif isequal(model,'P11b')
                model = Parncutt11_b_M_m;
            elseif isequal(model,'P94')
                model = Parncutt94_M_m;
            end
            model_train = model;
            if cross_validate==1
                model_train(valid_ind) = [];
                betahat = regress(K_data_train,[ones(22,1) model_train])';
            else
                betahat = regress(K_data_train,[ones(24,1) model_train])';
            end
            sim_data = betahat(1) + betahat(2)*model;
        elseif isequal(model,'Milne')
            %% Optimize nonlinear spectral similarity model
            % The parameters are: 
            % rolloff, sigma, omega, intercept, slope
            % Some useful start values to speed up optimization:
            param_vals = [0.6778 6.2078 0.4783 2.4316 4.7095];
            lb =         [0      1      0     -inf   -inf   ];
            ub =         [inf    inf    1      inf    inf   ];
            % Fixing the values of rho and sigma as optimized in other
            % experiment (a test of parameter sensitivity):
%             lb = [0.6    1        0    -inf  -inf ];
%             ub = [0.6    inf      1     inf   inf ];
%             lb = [0.42   10.28    0    -inf  -inf ];
%             ub = [0.42   10.28    1     inf   inf ];
%             lb = [0.42   10.28    0    -inf  -inf ];
%             ub = [0.42   10.28    1     inf   inf ];
            options = optimset('Algorithm','active-set','TolFun',1e-8, ...
                'Diagnostics','on','Display','iter-detailed', ...
                'MaxFunEvals',2400,'FinDiffType','central');
            [betahat,fval,exitflag,~,~,~,~] ...
                = fmincon(@Probe_tone_spectral_sim_model,param_vals,[],[],[],[],lb,ub,[],options)
            
            %% Modelled data and statistical tests
            betahat = [betahat(4) betahat(5) betahat(1:3)];
            sim_data = betahat(1) + betahat(2)*spec_sim;
        end

        %% Plots and analysis
        figure(11)
        scatter(K_data,sim_data)
        axis([1 7 1 7])
        set(gca,'PlotBoxAspectRatio', [1 1 1])
        
        if cross_validate==1;
            betahat_CV(fold,:,partition) = betahat;
            sim_data_valid = sim_data(valid_ind);
            K_data_CV_fat(:,fold,partition) = K_data_valid;
            sim_data_CV_fat(:,fold,partition) = sim_data_valid;
        end
    end %folds
    
    %% CV Score
    if cross_validate==1
        K_data_CV_part = K_data_CV_fat(:,:,partition);
        K_data_CV = K_data_CV_part(:);
        sim_data_CV_part = sim_data_CV_fat(:,:,partition);
        sim_data_CV = sim_data_CV_part(:);
        
        [R_maj_min,P_maj_min,RLO_maj_min,RUP_maj_min] = corrcoef(K_data_CV, sim_data_CV);
        maj_min_corr_CV(partition,:) = R_maj_min(1,2);
        RLO_maj_min_CV(partition,:) = RLO_maj_min(1,2);
        RUP_maj_min_CV(partition,:) = RUP_maj_min(1,2);
        
        SSTCV(partition,:) = sum((K_data_CV - mean(K_data_CV)).^2);
        SSRCV(partition,:) = sum((sim_data_CV - mean(K_data_CV)).^2);
        SSECV(partition,:) = sum((K_data_CV - sim_data_CV).^2);
        
        MSECV = SSECV/24;
        RMSECV = sqrt(MSECV);
    else
        [R,P,RLO,RUP] = corrcoef(K_data(1:12), sim_data(1:12),'alpha',0.05);
        maj_corr = R(1,2);
        maj_corr_CI = [RLO(1,2),RUP(1,2)];

        [R,P,RLO,RUP] = corrcoef(K_data(13:24), sim_data(13:24),'alpha',0.05);
        min_corr = R(1,2);
        min_corr_CI = [RLO(1,2),RUP(1,2)];

        [R,P,RLO,RUP] = corrcoef(K_data, sim_data,'alpha',0.05);
        maj_min_corr = R(1,2);
        maj_min_corr_CI = [RLO(1,2),RUP(1,2)];
                
        SST_all = sum((K_data - mean(K_data)).^2);
        SSR_all = sum((sim_data - mean(K_data)).^2);
        SSE_all = sum((K_data - sim_data).^2);

        SST_Maj = sum((K_data(1:12) - mean(K_data(1:12))).^2);
        SSR_Maj = sum((sim_data(1:12) - mean(K_data(1:12))).^2);
        SSE_Maj = sum((K_data(1:12) - sim_data(1:12)).^2);
        
        SST_min = sum((K_data(13:24) - mean(K_data(13:24))).^2);
        SSR_min = sum((sim_data(13:24) - mean(K_data(13:24))).^2);
        SSE_min = sum((K_data(13:24) - sim_data(13:24)).^2);

        MSR_all = SSR_all/24;
        MSE_all = SSE_all/24;
        RMSE_all = sqrt(MSE_all);

        MSR_Maj = SSR_Maj/12;
        MSE_Maj = SSE_Maj/12;
        RMSE_Maj = sqrt(MSE_Maj);
        
        MSR_min = SSR_min/12;
        MSE_min = SSE_min/12;
        RMSE_min = sqrt(MSE_min);
    end
end

if cross_validate==1
%     betahat_CV
    maj_min_corr_CV
%     RLO_maj_min_CV
%     RUP_maj_min_CV
    
    mean_betahat_CV = mean(mean(betahat_CV,1),3)
    mean_maj_min_corr_CV = mean(maj_min_corr_CV,1)
    mean_maj_min_CV_CI = [mean(RLO_maj_min_CV) mean(RUP_maj_min_CV)]
else
    betahat
    maj_min_corr
    maj_min_corr_CI    
    maj_corr
    maj_corr_CI
    min_corr
    min_corr_CI

% BIC_5 = 24*log(SSE_all) - 24*log(24) + log(24)*5 % BIC for model with k=5 (Milne b and c)
% BIC_3 = 24*log(SSE_all) - 24*log(24) + log(24)*3 % BIC for models with k=3 (Parncutt 94)
% BIC_2 = 24*log(SSE_all) - 24*log(24) + log(24)*2 % BIC for models with k=2 (standard linear models)
    
end

%% Residuals
Residuals = sim_data - K_data;
Resid_pred = [sim_data (Residuals)];
figure(12) % this plots residuals against predicted value to test homogeneity
scatter(Resid_pred(:,1),Resid_pred(:,2),'DisplayName','Resid_pred(:,1) vs. Resid_pred(:,2)','YDataSource','Resid_pred(:,2)');figure(gcf)

figure(13) % this gives a histogram of residuals to test normality
hist(Residuals)

figure(14) % Q-Q plot to test for normality of residuals
qqplot(normrnd(0,1,1000,1),Residuals)
set(gca,'PlotBoxAspectRatio', [1 1 1])

figure(15), subplot(2,2,1)
bar(K_data(1:12),'DisplayName','y_Maj');figure(gcf)
axis([0,13,1,7])
% set(gca,'XTick',1:q*100/period:q+1)
set(gca,'XTickLabel',0:13)
set(gca,'PlotBoxAspectRatio',[1.618034,1,1])

figure(15), subplot(2,2,2)
bar(K_data(13:24),'DisplayName','y_min');figure(gcf)
axis([0,13,1,7])
set(gca,'XTickLabel',0:12)
set(gca,'PlotBoxAspectRatio',[1.618034,1,1])

figure(15), subplot(2,2,3)
bar(sim_data(1:12),'DisplayName','model_Maj');figure(gcf)
axis([0,13,1,7])
set(gca,'XTickLabel',0:12)
set(gca,'PlotBoxAspectRatio',[1.618034,1,1])

figure(15), subplot(2,2,4)
bar(sim_data(13:24),'DisplayName','model_min');figure(gcf)
axis([0,13,1,7])
set(gca,'XTickLabel',0:12)
set(gca,'PlotBoxAspectRatio',[1.618034,1,1])

figure(16)
bar(K_data(1:12),'FaceColor',[0.2 0.2 0.2],'DisplayName','y_Maj');figure(gcf)
axis([0,13,1,7])
% set(gca,'XTick',1:q*100/period:q+1)
set(gca,'FontSize',15)
set(gca,'XTickLabel',0:13)
set(gca,'PlotBoxAspectRatio',[1.618034,1,1])
xlabel('Pitch class (relative to tonal context) of probe tone')
ylabel({'Mean rating of probe tone''s fit with' 'previously established major context'})

figure(17)
bar(K_data(13:24),'FaceColor',[0.2 0.2 0.2],'DisplayName','y_min');figure(gcf)
axis([0,13,1,7])
set(gca,'FontSize',15)
set(gca,'XTickLabel',0:12)
set(gca,'PlotBoxAspectRatio',[1.618034,1,1])
xlabel('Pitch class (relative to tonal context) of probe tone')
ylabel({'Mean rating of probe tone''s fit with' 'previously established minor context'})

% toc

%{
%% Additional bits and pieces

%% Krumhansl & Shepard scale data
% Approximation of the graph in Krumhansl & Shepard 1979, Group 1

% scale_data = [6.7 2.6 4.6 2.4 4 3.4 2.8 5 2.3 3.8 2.5 4.5];
% model = [0.68 0.05 0.45 0.06 0.44 0.43 0.15 0.58 0.06 0.5 0.15 0.38];
% 
% scale_corr = corr(scale_data', model')

%% Intercorrelation table
% These are the optimized non-cross-correlated models' predictions
Model_a = [0.727491394353878,0.0108110933032899,0.281818561060869,0.0318876097793087,0.537983242689168,0.239371536460207,0.0458417577889489,0.710436504884408,0.0489424992487780,0.246787896575210,0.0532581179039513,0.221395761444018,0.710436504884408,0.0458417577889489,0.239371536460207,0.537983242689168,0.0318876097793087,0.281818561060869,0.0108110933032899,0.727491394353877,0.221395761444018,0.0532581179039513,0.246787896575210,0.0489424992487780]';
Model_b = [0.824039916588250,0.0113166282429852,0.214571492479311,0.0397790474908824,0.465246932688757,0.251426611915019,0.0342116476493812,0.644572568117804,0.0709266214901930,0.181563020190466,0.0566621434001307,0.168051622326387,0.808237543792734,0.0323014133934872,0.186846139536973,0.455929841943957,0.0515158643013205,0.281384719658315,0.0134449745251776,0.667411494990442,0.184976303367907,0.0446975712476525,0.194668997960602,0.0551469554783941]';
Model_c = [0.870735249920137,0.00703913717543092,0.209100987315360,0.0274411181556309,0.400867918336043,0.286565610461689,0.0307673884055237,0.625575438301506,0.0530931325304822,0.165863164528487,0.0600578979262442,0.147317268598803,0.728215664189407,0.0465722114634279,0.154171047364734,0.633694808579838,0.0340182117296618,0.286595523634461,0.0119193712914539,0.553750190633269,0.249452040058836,0.0342211423890624,0.268408000706127,0.0448084634908144]';

All_data = [K_data Triad_M_m Krum_prev_M_m Krum_cons_M_m Smith_M_m Lerdahl_M_m Butler_M_m Parncutt89_M_m Parncutt88_M_m Parncutt88_b_M_m Parncutt94_M_m Model_a Model_b Model_c];
model_inter_corr = corr(All_data)

% Contexts
Maj_contexts = [y_Maj_I y_Maj_IV_V_I y_Maj_II_V_I y_Maj_VI_V_I];
Min_contexts = [y_min_I y_min_IV_V_I y_min_II_V_I y_min_VI_V_I];
maj_cont_corr = corr(Maj_contexts)
min_cont_corr = corr(Min_contexts)

% Lerdahl intercorrelations with each context
Data_L = [[y_Maj_I; y_min_I] [y_Maj_IV_V_I; y_min_IV_V_I] [y_Maj_II_V_I; y_min_II_V_I] [y_Maj_VI_V_I; y_min_VI_V_I] Lerdahl_M_m];
model_inter_corr_L = corr(Data_L)

% Parncutt 11a intercorrelations with each context
Data_P = [[y_Maj_I; y_min_I] [y_Maj_IV_V_I; y_min_IV_V_I] [y_Maj_II_V_I; y_min_II_V_I] [y_Maj_VI_V_I; y_min_VI_V_I] Parncutt88_M_m];
model_inter_corr_P = corr(Data_P)

% Milne 14c intercorrelations with each context
Data_M = [[y_Maj_I; y_min_I] [y_Maj_IV_V_I; y_min_IV_V_I] [y_Maj_II_V_I; y_min_II_V_I] [y_Maj_VI_V_I; y_min_VI_V_I] Model_c];
model_inter_corr_M = corr(Data_M)

% Lerdahl intercorrelations with each context
Data_Maj_L = [Lerdahl_M_m(1:12) y_Maj_I y_Maj_IV_V_I y_Maj_II_V_I y_Maj_VI_V_I];
model_inter_corr_Maj_L = corr(Data_Maj_L)

Data_min_L = [Lerdahl_M_m(13:24) y_min_I y_min_IV_V_I y_min_II_V_I y_min_VI_V_I];
model_inter_corr_min_L = corr(Data_min_L)

% Parncutt 11a intercorrelations with each context
Data_Maj_P = [Parncutt88_M_m(1:12) y_Maj_I y_Maj_IV_V_I y_Maj_II_V_I y_Maj_VI_V_I];
model_inter_corr_Maj_P = corr(Data_Maj_P)

Data_min_P = [Parncutt88_M_m(13:24) y_min_I y_min_IV_V_I y_min_II_V_I y_min_VI_V_I];
model_inter_corr_min_P = corr(Data_min_P)

% Milne intercorrelations with each context
Data_Maj_M = [Model_c(1:12) y_Maj_I y_Maj_IV_V_I y_Maj_II_V_I y_Maj_VI_V_I];
model_inter_corr_Maj_M = corr(Data_Maj_M)

Data_min_M = [Model_c(13:24) y_min_I y_min_IV_V_I y_min_II_V_I y_min_VI_V_I];
model_inter_corr_min_M = corr(Data_min_M)

%}


