% Generate simulated dataset to try out the toolbox
close all; clear all; clc;  % clear the memory etc

% simulate dataset with varying standards noise levels and the whole shabang

nrSubs = 10; % simulate 10 participants

standards = [0];
nrStandards = length(standards);

mu = 0;                                 % mean (PSE) is zero

sigCue1 = 15;                           % simulate one noise level for cue 1
sigCue2 = ((0.5:0.5:2.0)-0.25)*sigCue1; % simulate multiple noise levels for cue 2
nrNL = length(sigCue2);

conflicts = [-20,-10,10,20];            % simulate different conflicts between cues

fitlam = 0;                             % set not to assume lapse rate = 0
sampleRange = 30*[-1,1];

sigA = 3;                               % size of subject dependent noise

% subject dependend variations in terms of sigma (note base is not yet added)
ransubsiglst = rand(nrSubs,2);
ransubsiglst(:,1) = ransubsiglst(:,1) * 2 * sigA - sigA;
ransubsiglst(:,2) = ransubsiglst(:,2) * 2 * sigA - sigA;


% -------------------------------------------------------------------------
%   simulate constant stimuli
% -------------------------------------------------------------------------
nrcomp = 13;
nrreps = 15;
stepS = diff(sampleRange)/(nrcomp-1); % stepsize between comparisons
stimlst = sampleRange(1)+(0:(nrcomp-1))*stepS;
stimlst = stimlst' * ones(1,nrreps);
stimlst = stimlst(:);
nrsamples = length(stimlst);
allsubdata = [];
for j = 1:nrSubs,
    subdata = [];
    
    mu = 0;
    % simulate unimodel cue 1
    sigma = sigCue1 + ransubsiglst(j,1);
    proplst = normcdf(stimlst,mu,sigma);

    for s = 1:nrStandards
        responses = (rand(nrsamples,1) < proplst);
        subdata = [subdata;...
            [stimlst,responses,ones(nrsamples,1),...    % stimulus response lst (stimlst(:,1) is comparisn - standard)
             standards(s)*ones(nrsamples,1),...         % standard stimulus
             1*ones(nrsamples,1),0*ones(nrsamples,1),...% cues
             0*ones(nrsamples,1),0*ones(nrsamples,1),...% added conflict (relative to standard)
             0*ones(nrsamples,1),0*ones(nrsamples,1)... % noise levels
             ]];
    end

    % simulate unimodel cue 2
    for n = 1:nrNL,
        sigma = sigCue2(n) + ransubsiglst(j,2);
        proplst = normcdf(stimlst,mu,sigma);

        for s = 1:nrStandards
            responses = (rand(nrsamples,1) < proplst);
            subdata = [subdata;...
                [stimlst,responses,ones(nrsamples,1),...    % stimulus response lst (stimlst(:,1) is comparisn - standard)
                 standards(s)*ones(nrsamples,1),...         % standard stimulus
                 0*ones(nrsamples,1),1*ones(nrsamples,1),...% cues
                 0*ones(nrsamples,1),0*ones(nrsamples,1),...% added conflict (relative to standard)
                 0*ones(nrsamples,1),sigCue2(n)*ones(nrsamples,1)... % noise levels
                 ]];
        end
    end    

    
    % simulate bimodal
    for n = 1:nrNL,
        for c = 1:length(conflicts)
            ma = conflicts(c)/2;
            mu1 = mu + ma;
            mu2 = mu - ma;
            sigma1 = sigCue1    + ransubsiglst(j,1);
            sigma2 = sigCue2(n) + ransubsiglst(j,2);
            %
            mubi = ((mu1/sigma1^2)+(mu2/sigma2^2))/((1/sigma1^2)+(1/sigma2^2));
            sigmabi = sqrt(sigma1^2*sigma2^2/(sigma1^2+sigma2^2));
            proplst = normcdf(stimlst,mubi,sigmabi);

            for s = 1:nrStandards
                responses = (rand(nrsamples,1) < proplst);
                subdata = [subdata;...
                    [stimlst,responses,ones(nrsamples,1),...    % stimulus response lst (stimlst(:,1) is comparisn - standard)
                     standards(s)*ones(nrsamples,1),...         % standard stimulus
                      1*ones(nrsamples,1),  1*ones(nrsamples,1),...% cues
                     ma*ones(nrsamples,1),-ma*ones(nrsamples,1),...% added conflict (relative to standard)
                      0*ones(nrsamples,1), sigCue2(n)*ones(nrsamples,1)... % noise levels
                     ]];
            end
            
        end
    end    

    % analyse each individual similated participant
%    testres = MLE_SingleSubject(subdata,2,0,[1,1])
    
    % attach subject number in last column and add to overalsubmatrix
    allsubdata = [allsubdata;[subdata,j*ones(length(subdata(:,1)),1)]];
    
end

% analyse across simulated participants
[results,resultsSub] = MLE_Analysis(allsubdata,2,0,[],[0,0,1]);