%% This function uses spectral pitch similarity to model Krumhansl's probe
% tone data. It measures the spectral pitch distance between all chromatic
% pitch classes and a major and minor triad (with different
% root-weightings). The deviance is calculated as the sum of squared errors.

% It based on routines written for Milne, A.J., Sethares, W.A., 
% Laney, R., Sharp, D.B. (2011) Modelling the Similarity of Pitch
% Collections with Expectation Tensors, Journal of Mathematics and Music,
% 5(1), 1-20.

% Author: Andrew J. Milne
% Version: 1.1

function Deviance = Probe_tone_spectral_sim_model(param_vals)

global spec_sim weight_non_roots weight_m3 K_data_train valid_ind ...
    num_harmonics shift_by virtual len_shifts exp_or_sal

rolloff = param_vals(1);
sigma = param_vals(2);
omega = param_vals(3);
intercept = param_vals(4);
slope = param_vals(5);

%% Parameters
% "period" is the size in cents of the octave (or other period). The
% standard octave is 1200 cents. All entered pitch values are reduced to 
% lie in the interval [0, period).
% "q" = the resolution (granularity) of the pitch embedding (in divisions 
% of the period). A higher value (e.g., 1200) allows for more accurate
% modelling of perceptual uncertainty, and the distances between closely 
% spaced pitches, but increases calculation and memory requirements. If 
% modelling standard 12-EDO pitch classes, a value of 12 is adequate.
% "metric" = 'l_p' or 'cosine' 
% "p" = size of the norm when using l_p metrics.
% "sigma" = standard deviation, in cents, of the Gaussian probability
% function
% "win_len" = the length of the smoothing kernel (in standard deviations)
if virtual==1
    period = 6000;
    q = 6000;
else
    period = 1200;
    q = 1200;
end
metric = 'cosine';
p = 1;
win_len = 12;

%% Reference and generated triad parameters
% "low_limit" and "high_limit" sets the range (in cents) of the continuum
% triad relative to the reference triad 
% "inc" sets the size, in cents, of the increments between low_limit and 
% high_limit 
% "num_harmonics" sets the number of spectral pitches in each tone 
% "rolloff" sets the amplitude of each tone's harmonics to 1/n^rolloff, 
% where n is the harmonic number; rolloff = 0 gives a flat spectrum (all 
% saliences are 1; rolloff = 1 gives a sawtooth spectrum; rolloff = a very 
% big number results in all harmonics, except the fundamental, having zero 
% salience (though this is easier to achieve by setting num_harmonics to 1).
% "RefTriad" - can be either 'Major' or 'Minor' 
% "Gamma" - the amount of tempering of the reference triad's tones away 
% from just intonation (0 is tempered, 1 is just)
% "TuningContinuum" = 'Syntonic', 'Magic', or 'Hanson' - various ways of 
% generating the reference triad from two generators alpha and beta. To set
% the size of alpha and beta look inside the if loops in the Generate
% reference triad section.
low_limit = 0;
high_limit = 1200;
inc = 100;

shifts = low_limit:inc:high_limit-inc;
len_shifts = length(shifts);
Dist = zeros(1,len_shifts);

%% Generate reference triad's spectrum
SpectrumX = zeros(1, num_harmonics);
for j = 1:num_harmonics
    SpectrumX(j) = 1200*log2(j);
end

SpectrumX_w = zeros(1, num_harmonics);
for j = 1:num_harmonics
    SpectrumX_w(j) = 1/(j^rolloff);
end

SpectrumY_base = SpectrumX;
SpectrumY_w = SpectrumX_w;

if weight_non_roots==0 && weight_m3==0
    Ref_chord_wts = [1 1 1; 1 1 1];
elseif weight_non_roots==1 && weight_m3==0
    Ref_chord_wts = [1 omega omega; 1 omega omega];
elseif weight_non_roots==1 && weight_m3==1
    Ref_chord_wts = [1 omega omega; 1 1     omega];
end

for ref_mode = 0:1
    SpectrumX_chord = [SpectrumX 400-100*ref_mode+SpectrumX 700+SpectrumX];
    SpectrumX_w_chord ...
        = [Ref_chord_wts(ref_mode+1,1)*SpectrumX_w ...
           Ref_chord_wts(ref_mode+1,2)*SpectrumX_w ...
           Ref_chord_wts(ref_mode+1,3)*SpectrumX_w];
    
    SpectrumX_chord = SpectrumX_chord(:)';
    SpectrumX_w_chord = SpectrumX_w_chord(:)';
    
    %% Generate continuum tones and send to embeddings
    for i = 1:len_shifts
        shift_by = shifts(i);
        SpectrumY = SpectrumY_base + shift_by;
        Dist(i+ref_mode*len_shifts) ...
            = ExpectationMetrics(period,q,metric,p,sigma,win_len,...
            SpectrumX_chord,SpectrumX_w_chord,SpectrumY,SpectrumY_w,'abs',1,exp_or_sal);
    end
end

spec_sim = (1-Dist)';
spec_sim_train = spec_sim;
spec_sim_train(valid_ind) = [];

%% final model
model_data = intercept + slope.*spec_sim_train;

%% Calculate Deviance
Deviance = sum((model_data-K_data_train).^2);

