% Mambila downdrift - stochastic
%
% Sarah Harper's modification of Louis Goldstein's original script
% April 2 2020

% Model the results in Table 1 of Connell's article with stochastic
% variation in starting parameter values. Code set up to generate and save
% data from N runs of the model, so that generalizations about behavior of system 
% in the presence of variation/conclusions about the effect of variation in
% parameter values on the system can be evaluated.

% Idea:
% Over time, the internal target F0 of all four tones show decrease
% in F0, under the control of a 1st order linear declination system
% for each tone.
%
% The tone targets also repel each other with a force that is inversely
% propotional to the current f0 difference beteween the two tones.
%
% The precise strength of the repulsion and the strength of the declination
% force are randomized within a controlled range in this script. Initial
% targets for syllable 1 are set by random sampling from a gaussian
% distribution defined based on data from 
% http://www.mambila.info/Connell/downtrend/downtrend_abstract.html
%
% INPUT to script is the number N of simulations of the model that should be
% run (I've generally set N to 100)
%
% First OUTPUT of script is a struct with N rows (each row  = 1 run of model)
% and fields:
%   (1) declination = value of declination force in run
%   (2) repulsion = value of repulsive force in run
%   (3) ratio = ratio of declination/repulsion
%   (4) syll = number of syllables (default = 10; held constant across all
%   runs, but include as a field to facilitate plotting)
%   (5)-(8) = Tone targets simulated across 10 syllables (1st value =
%   randomly generated starting target value for that tone in that run)
%   (5 = T1, 6 = T2, 7 = T3, 8 = T4) 
%   (9)-(12) = Slope of tone trajectory for each tone in run
%   (9 = T1, 10 = T2, 11 = T3, 12 = T4)
%
% Second OUTPUT is a 2x4 matrix containing the regression coefficients (1st
% column) and intercepts (2nd column) for the relationship between
% declination/repulsion ratio and slope across all runs in simulation for
% each tone (T1 = row 1, T2 = row 2, T3 = row 3, T4 = row 4)

function [rundata, overall_corr] = Mambila_stochastic_simulations(N)

%Setting up simulation parameters (N = number of runs in simulation)
ntones = 4;

for n = 1:N
% Set the initial (syll 1) f0 value of each tone
% Value for each tone generated randomly from a normal distribution with
% mean and SD based on Connell data (note: SD used for each tone is half of
% SD in the data)
% Settings based on speaker JE:
    T(1,1) = normrnd(142.17, 4.4);
    T(2,1) = normrnd(135.72, 5.6);
    while T(2,1) >= T(1,1) - 2
        T(2,1) = normrnd(135.72, 5.6);
    end
    T(3,1) = normrnd(126.55, 2.4);
    while T(3,1) >= T(2,1) - 2
      T(3,1) = normrnd(126.55, 2.4);
    end
    T(4,1) = normrnd(118.63, 1.55);
    while T(4,1) >= T(3,1) - 2
        T(4,1) = normrnd(118.63, 1.55);
    end 
% Alternative settings based on speaker NM (can be used for 
% comparison/verification of generalizability):
%     T(1,1) = normrnd(160.19, 6.35);
%     T(2,1) = normrnd(151.13, 3.1);
%     while T(2,1) >= T(1,1) - 2
%         T(2,1) = normrnd(151.13, 3.1);
%     end
%     T(3,1) = normrnd(133.53, 4.2);
%     while T(3,1) >= T(2,1) - 2
%       T(3,1) = normrnd(133.53, 4.2);
%     end
%     T(4,1) = normrnd(124.5, 4.15);
%     while T(4,1) >= T(3,1) - 2
%         T(4,1) = normrnd(124.5, 4.15);
%     end 
%Set the goal value for the declination system for each tone
%Declination value (d) randomized across runs, but held constant for all tones
%within a run
    d = randi([5,15]); 
    for j = 1:ntones
        Tgoal(j) = T(j,1) - d; 
    end


    F = zeros(ntones, ntones); % F is the matrix of inter-tone repulsive forces

    k  = .2; % stiffness of the declination systems; increased from Louis's code b/c T2 was acting weird
    r = randi([5,20]); % scaling of the repulsive forces; like declination goal value, randomized across runs

    nsyll = 10;

    for i = 2:nsyll 
         for j = 1:ntones 
             for l = 1:ntones
                if ~(j==l)
                    F(j,l) = r./(T(j,i-1) - T(l,i-1)); % repulsive forces
                end
             end
                push = sum(F(j,:)); %sum(F,1); % push total amount of force on a given tone from all the others
                dT = -k*T(j, i-1) + k*Tgoal(j) + push; %(j);  % dynamics       
                T(j,i) = T(j, i-1) + dT;
         end 
    end
    %Add all data from this run to struct containing results of all runs in
    %simulation
    rundata(n).declination = d;
    rundata(n).repulsion = r;
    rundata(n).ratio = d/r;
    rundata(n).syll = 1:nsyll;
    rundata(n).T1 = T(1,:);
    rundata(n).T2 = T(2,:);
    rundata(n).T3 = T(3,:);
    rundata(n).T4 = T(4,:);
    b1 = polyfit(rundata(n).syll,rundata(n).T1,1);
    rundata(n).b1 = b1(1);
    b2 = polyfit(rundata(n).syll,rundata(n).T2,1);
    rundata(n).b2 = b2(1);
    b3 = polyfit(rundata(n).syll,rundata(n).T3,1);
    rundata(n).b3 = b3(1);
    b4 = polyfit(rundata(n).syll,rundata(n).T4,1);
    rundata(n).b4 = b4(1);
    
end

%find slope maxima and minima to set ylims in plots
ymx = max([rundata.b1]);
ymn = min([rundata.b4]);

%calculate regression coefficient and intercept for b ~ ratio across all
%runs and create matrix with data
tab = struct2table(rundata); %trying to get corr function to work with struct fields did not work
rho1 = corr(tab.ratio, tab.b1, 'Type', 'Spearman');
rho2 = corr(tab.ratio, tab.b2, 'Type', 'Spearman');
rho3 = corr(tab.ratio, tab.b3, 'Type', 'Spearman');
rho4 = corr(tab.ratio, tab.b4, 'Type', 'Spearman');
overall_corr = [rho1; rho2; rho3; rho4];

%plot slope data
tiledlayout(2,4);

nexttile([1 4])
scatter([rundata.ratio], [rundata.b1]);
title('Slope vs. declination/repulsion ratio')
ylim([ymn ymx])
ylabel('Slope (F0 Change)')
xlabel('Declination/Repulsion')
hold on
scatter([rundata.ratio], [rundata.b2]);
scatter([rundata.ratio], [rundata.b3]);
scatter([rundata.ratio], [rundata.b4]);
hold off

nexttile
scatter([rundata.ratio], [rundata.b1]);
title(['Tone 1 (rho = ' num2str(rho1(1)) ')'])
ylim([ymn ymx])
ylabel('Slope (F0 Change)')
xlabel('Declination/Repulsion')
nexttile
scatter([rundata.ratio], [rundata.b2]);
title(['Tone 2 (rho = ' num2str(rho2(1)) ')'])
ylabel('Slope (F0 Change)')
xlabel('Declination/Repulsion')
ylim([ymn ymx])
nexttile
scatter([rundata.ratio], [rundata.b3]);
title(['Tone 3 (rho = ' num2str(rho3(1)) ')'])
ylabel('Slope (F0 Change)')
xlabel('Declination/Repulsion')
ylim([ymn ymx])
nexttile
scatter([rundata.ratio], [rundata.b4]);
title(['Tone 4 (rho = ' num2str(rho4(1)) ')'])
ylabel('Slope (F0 Change)')
xlabel('Declination/Repulsion')
ylim([ymn ymx])