function [results,flagLargeJND] = MLE_SingleSubject(A, nrStim, fitlapses, comb_stand, plotopts, verbose)

% -------------------------------------------------------------------------
%   [results] = MLE_SingleSubject(A, nrStim, fitlapses, comb_stand, ploton, verbose)
%
%   Performs the MLE analysis on the data for a single participant. This is
%   done by separating single cue data from multicue data and then fitting 
%   a cumulative Gaussian (see FitCumulativeGauss) to
%   each part.
%
%   Input parameters:
%   A = The general data matrix. 
%       Rows denote data for different comparison stimuli.
%       Each row is of the following format:
%       [C-S, P(C>S), N, S, [cues on/off], [Conf], [noise level]]  
%       where  C: comparison
%              S: standard
%              P(): proportion
%              N: number of data points (repetitions)
%              Conf: inter-sensory conflict. 
%               
%       For instance, for a two-cue integration data set, a row like this:
%       [  -50, 0.2, 10, 0, 1, 1, -10, 10, 0, 2]
%
%       indicates a difference standard/comparison of -50 (1st col), 
%       20% 'yes' responses (2nd col) for this comparison 
%       in ten repetitions (3rd col), a standard value of 0 (4th col), 
%       a bimodal presentation (5th and 6th col are both set to 1),
%       an intersensory conflict of 20: 10 subtracted from cue 1 (7th col)
%       and 10 added to cue 2 (8th col), no noise on cue 1 (9th col) and 
%       the 2nd noise level on cue 2 (10th col)
%
%       A row like this:
%       [  20, 0.9, 10, 10, 0, 1, 0, 0, 0, 0]
%
%       indicates a difference standard/comparison of 20 (1st col), 
%       90% 'yes' responses (2nd col) for this comparison 
%       in ten repetitions (3rd col), a standard value of 10 (4th col), 
%       a unimodal presentation of cue 2 (5th col is 0, 6th col is 1),
%       with no conflict and no noise added.
% 
%   nrStim = the number of stimuli presented in one trial. For instance if
%       2 interval forced choice (2IFC) was used the JND represents sqrt(2)
%       times the sigma of a single stimulus. To take this correction into
%       account, the system needs to know the number of stimuli used in one
%       trial. I.e. nrStim = 1 for task that use only one stimulus (e.g.
%       for slant this could be which side of the single stimulus is more
%       in front). nrStim = 2 for tasks like the 2IFC in which a comparison
%       stimulus is compared to a standard (which one is more slanted).
%       Since most MLE studies involve the 2AFC/2IFC type of task the
%       default is nrStim = 2 (thus taking sigma = JND/sqrt(2)). For a
%       refence see Rohde, van Dam & Ernst (hopefully 2015).
%
%   fitlapses =  Whether or not to fit an additional parameter, the 
%       lapserate, to each curve. 1 = yes; 0 = no; default is 0.
%
%   comb_stand =  whether or not to force combining the standards in the
%       in the analyis. 1 = yes; 0 is no. If data for different standards 
%       are not forcefully combined the analysis is performed for each 
%       standard separately and a test is performed to check for 
%       differences between standards. comb_stand = 1 skips this test and
%       forces the combination.
%
%   plotopts = a list of two values. The first values indicates whether or 
%       not separate curves are being plotted (0 = off; 1 = on). The second
%       whether a summary of the results is plotted. default = [0,0]
%
%   verbose = whether or not to display messages on the evaluation
%       progress. 1 = on; 0 = off. default is 0.
%
%
%   The output is a structure results with the following fields:
%   singleCues = a substructure in which each field refers to one of the
%       cues used. For the structure of this substruct see below.
%   cueCombos = a substructure in which each field refers to one of the cue
%       combination conditions. See below.
%
%
%   Structure subfields for singleCues:
%   params:  stores the fitted mu,sigma and lapserate. If there are
%       multiple levels of noise in the design each row corresponds to a
%       different noise level.
%   params_STD: standard deviations for the fitted parameters (mu and 
%       sigma only). Each row corresponds to a different noise level for 
%       this cue.
%   params_CI: confidence intervals for each parameter.
%   noiselevels: The noise levels for this cue.
%            
%   Structure subfields for cueCombos:
%   params: stores the fitted mu,sigma and lapserate. For cueCombos params
%       is a 3D array with the first dimension representing the different
%       noise levels involved, the second dimension the levels of conflict
%       and the third the three different parameters.
%   params_STD: The standard deviation of the fitted parameters (mu and 
%       sigma only). The 3D structure is the same as for params.
%   params_CI: confidence intervals for each parameter.
%   prediction: The predicted mu and sigmas. The 3D array has the same
%       format as params and se.
%   prediction_STD: standard deviations for the predictions.
%   prediction_CI: confidence intervals for the predictions.
%   noiselevels: the noise levels involded. This is a 2D array with
%       dimensions noiselevels, cues.
%   noiselevelIDs: the indeces for the noise levels involded. This is a 2D 
%       array with dimensions noiselevels, cues. This array is useful for
%       finding the corresponding unisensory cues more easily
%   conflicts: the partial conflicts added to the standard for each cue.
%       This is a 3D array with dimensions noiselevels, conflicts, cues.
%   cueID: the cues-numbers included in this cue combination. This is of
%       particular use when more than 2 cues are included in the design.        
%   weights: the empirical weights for each cue in the combined conditions
%   weights_STD: standard deviations of the empirical weights
%   weights_CI: Confidence intervals of the empirical weights
%   predicted_weights: the predicted weights for each cue in the combined
%       conditions
%   predicted_weights_STD: standard deviations of the predicted weights
%   predicted_weights_CI: Confidence intervals of the predicted weights
%   sigma: the empirical sigma for the combined conditions across conflicts
%   sigma_STD: standard deviations of the empirical sigma across conflicts
%   sigma_CI: Confidence intervals of the empirical sigma across conflicts
%   predicted_sigma: the predicted sigma for the combined conditions
%   predicted_sigma_STD: standard deviations of the predicted sigma
%   predicted_sigma_CI: Confidence intervals of the predicted sigma
%
%
%   A second optional output is flagLargeJND. This output specifies whether
%   this dataset contains very large JNDs that could not be estimated
%   reliably. If the JND > half the measurement range for any of the fits
%   in this dataset the flag is set to 1. Otherwise it will be zero.
%
%   Note:
%   To measure optimal integration there are two things of importance:
%   1: The PSE (Point of Subjective Equality) for the combined percept
%   should be a weighted average of the single cues alone. The weights in
%   this case should correspond to the relative reliability.
%   2: The JND (Just Noticable Difference) should be better than either cue
%   alone according to
%   sigma_c^2 = sigma_1^2 * sigma_2^2 / ( sigma_1^2 + sigma_2^2 )
%   or since sigma relates to the JND by a constant also
%   JNDc^2 = JND1^2 * JND2^2 / ( JND1^2 + JND2^2 ).
%
%
%   by Loes van Dam, 2014
%   Loes.van_Dam@uni-bielefeld.de
%
% -------------------------------------------------------------------------

if nargin < 2 || isempty(nrStim),
    nrStim = 2;
end
if nrStim ~= 1 && nrStim ~= 2,
    error('nrStim does not seem to have an appriate value. Should be either 1 or 2.');
end

if nargin < 3 || isempty(fitlapses),
    fitlapses = 0;
end

if nargin < 4 || isempty(comb_stand),
    comb_stand = 0;
end

if nargin < 5 || isempty(plotopts),
    plotopts = [0,0];
end
if length(plotopts)<2,
    plotopts = [plotopts,0];
end

if nargin < 6 || isempty(verbose),
    verbose = 0;
end

sizeA = size(A);
% if verbose,
%     fprintf(...
%         ['Single participant analysis\n',...
%         'Will now try to infer number of cues from the size of A\n'...
%         '...\n']);
% end
nrcues = (sizeA(2)-4)/3;
if nrcues == round(nrcues),
    % succes
%     if verbose,
%         fprintf('Seems like nrcues = %d\n',nrcues);
%         fprintf('Will continue with this value\n');
%     end
else
    % fail
    if verbose,
        fprintf(...
            ['The matrix A does not have the appropriate format.\n',...
             'Did you include all the necessary variables?\n'...
            'Type "help MLE_SingleSubject" for more information\n']);
    end
    error([ 'Error: Single participant analysis. ',...
            'The matrix A does not seem to have the appropriate format.\n']);

end


% -------------------------------------------------------------------------
%   Extract nr of Standards, Conflicts, Noise levels
% -------------------------------------------------------------------------

cueID = 4+(1:nrcues);
conflictID = max(cueID)+(1:nrcues);
noiseID = max(conflictID)+(1:nrcues);

standard_conditions = unique(A(:,4));       % how many standard conditions were there

nrStandards = length(standard_conditions');

assumption_test = 1;

if nrStandards > 1 && comb_stand == 0,
    % -------------------------------------------------------------------------
    %   Do MLE analysis for each standard
    % -------------------------------------------------------------------------
    
    for base = 1:nrStandards,
        B = A(A(:,4) == standard_conditions(base),:); % get the trials corresponding to the base slant

        base_results{base} = MLE_SingleSet(B, nrStim, fitlapses, [0,0], 0); % force no plotting and no verbose for this first check
    end

    % -------------------------------------------------------------------------
    %   Compare the results for the different standards and test for
    %   differences. To reduce the number of tests we do this only for the
    %   single cue conditions.
    % -------------------------------------------------------------------------

    SCnoiselevels = 0;
    for cue = 1:nrcues,
        SCnoiselevels = SCnoiselevels + length(base_results{1}.singleCues{cue}.noiselevels);
    end
    
    alpha = 0.05;
    nrtests = sum(1:(nrStandards-1)) * 2 * SCnoiselevels; % for Bonferroni correction
    alpha = alpha/nrtests; % Bonferroni corrected alpha
    
    for cue = 1:nrcues,
        for j = 1:nrStandards,
            for k = (j+1):nrStandards,
                paramdiff = base_results{j}.singleCues{cue}.params(:,1:2) -...
                            base_results{k}.singleCues{cue}.params(:,1:2);
                paramSe   = sqrt(base_results{j}.singleCues{cue}.params_STD.^2 +...
                            base_results{k}.singleCues{cue}.params_STD.^2);
                        
                % calculate upper and lower confidence intervals
                upperCI   = norminv(1-alpha/2)*paramSe + paramdiff;
                lowerCI   = norminv(0+alpha/2)*paramSe + paramdiff;
                
                testvals = ~( ( (upperCI > 0) & (lowerCI < 0) ) | isnan(paramSe) ); % if 0 lies within the confidence interval it is okay otherwise our assumption is rejected.
                                
                % check if the null hypothesis of same values is rejected
                if sum(testvals(:)) > 0,
                    assumption_test = 0; % if so then our assumption is invalid
                    
                    fprintf([
                        'WARNING: the assumption of same parameters for different standards ',...
                        'was violated on an individual participant basis for cue = %d.\n'], cue);
                end                
            end
        end
    end
    
end

% -------------------------------------------------------------------------
%   Perform the analysis for the standards combined (i.e. pooled)
% -------------------------------------------------------------------------
if assumption_test == 1,
    [results,flagLargeJND] = MLE_SingleSet(A, nrStim, fitlapses, plotopts, verbose);
    results.assumption_standardssame = assumption_test;
else
    error(['Results for the unicue standards violates same parameters assumption. ',...
        'Change the analysis to compute the results for each standard separately.']);
end