%% This function fits the model of melodic affinity---described in Milne, 
% Laney, Sharp (2013)---to emprical data collected from 44 participants in
% March and April 2012. It is called by the routine
% Affinity_data_analysis.m. Other necessary functions are
% (convnfft by Bruno Luong, DERIVESTsuite by John D'Errico, myBinomTest.m by 
% Matthew Nelson), which can be downloaded from
% http://www.mathworks.com/matlabcentral/fileexchange/

% References Milne, Laney, Sharp (2013), Melodic affinity: Nature or
% nurture, Music Perception.

% Author: Andrew J. Milne
% Version: 1.0

function Deviance = Affinity_model_final(param_vals)

global win_len q period d_x_H x_pc_H d_y_H y_pc_H d_y_S y_pc_S ...
    constant variables simplify ...
    spsim_self vpsim_self harm_self ...
    successes_train trials_train valid_stims beta_hat ...
    spsim vpsim harm

%% Parameter values
par_sigma_S = param_vals(1);
par_rolloff_S = param_vals(2);
par_sigma_H = param_vals(3);
par_rolloff_H = param_vals(4);
b = param_vals(5:8)';
b_const = param_vals(5);
b_spsim = param_vals(6);
b_vpsim = param_vals(7);
b_harm = param_vals(8);
if simplify==1
    par_sigma_S = par_sigma_H;
    par_rolloff_S = par_rolloff_H;
end

%% Test for changed parameters
if ~exist('rolloff_S','var') || ~isequal(par_rolloff_S,rolloff_S)
    rolloff_S_chg = 1;
else
    rolloff_S_chg = 0;
end
rolloff_S = par_rolloff_S;

if ~exist('rolloff_H','var') || ~isequal(par_rolloff_H,rolloff_H)
    rolloff_H_chg = 1;
else
    rolloff_H_chg = 0;
end
rolloff_H = par_rolloff_H;

if ~exist('sigma_H','var') || ~isequal(par_sigma_H,sigma_H)
    sigma_H_chg = 1;
else
    sigma_H_chg = 0;
end
sigma_H = par_sigma_H;

if ~exist('sigma_S','var') || ~isequal(par_sigma_S,sigma_S)
    sigma_S_chg = 1;
else
    sigma_S_chg = 0;
end
sigma_S = par_sigma_S;

%% Spectral saliences for x (template) and y (test tones)
if rolloff_H_chg==1
%    x_w_H = 1./(1:12).^rolloff_H;
    x_w_H = ((1:12).^-rolloff_H).*min(1,12^rolloff_H);
    y_w_H = x_w_H;
end

%% HARMONICITY MODEL

%% Embed and smooth pitches
if sigma_H_chg==1
    windowlength_H = ceil(win_len*sigma_H*q/period);
    n_H = -(windowlength_H-1)/2 : (windowlength_H-1)/2;
    gwin_H = exp(-(n_H.^2) / (2 * (sigma_H*q/period)^2));
    gwin_H = gwin_H/norm(gwin_H,1);
end

%% Pitch class response matrix X
if sigma_H_chg==1 || rolloff_H==1
    X_pcr_ij_H = zeros(d_x_H,q);
    for i = 1:d_x_H
        X_pcr_ij_H(i,x_pc_H(i)+1) = x_w_H(i);
        if length(gwin_H) > 1
            X_pcr_ij_H(i,:) = real(ifft(...
                fft(X_pcr_ij_H(i,:), q).*conj(fft(gwin_H, q)), q));
        end
    end
    X_pcr_ij_H(X_pcr_ij_H < 0.0000001) = 0;
    
    %% Salience vector x
    x_1_s_j_H = 1 - prod(1 - X_pcr_ij_H, 1);
    x_s_j_H = circshift(x_1_s_j_H,[0,floor(windowlength_H/2)-1]);
    %x_s_j = x_1_s_j;
    
    %figure(1)
    %plot(x_s_j)

    %% Pitch class response matrix Y
    Y_pcr_ij_H = zeros(d_y_H,q,11);
    for k = 1:11
        for i = 1:d_y_H
            Y_pcr_ij_H(i,y_pc_H(k,i)+1,k) = y_w_H(i);
            if length(gwin_H) > 1
                Y_pcr_ij_H(i,:,k) ...
            = real(ifft(fft(Y_pcr_ij_H(i,:,k), q).*conj(fft(gwin_H, q)), q));
            end
        end
    end
    Y_pcr_ij_H(Y_pcr_ij_H < 0.0000001) = 0;

    %% Salience vectors y (for each n-TET)
    y_1_s_j_H = 1 - prod(1 - Y_pcr_ij_H, 1);
    y_s_j_H = circshift(y_1_s_j_H,[0,floor(windowlength_H/2)-1]);
    %y_s_j = y_1_s_j;
    y_s_j_H = squeeze(permute(y_s_j_H, [3 2 1]));
    
    %figure(2)
    %surface(y_s_j)

    %% Max of cross correlations for each of the 11 different spectra
    % convnfft comes from Bruno Luong (it's faster than conv)
    virtp=zeros(11,q*2-1);
    for k=1:11
        virtp(k,:) = convnfft(x_s_j_H,fliplr(y_s_j_H(k,:)),'full');
    end
    virtp;
    maxes = max(virtp,[],2)';

end

%% SPECTRAL PITCH DISTANCE MODEL

%% Spectral saliences for y (test tones)
if rolloff_S_chg==1
%    y_w_S = 1./(1:12).^rolloff_S;
    y_w_S = ((1:12).^-rolloff_S).*min(1,12^rolloff_S);
end

%% Embed and smooth pitches
if sigma_S_chg==1
    windowlength_S = ceil(win_len*sigma_S*q/period);
    n_S = -(windowlength_S-1)/2 : (windowlength_S-1)/2;
    gwin_S = exp(-(n_S.^2) / (2 * (sigma_S*q/period)^2));
    gwin_S = gwin_S/norm(gwin_S,1);
end

%% Pitch class response matrix Y
if sigma_S_chg==1 || rolloff_S_chg==1
Y_pcr_ij_S = zeros(d_y_S,q,11);
for k = 1:11
    for i = 1:d_y_S
        Y_pcr_ij_S(i,y_pc_S(k,i)+1,k) = y_w_S(i);
        if length(gwin_S) > 1
            Y_pcr_ij_S(i,:,k) = real(ifft(...
                fft(Y_pcr_ij_S(i,:,k), q).*conj(fft(gwin_S, q)), q));
        end
    end
end
Y_pcr_ij_S(Y_pcr_ij_S < 0.0000001) = 0;

%% Salience vectors y (for each n-TET)
y_1_s_j_S = 1 - prod(1 - Y_pcr_ij_S, 1);
y_s_j_S = circshift(y_1_s_j_S,[0,floor(windowlength_S/2)-1]);
%y_s_j = y_1_s_j;
y_s_j_S = squeeze(permute(y_s_j_S, [3 2 1]));

%% Calculate the spectral pitch distance for all upwards intervals 
% (columns) in each tuning for each spectral tuning (rows)

% The faster code for circular shift is taken Matt J:
% http://www.mathworks.com/matlabcentral/newsreader/view_thread/279310
spec_dist_3 = zeros(11,3);
step_size = 400;
steps = round(0:step_size:2*step_size);
for j=1:11
    y_s_S_tun = y_s_j_S(j,:);
    y_s_S_tun = y_s_S_tun(ones(3,1),:)';
    col = y_s_S_tun(:,1);
    for i=1:3
        step_i = steps(i);
        col(:) = y_s_S_tun(:,i);
        y_s_S_tun(1:end-step_i,i) = col(step_i+1:end);
        y_s_S_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    y_s_S_shft = y_s_S_tun';
    y_s_S_shft_1 = y_s_S_shft(1,:);
    for i=1:3
        y_s_S_shft_i = y_s_S_shft(i,:);
        spec_dist_3(j,i)...
            = 1-(y_s_S_shft_1*y_s_S_shft_i' ...
            / sqrt((y_s_S_shft_1*y_s_S_shft_1')*(y_s_S_shft_i*y_s_S_shft_i')));
    end
end
spec_dist_3;

spec_dist_4 = zeros(11,4);
step_size = 300;
steps = round(0:step_size:3*step_size);
for j=1:11
    y_s_S_tun = y_s_j_S(j,:);
    y_s_S_tun = y_s_S_tun(ones(4,1),:)';
    col = y_s_S_tun(:,1);
    for i=1:4
        step_i = steps(i);
        col(:) = y_s_S_tun(:,i);
        y_s_S_tun(1:end-step_i,i) = col(step_i+1:end);
        y_s_S_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    y_s_S_shft = y_s_S_tun';
    y_s_S_shft_1 = y_s_S_shft(1,:);
    for i=1:4
        y_s_S_shft_i = y_s_S_shft(i,:);
        spec_dist_4(j,i)...
            = 1-(y_s_S_shft_1*y_s_S_shft_i' ...
            / sqrt((y_s_S_shft_1*y_s_S_shft_1')*(y_s_S_shft_i*y_s_S_shft_i')));
    end
end
spec_dist_4;

spec_dist_5 = zeros(11,5);
step_size = 240;
steps = round(0:step_size:4*step_size);
for j=1:11
    y_s_S_tun = y_s_j_S(j,:);
    y_s_S_tun = y_s_S_tun(ones(5,1),:)';
    col = y_s_S_tun(:,1);
    for i=1:5
        step_i = steps(i);
        col(:) = y_s_S_tun(:,i);
        y_s_S_tun(1:end-step_i,i) = col(step_i+1:end);
        y_s_S_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    y_s_S_shft = y_s_S_tun';
    y_s_S_shft_1 = y_s_S_shft(1,:);
    for i=1:5
        y_s_S_shft_i = y_s_S_shft(i,:);
        spec_dist_5(j,i)...
            = 1-(y_s_S_shft_1*y_s_S_shft_i' ...
            / sqrt((y_s_S_shft_1*y_s_S_shft_1')*(y_s_S_shft_i*y_s_S_shft_i')));
    end
end
spec_dist_5;

spec_dist_7 = zeros(11,7);
step_size = 1200/7;
steps = round(0:step_size:6*step_size);
for j=1:11
    y_s_S_tun = y_s_j_S(j,:);
    y_s_S_tun = y_s_S_tun(ones(7,1),:)';
    col = y_s_S_tun(:,1);
    for i=1:7
        step_i = steps(i);
        col(:) = y_s_S_tun(:,i);
        y_s_S_tun(1:end-step_i,i) = col(step_i+1:end);
        y_s_S_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    y_s_S_shft = y_s_S_tun';
    y_s_S_shft_1 = y_s_S_shft(1,:);
    for i=1:7
        y_s_S_shft_i = y_s_S_shft(i,:);
        spec_dist_7(j,i)...
            = 1-(y_s_S_shft_1*y_s_S_shft_i' ...
            / sqrt((y_s_S_shft_1*y_s_S_shft_1')*(y_s_S_shft_i*y_s_S_shft_i')));
    end
end
spec_dist_7;

spec_dist_10 = zeros(11,10);
step_size = 120;
steps = round(0:step_size:9*step_size);
for j=1:11
    y_s_S_tun = y_s_j_S(j,:);
    y_s_S_tun = y_s_S_tun(ones(10,1),:)';
    col = y_s_S_tun(:,1);
    for i=1:10
        step_i = steps(i);
        col(:) = y_s_S_tun(:,i);
        y_s_S_tun(1:end-step_i,i) = col(step_i+1:end);
        y_s_S_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    y_s_S_shft = y_s_S_tun';
    y_s_S_shft_1 = y_s_S_shft(1,:);
    for i=1:10
        y_s_S_shft_i = y_s_S_shft(i,:);
        spec_dist_10(j,i)...
            = 1-(y_s_S_shft_1*y_s_S_shft_i' ...
            / sqrt((y_s_S_shft_1*y_s_S_shft_1')*(y_s_S_shft_i*y_s_S_shft_i')));
    end
end
spec_dist_10;

spec_dist_11 = zeros(11,11);
step_size = 1200/11;
steps = round(0:step_size:10*step_size);
for j=1:11
    y_s_S_tun = y_s_j_S(j,:);
    y_s_S_tun = y_s_S_tun(ones(11,1),:)';
    col = y_s_S_tun(:,1);
    for i=1:11
        step_i = steps(i);
        col(:) = y_s_S_tun(:,i);
        y_s_S_tun(1:end-step_i,i) = col(step_i+1:end);
        y_s_S_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    y_s_S_shft = y_s_S_tun';
    y_s_S_shft_1 = y_s_S_shft(1,:);
    for i=1:11
        y_s_S_shft_i = y_s_S_shft(i,:);
        spec_dist_11(j,i)...
            = 1-(y_s_S_shft_1*y_s_S_shft_i' ...
            / sqrt((y_s_S_shft_1*y_s_S_shft_1')*(y_s_S_shft_i*y_s_S_shft_i')));
    end
end
spec_dist_11;

spec_dist_12 = zeros(11,12);
step_size = 100;
steps = round(0:step_size:11*step_size);
for j=1:11
    y_s_S_tun = y_s_j_S(j,:);
    y_s_S_tun = y_s_S_tun(ones(12,1),:)';
    col = y_s_S_tun(:,1);
    for i=1:12
        step_i = steps(i);
        col(:) = y_s_S_tun(:,i);
        y_s_S_tun(1:end-step_i,i) = col(step_i+1:end);
        y_s_S_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    y_s_S_shft = y_s_S_tun';
    y_s_S_shft_1 = y_s_S_shft(1,:);
    for i=1:12
        y_s_S_shft_i = y_s_S_shft(i,:);
        spec_dist_12(j,i)...
            = 1-(y_s_S_shft_1*y_s_S_shft_i' ...
            / sqrt((y_s_S_shft_1*y_s_S_shft_1')*(y_s_S_shft_i*y_s_S_shft_i')));
    end
end
spec_dist_12;

spec_dist_13 = zeros(11,13);
step_size = 1200/13;
steps = round(0:step_size:12*step_size);
for j=1:11
    y_s_S_tun = y_s_j_S(j,:);
    y_s_S_tun = y_s_S_tun(ones(13,1),:)';
    col = y_s_S_tun(:,1);
    for i=1:13
        step_i = steps(i);
        col(:) = y_s_S_tun(:,i);
        y_s_S_tun(1:end-step_i,i) = col(step_i+1:end);
        y_s_S_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    y_s_S_shft = y_s_S_tun';
    y_s_S_shft_1 = y_s_S_shft(1,:);
    for i=1:13
        y_s_S_shft_i = y_s_S_shft(i,:);
        spec_dist_13(j,i)...
            = 1-(y_s_S_shft_1*y_s_S_shft_i' ...
            / sqrt((y_s_S_shft_1*y_s_S_shft_1')*(y_s_S_shft_i*y_s_S_shft_i')));
    end
end
spec_dist_13;

spec_dist_15 = zeros(11,15);
step_size = 80;
steps = round(0:step_size:14*step_size);
for j=1:11
    y_s_S_tun = y_s_j_S(j,:);
    y_s_S_tun = y_s_S_tun(ones(15,1),:)';
    col = y_s_S_tun(:,1);
    for i=1:15
        step_i = steps(i);
        col(:) = y_s_S_tun(:,i);
        y_s_S_tun(1:end-step_i,i) = col(step_i+1:end);
        y_s_S_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    y_s_S_shft = y_s_S_tun';
    y_s_S_shft_1 = y_s_S_shft(1,:);
    for i=1:15
        y_s_S_shft_i = y_s_S_shft(i,:);
        spec_dist_15(j,i)...
            = 1-(y_s_S_shft_1*y_s_S_shft_i' ...
            / sqrt((y_s_S_shft_1*y_s_S_shft_1')*(y_s_S_shft_i*y_s_S_shft_i')));
    end
end
spec_dist_15;

spec_dist_16 = zeros(11,16);
step_size = 75;
steps = round(0:step_size:15*step_size);
for j=1:11
    y_s_S_tun = y_s_j_S(j,:);
    y_s_S_tun = y_s_S_tun(ones(16,1),:)';
    col = y_s_S_tun(:,1);
    for i=1:16
        step_i = steps(i);
        col(:) = y_s_S_tun(:,i);
        y_s_S_tun(1:end-step_i,i) = col(step_i+1:end);
        y_s_S_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    y_s_S_shft = y_s_S_tun';
    y_s_S_shft_1 = y_s_S_shft(1,:);
    for i=1:16
        y_s_S_shft_i = y_s_S_shft(i,:);
        spec_dist_16(j,i)...
            = 1-(y_s_S_shft_1*y_s_S_shft_i' ...
            / sqrt((y_s_S_shft_1*y_s_S_shft_1')*(y_s_S_shft_i*y_s_S_shft_i')));
    end
end
spec_dist_16;

spec_dist_17 = zeros(11,17);
step_size = 1200/17;
steps = round(0:step_size:16*step_size);
for j=1:11
    y_s_S_tun = y_s_j_S(j,:);
    y_s_S_tun = y_s_S_tun(ones(17,1),:)';
    col = y_s_S_tun(:,1);
    for i=1:17
        step_i = steps(i);
        col(:) = y_s_S_tun(:,i);
        y_s_S_tun(1:end-step_i,i) = col(step_i+1:end);
        y_s_S_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    y_s_S_shft = y_s_S_tun';
    y_s_S_shft_1 = y_s_S_shft(1,:);
    for i=1:17
        y_s_S_shft_i = y_s_S_shft(i,:);
        spec_dist_17(j,i)...
            = 1-(y_s_S_shft_1*y_s_S_shft_i' ...
            / sqrt((y_s_S_shft_1*y_s_S_shft_1')*(y_s_S_shft_i*y_s_S_shft_i')));
    end
end
spec_dist_17;

%% Probs of intervals (-1200, +1200) cents due to generator moves
% Probability for each gen move
gen_moves = [-6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6]';
tri_dist_gen = (1-abs(gen_moves)/6.5)/sum(1-abs(gen_moves)/6.5);
% plot(tri_dist_gen)

% Calculate the scale degree each gen move maps to
 scale_deg_3 = mod(gen_moves,3);
 scale_deg_4 = mod(gen_moves,4);
 scale_deg_5 = mod(3*gen_moves,5);
 scale_deg_7 = mod(4*gen_moves,7);
scale_deg_10 = mod(3*gen_moves,10);
scale_deg_11 = mod(3*gen_moves,11);
scale_deg_12 = mod(7*gen_moves,12);
scale_deg_13 = mod(4*gen_moves,13);
scale_deg_15 = mod(4*gen_moves,15);
scale_deg_16 = mod(5*gen_moves,16);
scale_deg_17 = mod(10*gen_moves,17);

% This calculates the probability of each interval over -1200 to 
% +1200 cents as a function of gen moves

scale_deg_gen_3 = zeros(13,3);
for i=1:3
    scale_deg_gen_3(scale_deg_3==i-1,i) = 1;
end
scale_deg_gen_3;
scale_int_gen_3 = sum(bsxfun(@times, tri_dist_gen, scale_deg_gen_3));
full_int_gen_3 = [scale_int_gen_3(1,2:end) scale_int_gen_3];

scale_deg_gen_4 = zeros(13,4);
for i=1:4
    scale_deg_gen_4(scale_deg_4==i-1,i) = 1;
end
scale_deg_gen_4;
scale_int_gen_4 = sum(bsxfun(@times, tri_dist_gen, scale_deg_gen_4));
full_int_gen_4 = [scale_int_gen_4(1,2:end) scale_int_gen_4];

scale_deg_gen_5 = zeros(13,5);
for i=1:5
    scale_deg_gen_5(scale_deg_5==i-1,i) = 1;
end
scale_deg_gen_5;
scale_int_gen_5 = sum(bsxfun(@times, tri_dist_gen, scale_deg_gen_5));
full_int_gen_5 = [scale_int_gen_5(1,2:end) scale_int_gen_5];

scale_deg_gen_7 = zeros(13,7);
for i=1:7
    scale_deg_gen_7(scale_deg_7==i-1,i) = 1;
end
scale_deg_gen_7;
scale_int_gen_7 = sum(bsxfun(@times, tri_dist_gen, scale_deg_gen_7));
full_int_gen_7 = [scale_int_gen_7(1,2:end) scale_int_gen_7];

scale_deg_gen_10 = zeros(13,10);
for i=1:10
    scale_deg_gen_10(scale_deg_10==i-1,i) = 1;
end
scale_deg_gen_10;
scale_int_gen_10 = sum(bsxfun(@times, tri_dist_gen, scale_deg_gen_10));
full_int_gen_10 = [scale_int_gen_10(1,2:end) scale_int_gen_10];

scale_deg_gen_11 = zeros(13,11);
for i=1:11
    scale_deg_gen_11(scale_deg_11==i-1,i) = 1;
end
scale_deg_gen_11;
scale_int_gen_11 = sum(bsxfun(@times, tri_dist_gen, scale_deg_gen_11));
full_int_gen_11 = [scale_int_gen_11(1,2:end) scale_int_gen_11];

scale_deg_gen_12 = zeros(13,12);
for i=1:12
    scale_deg_gen_12(scale_deg_12==i-1,i) = 1;
end
scale_deg_gen_12;
scale_int_gen_12 = sum(bsxfun(@times, tri_dist_gen, scale_deg_gen_12));
full_int_gen_12 = [scale_int_gen_12(1,2:end) scale_int_gen_12];

scale_deg_gen_13 = zeros(13,13);
for i=1:13
    scale_deg_gen_13(scale_deg_13==i-1,i) = 1;
end
scale_deg_gen_13;
scale_int_gen_13 = sum(bsxfun(@times, tri_dist_gen, scale_deg_gen_13));
full_int_gen_13 = [scale_int_gen_13(1,2:end) scale_int_gen_13];

scale_deg_gen_15 = zeros(13,15);
for i=1:15
    scale_deg_gen_15(scale_deg_15==i-1,i) = 1;
end
scale_deg_gen_15;
scale_int_gen_15 = sum(bsxfun(@times, tri_dist_gen, scale_deg_gen_15));
full_int_gen_15 = [scale_int_gen_15(1,2:end) scale_int_gen_15];

scale_deg_gen_16 = zeros(13,16);
for i=1:16
    scale_deg_gen_16(scale_deg_16==i-1,i) = 1;
end
scale_deg_gen_16;
scale_int_gen_16 = sum(bsxfun(@times, tri_dist_gen, scale_deg_gen_16));
full_int_gen_16 = [scale_int_gen_16(1,2:end) scale_int_gen_16];

scale_deg_gen_17 = zeros(13,17);
for i=1:17
    scale_deg_gen_17(scale_deg_17==i-1,i) = 1;
end
scale_deg_gen_17;
scale_int_gen_17 = sum(bsxfun(@times, tri_dist_gen, scale_deg_gen_17));
full_int_gen_17 = [scale_int_gen_17(1,2:end) scale_int_gen_17];

%% Probs of intervals (-1200, +1200) cents due to interval size
scale_int_size_3 = zeros(1,3);
for i=1:3
    scale_int_size_3(1,i) = 1-abs(400*(i-1))/1000;
end
scale_int_size_3(scale_int_size_3<0)=0;
full_int_size_3 = [fliplr(scale_int_size_3(1,2:end)) scale_int_size_3];
full_int_size_3 = full_int_size_3/sum(full_int_size_3);

scale_int_size_4 = zeros(1,4);
for i=1:4
    scale_int_size_4(1,i) = 1-abs(300*(i-1))/1000;
end
scale_int_size_4(scale_int_size_4<0)=0;
full_int_size_4 = [fliplr(scale_int_size_4(1,2:end)) scale_int_size_4];
full_int_size_4 = full_int_size_4/sum(full_int_size_4);

scale_int_size_5 = zeros(1,5);
for i=1:5
    scale_int_size_5(1,i) = 1-abs(240*(i-1))/1000;
end
scale_int_size_5(scale_int_size_5<0)=0;
full_int_size_5 = [fliplr(scale_int_size_5(1,2:end)) scale_int_size_5];
full_int_size_5 = full_int_size_5/sum(full_int_size_5);

scale_int_size_7 = zeros(1,7);
for i=1:7
    scale_int_size_7(1,i) = 1-abs((1200/7)*(i-1))/1000;
end
scale_int_size_7(scale_int_size_7<0)=0;
full_int_size_7 = [fliplr(scale_int_size_7(1,2:end)) scale_int_size_7];
full_int_size_7 = full_int_size_7/sum(full_int_size_7);

scale_int_size_10 = zeros(1,10);
for i=1:10
    scale_int_size_10(1,i) = 1-abs(120*(i-1))/1000;
end
scale_int_size_10(scale_int_size_10<0)=0;
full_int_size_10 = [fliplr(scale_int_size_10(1,2:end)) scale_int_size_10];
full_int_size_10 = full_int_size_10/sum(full_int_size_10);

scale_int_size_11 = zeros(1,11);
for i=1:11
    scale_int_size_11(1,i) = 1-abs((1200/11)*(i-1))/1000;
end
scale_int_size_11(scale_int_size_11<0)=0;
full_int_size_11 = [fliplr(scale_int_size_11(1,2:end)) scale_int_size_11];
full_int_size_11 = full_int_size_11/sum(full_int_size_11);

scale_int_size_12 = zeros(1,12);
for i=1:12
    scale_int_size_12(1,i) = 1-abs(100*(i-1))/1000;
end
scale_int_size_12(scale_int_size_12<0)=0;
full_int_size_12 = [fliplr(scale_int_size_12(1,2:end)) scale_int_size_12];
full_int_size_12 = full_int_size_12/sum(full_int_size_12);

scale_int_size_13 = zeros(1,13);
for i=1:13
    scale_int_size_13(1,i) = 1-abs((1200/13)*(i-1))/1000;
end
scale_int_size_13(scale_int_size_13<0)=0;
full_int_size_13 = [fliplr(scale_int_size_13(1,2:end)) scale_int_size_13];
full_int_size_13 = full_int_size_13/sum(full_int_size_13);

scale_int_size_15 = zeros(1,15);
for i=1:15
    scale_int_size_15(1,i) = 1-abs(80*(i-1))/1000;
end
scale_int_size_15(scale_int_size_15<0)=0;
full_int_size_15 = [fliplr(scale_int_size_15(1,2:end)) scale_int_size_15];
full_int_size_15 = full_int_size_15/sum(full_int_size_15);

scale_int_size_16 = zeros(1,16);
for i=1:16
    scale_int_size_16(1,i) = 1-abs(75*(i-1))/1000;
end
scale_int_size_16(scale_int_size_16<0)=0;
full_int_size_16 = [fliplr(scale_int_size_16(1,2:end)) scale_int_size_16];
full_int_size_16 = full_int_size_16/sum(full_int_size_16);

scale_int_size_17 = zeros(1,17);
for i=1:17
    scale_int_size_17(1,i) = 1-abs((1200/17)*(i-1))/1000;
end
scale_int_size_17(scale_int_size_17<0)=0;
full_int_size_17 = [fliplr(scale_int_size_17(1,2:end)) scale_int_size_17];
full_int_size_17 = full_int_size_17/sum(full_int_size_17);

%% Probs of intervals (-1200, +1200) cents due to gen move and int size
full_int_prob_3 = full_int_gen_3 .* full_int_size_3;
full_int_prob_3 = full_int_prob_3/sum(full_int_prob_3);

full_int_prob_4 = full_int_gen_4 .* full_int_size_4;
full_int_prob_4 = full_int_prob_4/sum(full_int_prob_4);

full_int_prob_5 = full_int_gen_5 .* full_int_size_5;
full_int_prob_5 = full_int_prob_5/sum(full_int_prob_5);

full_int_prob_7 = full_int_gen_7 .* full_int_size_7;
full_int_prob_7 = full_int_prob_7/sum(full_int_prob_7);

full_int_prob_10 = full_int_gen_10 .* full_int_size_10;
full_int_prob_10 = full_int_prob_10/sum(full_int_prob_10);

full_int_prob_11 = full_int_gen_11 .* full_int_size_11;
full_int_prob_11 = full_int_prob_11/sum(full_int_prob_11);

full_int_prob_12 = full_int_gen_12 .* full_int_size_12;
full_int_prob_12 = full_int_prob_12/sum(full_int_prob_12);

full_int_prob_13 = full_int_gen_13 .* full_int_size_13;
full_int_prob_13 = full_int_prob_13/sum(full_int_prob_13);

full_int_prob_15 = full_int_gen_15 .* full_int_size_15;
full_int_prob_15 = full_int_prob_15/sum(full_int_prob_15);

full_int_prob_16 = full_int_gen_16 .* full_int_size_16;
full_int_prob_16 = full_int_prob_16/sum(full_int_prob_16);

full_int_prob_17 = full_int_gen_17 .* full_int_size_17;
full_int_prob_17 = full_int_prob_17/sum(full_int_prob_17);

%% Probs of intervals expressed as positive intervals
scale_int_prob_3...
    = [full_int_prob_3(ceil(length(full_int_prob_3)/2)) fliplr(full_int_prob_3(1:floor(length(full_int_prob_3)/2))) ...
    + full_int_prob_3(1,1+ceil(length(full_int_prob_3)/2):end)];

scale_int_prob_4...
    = [full_int_prob_4(ceil(length(full_int_prob_4)/2)) fliplr(full_int_prob_4(1:floor(length(full_int_prob_4)/2))) ...
    + full_int_prob_4(1,1+ceil(length(full_int_prob_4)/2):end)];

scale_int_prob_5...
    = [full_int_prob_5(ceil(length(full_int_prob_5)/2)) fliplr(full_int_prob_5(1:floor(length(full_int_prob_5)/2))) ...
    + full_int_prob_5(1,1+ceil(length(full_int_prob_5)/2):end)];

scale_int_prob_7...
    = [full_int_prob_7(ceil(length(full_int_prob_7)/2)) fliplr(full_int_prob_7(1:floor(length(full_int_prob_7)/2))) ...
    + full_int_prob_7(1,1+ceil(length(full_int_prob_7)/2):end)];

scale_int_prob_10...
    = [full_int_prob_10(ceil(length(full_int_prob_10)/2)) fliplr(full_int_prob_10(1:floor(length(full_int_prob_10)/2))) ...
    + full_int_prob_10(1,1+ceil(length(full_int_prob_10)/2):end)];

scale_int_prob_11...
    = [full_int_prob_11(ceil(length(full_int_prob_11)/2)) fliplr(full_int_prob_11(1:floor(length(full_int_prob_11)/2))) ...
    + full_int_prob_11(1,1+ceil(length(full_int_prob_11)/2):end)];

scale_int_prob_12...
    = [full_int_prob_12(ceil(length(full_int_prob_12)/2)) fliplr(full_int_prob_12(1:floor(length(full_int_prob_12)/2))) ...
    + full_int_prob_12(1,1+ceil(length(full_int_prob_12)/2):end)];

scale_int_prob_13...
    = [full_int_prob_13(ceil(length(full_int_prob_13)/2)) fliplr(full_int_prob_13(1:floor(length(full_int_prob_13)/2))) ...
    + full_int_prob_13(1,1+ceil(length(full_int_prob_13)/2):end)];

scale_int_prob_15...
    = [full_int_prob_15(ceil(length(full_int_prob_15)/2)) fliplr(full_int_prob_15(1:floor(length(full_int_prob_15)/2))) ...
    + full_int_prob_15(1,1+ceil(length(full_int_prob_15)/2):end)];

scale_int_prob_16...
    = [full_int_prob_16(ceil(length(full_int_prob_16)/2)) fliplr(full_int_prob_16(1:floor(length(full_int_prob_16)/2))) ...
    + full_int_prob_16(1,1+ceil(length(full_int_prob_16)/2):end)];

scale_int_prob_17...
    = [full_int_prob_17(ceil(length(full_int_prob_17)/2)) fliplr(full_int_prob_17(1:floor(length(full_int_prob_17)/2))) ...
    + full_int_prob_17(1,1+ceil(length(full_int_prob_17)/2):end)];

%% Expected spectral pitch similarity (all spectral tunings all 
% underlying tunings).
if sigma_S_chg==1 || rolloff_S==1
    expected_spsim_3 ...
        = sum(bsxfun(@times,scale_int_prob_3,1-spec_dist_3),2);
    expected_spsim_4 ...
        = sum(bsxfun(@times,scale_int_prob_4,1-spec_dist_4),2);
    expected_spsim_5 ...
        = sum(bsxfun(@times,scale_int_prob_5,1-spec_dist_5),2);
    expected_spsim_7 ...
        = sum(bsxfun(@times,scale_int_prob_7,1-spec_dist_7),2);
    expected_spsim_10 ...
        = sum(bsxfun(@times,scale_int_prob_10,1-spec_dist_10),2);
    expected_spsim_11 ...
        = sum(bsxfun(@times,scale_int_prob_11,1-spec_dist_11),2);
    expected_spsim_12 ...
        = sum(bsxfun(@times,scale_int_prob_12,1-spec_dist_12),2);
    expected_spsim_13 ...
        = sum(bsxfun(@times,scale_int_prob_13,1-spec_dist_13),2);
    expected_spsim_15 ...
        = sum(bsxfun(@times,scale_int_prob_15,1-spec_dist_15),2);
    expected_spsim_16 ...
        = sum(bsxfun(@times,scale_int_prob_16,1-spec_dist_16),2);
    expected_spsim_17 ...
        = sum(bsxfun(@times,scale_int_prob_17,1-spec_dist_17),2);
    
    expected_spsim_all ...
        = [expected_spsim_3  expected_spsim_4  expected_spsim_5 ...
        expected_spsim_7  expected_spsim_10 expected_spsim_11 ...
        expected_spsim_12 expected_spsim_13 expected_spsim_15 ...
        expected_spsim_16 expected_spsim_17]';
end

%% Calculate the virtual pitch distance for all upwards intervals 
% (columns) in each tuning for each spectral tuning (rows)

% The faster code for circular shift is taken Matt J:
% http://www.mathworks.com/matlabcentral/newsreader/view_thread/279310
virt_dist_3 = zeros(11,3);
step_size = 400;
steps = round(0:step_size:2*step_size);
for j=1:11
    virtp_tun = virtp(j,:);
    virtp_tun = virtp_tun(ones(3,1),:)';
    col = virtp_tun(:,1);
    for i=1:3
        step_i = steps(i);
        col(:) = virtp_tun(:,i);
        virtp_tun(1:end-step_i,i) = col(step_i+1:end);
        virtp_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    virtp_shft = virtp_tun';
    virtp_shft_1 = virtp_shft(1,:);
    for i=1:3
        virtp_shft_i = virtp_shft(i,:);
        virt_dist_3(j,i)...
            = 1-(virtp_shft_1*virtp_shft_i' ...
            / sqrt((virtp_shft_1*virtp_shft_1')*(virtp_shft_i*virtp_shft_i')));
    end
end
virt_dist_3;

virt_dist_4 = zeros(11,4);
step_size = 300;
steps = round(0:step_size:3*step_size);
for j=1:11
    virtp_tun = virtp(j,:);
    virtp_tun = virtp_tun(ones(4,1),:)';
    col = virtp_tun(:,1);
    for i=1:4
        step_i = steps(i);
        col(:) = virtp_tun(:,i);
        virtp_tun(1:end-step_i,i) = col(step_i+1:end);
        virtp_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    virtp_shft = virtp_tun';
    virtp_shft_1 = virtp_shft(1,:);
    for i=1:4
        virtp_shft_i = virtp_shft(i,:);
        virt_dist_4(j,i)...
            = 1-(virtp_shft_1*virtp_shft_i' ...
            / sqrt((virtp_shft_1*virtp_shft_1')*(virtp_shft_i*virtp_shft_i')));
    end
end
virt_dist_4;

virt_dist_5 = zeros(11,5);
step_size = 240;
steps = round(0:step_size:4*step_size);
for j=1:11
    virtp_tun = virtp(j,:);
    virtp_tun = virtp_tun(ones(5,1),:)';
    col = virtp_tun(:,1);
    for i=1:5
        step_i = steps(i);
        col(:) = virtp_tun(:,i);
        virtp_tun(1:end-step_i,i) = col(step_i+1:end);
        virtp_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    virtp_shft = virtp_tun';
    virtp_shft_1 = virtp_shft(1,:);
    for i=1:5
        virtp_shft_i = virtp_shft(i,:);
        virt_dist_5(j,i)...
            = 1-(virtp_shft_1*virtp_shft_i' ...
            / sqrt((virtp_shft_1*virtp_shft_1')*(virtp_shft_i*virtp_shft_i')));
    end
end
virt_dist_5;

virt_dist_7 = zeros(11,7);
step_size = 1200/7;
steps = round(0:step_size:6*step_size);
for j=1:11
    virtp_tun = virtp(j,:);
    virtp_tun = virtp_tun(ones(7,1),:)';
    col = virtp_tun(:,1);
    for i=1:7
        step_i = steps(i);
        col(:) = virtp_tun(:,i);
        virtp_tun(1:end-step_i,i) = col(step_i+1:end);
        virtp_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    virtp_shft = virtp_tun';
    virtp_shft_1 = virtp_shft(1,:);
    for i=1:7
        virtp_shft_i = virtp_shft(i,:);
        virt_dist_7(j,i)...
            = 1-(virtp_shft_1*virtp_shft_i' ...
            / sqrt((virtp_shft_1*virtp_shft_1')*(virtp_shft_i*virtp_shft_i')));
    end
end
virt_dist_7;

virt_dist_10 = zeros(11,10);
step_size = 120;
steps = round(0:step_size:9*step_size);
for j=1:11
    virtp_tun = virtp(j,:);
    virtp_tun = virtp_tun(ones(10,1),:)';
    col = virtp_tun(:,1);
    for i=1:10
        step_i = steps(i);
        col(:) = virtp_tun(:,i);
        virtp_tun(1:end-step_i,i) = col(step_i+1:end);
        virtp_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    virtp_shft = virtp_tun';
    virtp_shft_1 = virtp_shft(1,:);
    for i=1:10
        virtp_shft_i = virtp_shft(i,:);
        virt_dist_10(j,i)...
            = 1-(virtp_shft_1*virtp_shft_i' ...
            / sqrt((virtp_shft_1*virtp_shft_1')*(virtp_shft_i*virtp_shft_i')));
    end
end
virt_dist_10;

virt_dist_11 = zeros(11,11);
step_size = 1200/11;
steps = round(0:step_size:10*step_size);
for j=1:11
    virtp_tun = virtp(j,:);
    virtp_tun = virtp_tun(ones(11,1),:)';
    col = virtp_tun(:,1);
    for i=1:11
        step_i = steps(i);
        col(:) = virtp_tun(:,i);
        virtp_tun(1:end-step_i,i) = col(step_i+1:end);
        virtp_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    virtp_shft = virtp_tun';
    virtp_shft_1 = virtp_shft(1,:);
    for i=1:11
        virtp_shft_i = virtp_shft(i,:);
        virt_dist_11(j,i)...
            = 1-(virtp_shft_1*virtp_shft_i' ...
            / sqrt((virtp_shft_1*virtp_shft_1')*(virtp_shft_i*virtp_shft_i')));
    end
end
virt_dist_11;

virt_dist_12 = zeros(11,12);
step_size = 100;
steps = round(0:step_size:11*step_size);
for j=1:11
    virtp_tun = virtp(j,:);
    virtp_tun = virtp_tun(ones(12,1),:)';
    col = virtp_tun(:,1);
    for i=1:12
        step_i = steps(i);
        col(:) = virtp_tun(:,i);
        virtp_tun(1:end-step_i,i) = col(step_i+1:end);
        virtp_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    virtp_shft = virtp_tun';
    virtp_shft_1 = virtp_shft(1,:);
    for i=1:12
        virtp_shft_i = virtp_shft(i,:);
        virt_dist_12(j,i)...
            = 1-(virtp_shft_1*virtp_shft_i' ...
            / sqrt((virtp_shft_1*virtp_shft_1')*(virtp_shft_i*virtp_shft_i')));
    end
end
virt_dist_12;

virt_dist_13 = zeros(11,13);
step_size = 1200/13;
steps = round(0:step_size:12*step_size);
for j=1:11
    virtp_tun = virtp(j,:);
    virtp_tun = virtp_tun(ones(13,1),:)';
    col = virtp_tun(:,1);
    for i=1:13
        step_i = steps(i);
        col(:) = virtp_tun(:,i);
        virtp_tun(1:end-step_i,i) = col(step_i+1:end);
        virtp_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    virtp_shft = virtp_tun';
    virtp_shft_1 = virtp_shft(1,:);
    for i=1:13
        virtp_shft_i = virtp_shft(i,:);
        virt_dist_13(j,i)...
            = 1-(virtp_shft_1*virtp_shft_i' ...
            / sqrt((virtp_shft_1*virtp_shft_1')*(virtp_shft_i*virtp_shft_i')));
    end
end
virt_dist_13;

virt_dist_15 = zeros(11,15);
step_size = 80;
steps = round(0:step_size:14*step_size);
for j=1:11
    virtp_tun = virtp(j,:);
    virtp_tun = virtp_tun(ones(15,1),:)';
    col = virtp_tun(:,1);
    for i=1:15
        step_i = steps(i);
        col(:) = virtp_tun(:,i);
        virtp_tun(1:end-step_i,i) = col(step_i+1:end);
        virtp_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    virtp_shft = virtp_tun';
    virtp_shft_1 = virtp_shft(1,:);
    for i=1:15
        virtp_shft_i = virtp_shft(i,:);
        virt_dist_15(j,i)...
            = 1-(virtp_shft_1*virtp_shft_i' ...
            / sqrt((virtp_shft_1*virtp_shft_1')*(virtp_shft_i*virtp_shft_i')));
    end
end
virt_dist_15;

virt_dist_16 = zeros(11,16);
step_size = 75;
steps = round(0:step_size:15*step_size);
for j=1:11
    virtp_tun = virtp(j,:);
    virtp_tun = virtp_tun(ones(16,1),:)';
    col = virtp_tun(:,1);
    for i=1:16
        step_i = steps(i);
        col(:) = virtp_tun(:,i);
        virtp_tun(1:end-step_i,i) = col(step_i+1:end);
        virtp_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    virtp_shft = virtp_tun';
    virtp_shft_1 = virtp_shft(1,:);
    for i=1:16
        virtp_shft_i = virtp_shft(i,:);
        virt_dist_16(j,i)...
            = 1-(virtp_shft_1*virtp_shft_i' ...
            / sqrt((virtp_shft_1*virtp_shft_1')*(virtp_shft_i*virtp_shft_i')));
    end
end
virt_dist_16;

virt_dist_17 = zeros(11,17);
step_size = 1200/17;
steps = round(0:step_size:16*step_size);
for j=1:11
    virtp_tun = virtp(j,:);
    virtp_tun = virtp_tun(ones(17,1),:)';
    col = virtp_tun(:,1);
    for i=1:17
        step_i = steps(i);
        col(:) = virtp_tun(:,i);
        virtp_tun(1:end-step_i,i) = col(step_i+1:end);
        virtp_tun(end+1-step_i:end,i) = col(1:step_i);
    end
    virtp_shft = virtp_tun';
    virtp_shft_1 = virtp_shft(1,:);
    for i=1:17
        virtp_shft_i = virtp_shft(i,:);
        virt_dist_17(j,i)...
            = 1-(virtp_shft_1*virtp_shft_i' ...
            / sqrt((virtp_shft_1*virtp_shft_1')*(virtp_shft_i*virtp_shft_i')));
    end
end
virt_dist_17;

%% Expected virtual pitch similarity (all spectral tunings all 
% underlying tunings).
if sigma_S_chg==1 || rolloff_S==1
    expected_vpsim_3 ...
        = sum(bsxfun(@times,scale_int_prob_3,1-virt_dist_3),2);
    expected_vpsim_4 ...
        = sum(bsxfun(@times,scale_int_prob_4,1-virt_dist_4),2);
    expected_vpsim_5 ...
        = sum(bsxfun(@times,scale_int_prob_5,1-virt_dist_5),2);
    expected_vpsim_7 ...
        = sum(bsxfun(@times,scale_int_prob_7,1-virt_dist_7),2);
    expected_vpsim_10 ...
        = sum(bsxfun(@times,scale_int_prob_10,1-virt_dist_10),2);
    expected_vpsim_11 ...
        = sum(bsxfun(@times,scale_int_prob_11,1-virt_dist_11),2);
    expected_vpsim_12 ...
        = sum(bsxfun(@times,scale_int_prob_12,1-virt_dist_12),2);
    expected_vpsim_13 ...
        = sum(bsxfun(@times,scale_int_prob_13,1-virt_dist_13),2);
    expected_vpsim_15 ...
        = sum(bsxfun(@times,scale_int_prob_15,1-virt_dist_15),2);
    expected_vpsim_16 ...
        = sum(bsxfun(@times,scale_int_prob_16,1-virt_dist_16),2);
    expected_vpsim_17 ...
        = sum(bsxfun(@times,scale_int_prob_17,1-virt_dist_17),2);
    
    expected_vpsim_all ...
        = [expected_vpsim_3  expected_vpsim_4  expected_vpsim_5 ...
        expected_vpsim_7  expected_vpsim_10 expected_vpsim_11 ...
        expected_vpsim_12 expected_vpsim_13 expected_vpsim_15 ...
        expected_vpsim_16 expected_vpsim_17]';
end

end

%% Put variables and data into vectors (including or excluding 
% self-comparisons) and standardizing or not (ind. variables only).
spsim_matr = zeros(11,11);
for i=1:11
    for j=1:11
        spsim_matr(i,j) = expected_spsim_all(i,i) - expected_spsim_all(i,j);
    end
end
% spsim_matr = log((spsim_matr+1)./2) - log(1-((spsim_matr+1)./2)); this
% transform makes no sinificant difference to final fit

%%%%%%%%%%%%%%%%%%%%%%
% spsim_self is vector of spectral similarities including self-comparisons
spsim_self = spsim_matr(:);
%%%%%%%%%%%%%%%%%%%%%%

spsim = spsim_self;
spsim(1:12:121) = [];
%%%%%%%%%%%%%%%%%%%%%%%
% spsim is vector of spectral similarities without self-comparisons
spsim;
%%%%%%%%%%%%%%%%%%%%%%%

spsim(valid_stims) = []; % This is the spsim for the training set (when doing 
% cross-validation)


vpsim_matr = zeros(11,11);
for i=1:11
    for j=1:11
        vpsim_matr(i,j) = expected_vpsim_all(i,i) - expected_vpsim_all(i,j);
    end
end
% vpsim_matr = log((vpsim_matr+1)./2) - log(1-((vpsim_matr+1)./2)); this
% transform makes no sinificant difference to final fit

%%%%%%%%%%%%%%%%%%%%%%
% vpsim_self is vector of spectral similarities including self-comparisons
vpsim_self = vpsim_matr(:);
%%%%%%%%%%%%%%%%%%%%%%

vpsim = vpsim_self;
vpsim(1:12:121) = [];
%%%%%%%%%%%%%%%%%%%%%%%
% vpsim is vector of spectral similarities without self-comparisons
vpsim;
%%%%%%%%%%%%%%%%%%%%%%%

vpsim(valid_stims) = []; % This is the vpsim for the training set (when doing 
% cross-validation)


harm_matr = zeros(11,11);
for i=1:11
    for j=1:11
        harm_matr(i,j) = maxes(i) - maxes(j);
    end
end
% harm_matr = log((harm_matr+1)./2) - log(1-((harm_matr+1)./2)); this
% transform makes no sinificant difference to final fit

%%%%%%%%%%%%%%%%%%%%%%
% harm_self is vector of harmonicities including self-comparisons
harm_self = harm_matr(:);
%%%%%%%%%%%%%%%%%%%%%%

harm = harm_self;
harm(1:12:121) = [];
%%%%%%%%%%%%%%%%%%%%%%%
% vpsim is vector of harmonicities without self-comparisons
harm;
%%%%%%%%%%%%%%%%%%%%%%%

harm(valid_stims) = []; % This is the harm for the training set (when doing 
% cross-validation)

%% Logistic values and likelihood calculation
if constant==1
    if variables==0
        yfit ...
        = glmval([b_const b_spsim b_harm]',[spsim harm],'logit','size',trials_train);
    elseif variables==1
        yfit ...
        = glmval([b_const b_vpsim b_harm]',[vpsim harm],'logit','size',trials_train);
    elseif variables==2
        yfit ...
        = glmval([b_const b_spsim b_vpsim b_harm]',[spsim vpsim harm],'logit','size',trials_train);
    end        
else
    if variables==0
        yfit ...
        = glmval([b_spsim b_harm]',[spsim harm],'logit','size',trials_train,'constant','off');
    elseif variables==1        
        yfit ...
        = glmval([b_vpsim b_harm]',[vpsim harm],'logit','size',trials_train,'constant','off');
    elseif variables==2
        yfit ...
        = glmval([b_spsim b_vpsim b_harm]',[spsim vpsim harm],'logit','size',trials_train,'constant','off');
    end
end
Deviance = -2*sum(log(binopdf(successes_train,trials_train,yfit./trials_train)));

if variables==0
    beta_hat = [b_const b_spsim b_harm]';
elseif variables==1
    beta_hat = [b_const b_vpsim b_harm]';
elseif variables==2
    beta_hat = [b_const b_spsim b_vpsim b_harm]';    
end

