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

% -------------------------------------------------------------------------
%   [results] = MLE_SingleSet(A, nrStim, fitlapses, ploton, verbose)
%
%   Performs the MLE analysis on the data for a single set of results (e.g.
%   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% 'comparison chosen' 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% 'comparison chosen' 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.
%
%   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(plotopts),
    plotopts = [0,0];
end
if length(plotopts)<2,
    plotopts = [plotopts,0];
end

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

sizeA = size(A);
if verbose,
    fprintf(...
        ['Single dataset 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\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_Analysis" for more information\n']);
    end
    error([ 'Error: Single participant analysis. ',...
            'The matrix A does not seem to have the appropriate format.\n']);

end

measRange = diff([min(A(:,1)),max(A(:,1))])/1.0;
flagLargeJND = 0;

% -------------------------------------------------------------------------
%   specify nr of participants in this dataset
% -------------------------------------------------------------------------

results.nrSub = [1,1];

% -------------------------------------------------------------------------
%   Data in single trials or [successes, #trials] (for plotting purposes)
% -------------------------------------------------------------------------

singletrials = (sum(A(:,3)) == length(A(:,3)));

% -------------------------------------------------------------------------
%   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
noiseLevels = unique(A(:,noiseID),'rows');

nrStandards = length(standard_conditions');
nrNoises = length(noiseLevels(:,1));

% -------------------------------------------------------------------------
%   If requested: Get overall lapserate across all data
% -------------------------------------------------------------------------

if fitlapses == 1,
    [mu,sigma,lapserate] = FitCumulativeGauss(A(:,1:3));  % fit the psychometric function using GLM regression
else
    lapserate = 0;
end


% -------------------------------------------------------------------------
%   Set how to compute confidence intervals:
%   using norminv (assuming normal distribution): CI_NormT = 0
%   using tinv (assuming student t distribution): CI_NormT = 1
% -------------------------------------------------------------------------

CI_NormT = 0;

% -------------------------------------------------------------------------
%   Do MLE analysis
%
%   Note that here we do not separate the data for different standards.
%   This is because depending on the design the standards will either be
%   pooled together to form one big dataset (if assumed there is no
%   difference in the standard) or performed separately for each standard
%   (if assumed there is an effect of standard). Therefore, here the data
%   is only splitted for noise and conflict levels. If a split for
%   different standards is desired, this should be done in the function
%   calling this function. Alternatively, they can also be added in the
%   design as addition "noise levels".
%
% -------------------------------------------------------------------------

B = A;

% -------------------------------------------------------------------------
%   Sort out the trials for single cue conditions
% -------------------------------------------------------------------------

SC = B(sum(B(:,cueID)')' == 1,:);                               % get all the single cue conditions

for cue = 1:nrcues,
    C = SC(SC(:,cueID(cue))==1,:);                              % get all data for that particular cue

    CNoise = unique(C(:,noiseID(cue)));                         % get the noise levels for that cue

    for NL = 1:length(CNoise),                                  % loop through the different noise levels
        CN = C(C(:,noiseID(cue))==CNoise(NL),:);                % get all data for the noise level within the single cue

        [mu,sigma,delme,seMuSig] = FitCumulativeGauss(CN(:,1:3),lapserate);     % fit the psychometric function using GLM regression using the lapse rate 'lapserate'
        if abs(sigma) > measRange,
            if verbose,
                fprintf(['WARNING: the JND is larger than the measurement range used.\n'...
                    'Condition: Single cue = %d Noise = %f\n'],cue,CNoise(NL));
            end
            flagLargeJND = 1;
        end
        % correct the sigma for the number of stimuli used in each trial
        sigma = sigma/sqrt(nrStim);
        seMuSig(2) = seMuSig(2)/sqrt(nrStim);
        % put the fitted mean and sigma of the Gaussian distribution in
        % the results list: added also are the lapse rate used for the fit and the noise level
        
        results.singleCues{cue}.params(NL,:)        = [mu,sigma,lapserate];  
        results.singleCues{cue}.noiselevels(NL)     = CNoise(NL);
        results.singleCues{cue}.params_STD(NL,:)           = seMuSig';
        if CI_NormT == 0,
            results.singleCues{cue}.params_CI(NL,:,:)          = [mu,sigma]'*ones(1,2)+norminv(0.975)*seMuSig'*[-1,1];
        else
            results.singleCues{cue}.params_CI(NL,:,:)          = [mu,sigma]'*ones(1,2)+tinv(0.975,1)*seMuSig'*[-1,1];
        end
        
        % -----------------------------------------------------------------
        %   Sort the data for plotting purposes
        % -----------------------------------------------------------------
        if singletrials,
            IVs = unique(CN(:,1));
            for k = 1:length(IVs),
                IVtrials = CN(CN(:,1) == IVs(k),2);
                data(k,:) = [IVs(k),sum(IVtrials == 1),length(IVtrials)];
            end
        else
            data = CN(:,1:3);
        end
        results.singleCues{cue}.rawdata{NL} = data;

    end

end

% -------------------------------------------------------------------------
%   Sort out the trials for all possible cue combinations
% -------------------------------------------------------------------------


CC = B(sum(B(:,cueID)')'>1,:);          % get multimodal trials to find conflict conditions

cueCombos = unique(CC(:,cueID),'rows'); % make a list of the different cue combinations
nrCombos = length(cueCombos(:,1));      % the number of cue combinations


maxNrNoises = 0;

for cueC = 1:nrCombos,

    cueC_ID = cueID(cueCombos(cueC,:) == 1);
    cueN_ID = cueID(cueCombos(cueC,:) == 0);

    C = CC(sum(CC(:,cueC_ID)')' == sum(cueCombos(cueC,:)') & sum(CC(:,cueN_ID)')' == 0,:); % get trials for that particular cue combination

    % set elements to look at further particulars for the cue
    % combination
    noiseC_ID = noiseID(cueCombos(cueC,:) == 1);
    conflictC_ID = conflictID(cueCombos(cueC,:) == 1);


    % get levels of noise
    CNoise = unique(C(:,noiseC_ID),'rows');    % get the noise levels for that cue combination

    if length(CNoise(:,1)) > maxNrNoises,
        maxNrNoises = length(CNoise(:,1));     % for plotting purposes sort out the max number of noise combos
    end
    
    results.cueCombos{cueC}.noiselevels = CNoise;

    for NL = 1:length(CNoise(:,1)),          % loop through the different noise levels

        evalstring = [];
        for ce = 1:length(cueC_ID),
            evalstring = [evalstring, ['C(:,',num2str(noiseC_ID(ce)),') == CNoise(NL,',num2str(cueC_ID(ce)-4),')']];
            if ce < length(cueC_ID)
                evalstring = [evalstring, '&'];
            end
        end

        eval(['CN = C(',evalstring,',:);']); % get data for that combination of cues and noise levels

        % get levels of conflict
        Cconf = unique(CN(:,conflictC_ID),'rows');


        % -------------------------------------------------------------------------
        %   find the corresponding unisensory condition
        % -------------------------------------------------------------------------

        for cue = 1:nrcues,
            onecuestuff(cue,:)    = results.singleCues{cue}.params( results.singleCues{cue}.noiselevels==CNoise(NL,cue),:);
            onecuestuffStd(cue,:) = results.singleCues{cue}.params_STD(    results.singleCues{cue}.noiselevels==CNoise(NL,cue),:);
            results.cueCombos{cueC}.noiselevelIDs(NL,cue) = find(results.singleCues{cue}.noiselevels==CNoise(NL,cue)); 
        end

        forRegression = [];
        for CF = 1:length(Cconf(:,1)),

            % -------------------------------------------------------------------------
            %   make the predictions
            % -------------------------------------------------------------------------

            predicted_JND = sqrt(1/sum(1./onecuestuff(:,2).^2));
            predicted_PSE = sum(Cconf(CF,:)'./(onecuestuff(:,2).^2))/sum(1./onecuestuff(:,2).^2);
            predicted_weights = 1./(onecuestuff(:,2).^2)/sum(1./onecuestuff(:,2).^2);

            % for each of these compute the standard devation using the
            % delta method.            
            predicted_JND_std = sqrt( (1/sum(1./onecuestuff(:,2).^2)^3) * sum( (1./onecuestuff(:,2).^6).*(onecuestuffStd(:,2).^2)  ) );
            predicted_weights_std = sqrt( (predicted_weights.^2)*...
                        sum( (2*predicted_weights.*onecuestuffStd(:,2)./onecuestuff(:,2) ).^2) +...
                        (1-2*predicted_weights).*(2*predicted_weights./onecuestuff(:,2) ).^2 .*onecuestuffStd(:,2).^2  );
            scaledE = Cconf(CF,:)'.*predicted_weights;
            predicted_PSE_std = sqrt( ...
                        sum(sum(scaledE*scaledE'))*sum( (2*predicted_weights.*onecuestuffStd(:,2)./onecuestuff(:,2) ).^2) +...
                     -2*sum(sum(scaledE*(scaledE.*predicted_weights.*(2.*onecuestuffStd(:,2)./onecuestuff(:,2) ).^2)'))+...                    
                        sum(scaledE.^2.*(2./onecuestuff(:,2) ).^2 .*onecuestuffStd(:,2).^2 ) );
                            
            evalstring = [];
            for ce = 1:length(cueC_ID),
                evalstring = [evalstring, ['CN(:,',num2str(conflictC_ID(ce)),') == Cconf(CF,',num2str(cueC_ID(ce)-4),')']];
                if ce < length(cueC_ID)
                    evalstring = [evalstring, '&'];
                end
            end

            eval(['CNCF = CN(',evalstring,',:);']); % get data for that combination of cues, noise levels and conflict

            % do the fitting
            [mu,sigma,delme,seMuSig] = FitCumulativeGauss(CNCF(:,1:3),lapserate);     % fit the psychometric function using GLM regression using the lapse rate 'lapserate'
            if abs(sigma) > measRange,
                if verbose,
                    strd = [];
                    strf = [];
                    for ce = 1:length(cueC_ID),
                        strd = [strd,' %d'];
                        strf = [strf,' %f'];
                    end                
                    fprintf(['WARNING: the JND is larger than the measurement range used.\n'...
                        'Condition: cue combo = (',strd,' ) Noise = (',strf,' ) Conflict: (',strf,' )\n'],cueC_ID,CNoise(NL,:),Cconf(CF,:));
                end
                flagLargeJND = 1;
            end
            % correct sigma with respect to the number of stimuli in each
            % trial
            sigma = sigma/sqrt(nrStim);
            seMuSig(2) = seMuSig(2)/sqrt(nrStim);
            % put the fitted mean and sigma of the Gaussian distribution in
            % the results list: added also are the lapse rate used for the fit and the noise level
            results.cueCombos{cueC}.params(NL,CF,:) = [mu,sigma,lapserate];
            results.cueCombos{cueC}.conflicts(NL,CF,:) = Cconf(CF,:);
            results.cueCombos{cueC}.params_STD(NL,CF,:) = seMuSig';
            if CI_NormT == 0,
                results.cueCombos{cueC}.params_CI(NL,CF,:,:) = [mu,sigma]'*ones(1,2)+norminv(0.975)*seMuSig'*[-1,1];
            else
                results.cueCombos{cueC}.params_CI(NL,CF,:,:) = [mu,sigma]'*ones(1,2)+tinv(0.975,1)*seMuSig'*[-1,1];
            end
            results.cueCombos{cueC}.cueID = cueC_ID-4;
            

            results.cueCombos{cueC}.prediction(NL,CF,:) = [predicted_PSE,predicted_JND];
            results.cueCombos{cueC}.prediction_STD(NL,CF,:) = [predicted_PSE_std,predicted_JND_std];
            if CI_NormT == 0,
                results.cueCombos{cueC}.prediction_CI(NL,CF,:,:) = [predicted_PSE;    predicted_JND     ]*ones(1,2)+...
                                                     norminv(0.975)*[predicted_PSE_std;predicted_JND_std ]*[-1,1];
                results.cueCombos{cueC}.predicted_weights_CI(NL,:,:) = predicted_weights*ones(1,2)+norminv(0.975)*predicted_weights_std*[-1,1];
            else
                results.cueCombos{cueC}.prediction_CI(NL,CF,:,:) = [predicted_PSE;    predicted_JND     ]*ones(1,2)+...
                                                     tinv(0.975,1)*[predicted_PSE_std;predicted_JND_std ]*[-1,1];
                results.cueCombos{cueC}.predicted_weights_CI(NL,:,:) = predicted_weights*ones(1,2)+tinv(0.975,1)*predicted_weights_std*[-1,1];
            end
            results.cueCombos{cueC}.predicted_weights(NL,:) = predicted_weights;
            results.cueCombos{cueC}.predicted_weights_STD(NL,:) = predicted_weights_std;


            % -------------------------------------------------------------------------
            %   Sort the data for plotting and plot results
            % -------------------------------------------------------------------------

            if singletrials,
                IVs = unique(CNCF(:,1));
                for k = 1:length(IVs),
                    IVtrials = CNCF(CNCF(:,1) == IVs(k),2);
                    data(k,:) = [IVs(k),sum(IVtrials == 1),length(IVtrials)];
                end
            else
                data = CNCF(:,1:3);
            end
            results.cueCombos{cueC}.rawdata{NL,CF} = data;

        end
        
        
        % -------------------------------------------------------------------------
        % mostly for plotting: compute average JND across conflicts but
        % within the same noise level. The JND should be the same for all
        % conflicts and this makes it easier for plotting later.
        % -------------------------------------------------------------------------
        nrconflicts = size(results.cueCombos{cueC}.params_STD(:,:,2),2);
        tempJND = squeeze(results.cueCombos{cueC}.params(NL,:,2));
        results.cueCombos{cueC}.sigma(NL,1) = mean(tempJND);
        if nrconflicts > 1,
            % compute emperical standard deviation
            results.cueCombos{cueC}.sigma_STD(NL,1) = std(tempJND);
        else
            % compute standard deviation through error propagation.
            results.cueCombos{cueC}.sigma_STD(NL,1) = sqrt(sum(squeeze(results.cueCombos{cueC}.std(NL,:,2)).^2))/nrconflicts;
        end
        if CI_NormT == 0,
            results.cueCombos{cueC}.sigma_CI(NL,:) = mean(tempJND)'*ones(1,2)+...
                norminv(0.975)*std(tempJND)*[-1,1]/sqrt(nrconflicts);
        else
            if nrconflicts == 1,
                results.cueCombos{cueC}.sigma_CI(NL,:) = mean(tempJND)'*ones(1,2)+...
                    tinv(0.975,nrconflicts)*std(tempJND)*[-1,1];
            else
                results.cueCombos{cueC}.sigma_CI(NL,:) = mean(tempJND)'*ones(1,2)+...
                    tinv(0.975,nrconflicts-1)*std(tempJND)*[-1,1]/sqrt(nrconflicts);
            end
        end

        % same for predictions although here they are already the same so
        % just take first one ;-)
        results.cueCombos{cueC}.predicted_sigma(NL,1) = squeeze(results.cueCombos{cueC}.prediction(NL,1,2));
        results.cueCombos{cueC}.predicted_sigma_STD(NL,1) = squeeze(results.cueCombos{cueC}.prediction_STD(NL,1,2));
        results.cueCombos{cueC}.predicted_sigma_CI(NL,:)  = squeeze(results.cueCombos{cueC}.prediction_CI(NL,1,2,:));

        
        % -------------------------------------------------------------------------
        % determine emperical cue weights from the PSEs
        % -------------------------------------------------------------------------
        absconf = abs(diff(Cconf'));
        if length(cueC_ID)>2,
            absconf = sum(absconf);
        end
        nrconf = sum(absconf ~= 0);         % zero conflict conditions can not provide info about the weights from the pse
        % -------------------------------------------------------------------------
        % first check if determining weights is at all possible
        %
        % To solve a set of equations (which is what a regression does in a
        % way), you need at least as many equations as there variables to
        % estimate. Each conflict can in principle represent such an
        % equation if the conflicts are chosen wisely. Added to that is the
        % equation sum(weights) = 1.
        % -------------------------------------------------------------------------
        if nrconf < length(cueC_ID)-1,      % there is not enough data to do the regression
            warning('There are not enough conflicts to unambiguously determine the cue weights from the PSEs. Skipping this part.');
            results.cueCombos{cueC}.weights(NL,:) = NaN*ones(1,length(cueC_ID)); 
            results.cueCombos{cueC}.weights_STD(NL,:) = NaN*ones(1,length(cueC_ID));
            results.cueCombos{cueC}.weights_CI(NL,:,:) = NaN*ones(2,length(cueC_ID));
        else
            % get the variables for the regression
            forRegression = [];
            for CF = 1:length(Cconf(:,1)),
                % check if nonzero and if so add to the list
                pse = results.cueCombos{cueC}.params(NL,CF,1);
                cuevals = Cconf(CF,:);
                if absconf(CF) ~= 0,
                    forRegression = [forRegression;[pse-cuevals(end),cuevals(1:(end-1))-cuevals(end)]]; %add the regression variables
                end
            end
            try
                statsweight = regstats(forRegression(:,1),forRegression(:,2:end),eye(length(forRegression(1,:))-1));
                weights = statsweight.beta;
                weights(end+1) = 1 - sum(weights);                          % add the weight for the last cue
                weightsstd = diag(sqrt(statsweight.covb));                  % get standard deviation
                weightsstd(end+1) = sqrt(sum(statsweight.covb(:)));         % add last element of standard deviation

                results.cueCombos{cueC}.weights(NL,:) = weights;
                results.cueCombos{cueC}.weights_STD(NL,:) = weightsstd;
                
                
                sw = size(weights);
                if sw(2)> sw(1),
                    weights = weights';
                    weightsstd = weightsstd';
                end
                
                if CI_NormT == 0,
                    results.cueCombos{cueC}.weights_CI(NL,:,:) = weights*ones(1,2)+norminv(0.975)*weightsstd*[-1,1];
                else
                    results.cueCombos{cueC}.weights_CI(NL,:,:) = weights*ones(1,2)+tinv(0.975,1)*weightsstd*[-1,1];
                end
            catch
                results.cueCombos{cueC}.weights(NL,:) = NaN*ones(1,length(cueC_ID));
                results.cueCombos{cueC}.weights_STD(NL,:) = NaN*ones(1,length(cueC_ID)); 
                results.cueCombos{cueC}.weights_CI(NL,:,:) = NaN*ones(2,length(cueC_ID));
            end
        end
        
    end
end

% -------------------------------------------------------------------------
%   Here we will plot the results
%
%   There are two plotoptions in plotopts specified below
%   1st: if 1 the individual curbes are plotted for each fit in several
%   subplots.
%   The organisation of the subplots is such that each pair of two rows
%   represents a particular cue-combination in the dataset. 
%   The top of these two represents the fits for single cues, the second
%   row for cue combination conditions.
%   The colums of the figure represent te different noise levels for that
%   combination. Thus if 1 noise level in one cue is combined with several
%   levels for the other cue, it will be the case that certain fits are
%   repeatedly represented (the one for the cue with only one noise level
%
%   2nd plotoption: if 1 a summary plot will be made (see below for more
%   details.
% -------------------------------------------------------------------------


colstring = 'brmgcyk';

if plotopts(1),
    nrplotcol = maxNrNoises;
    nrplotrow = nrCombos*2;

    plotrange = [min(B(:,1)),max(B(:,1)),0,1];
    IVs = unique(B(:,1));
    IVscale = IVs(1)+(0:100)*(IVs(end)-IVs(1))/100;

    figure();
    for cueC = 1:nrCombos,
        noiselevels = results.cueCombos{cueC}.noiselevels;
        nrNoises = length(noiselevels(:,1));
        nrCueInCombo = length(results.cueCombos{cueC}.cueID);
        for NL = 1:nrNoises,
            
            % make the subplot for single cue conditions
            subplot(nrplotrow,nrplotcol,(cueC-1)*2*nrplotcol+NL);
            curve_count = 0;
            plotstr = [];
            for cue = 1:nrCueInCombo,
                cueID = results.cueCombos{cueC}.cueID(cue);
                cueNoiseID = (results.singleCues{cue}.noiselevels==noiselevels(NL,cue));
                cp = results.singleCues{cue}.params(cueNoiseID,:);
                
                curve_count = curve_count + 1;
                if curve_count > length(colstring),
                    curve_count = 1;
                end

                data = results.singleCues{cue}.rawdata{cueNoiseID};
                plot(data(:,1),data(:,2)./data(:,3),['.',colstring(curve_count)]); hold on;
                
                h(cue) = plot(IVscale,cp(3)+(1-2*cp(3))*normcdf(IVscale,cp(1),sqrt(nrStim)*cp(2)),['-',colstring(curve_count)]);
                axis(plotrange);                
                
                plotstr = [plotstr,'(C: ', num2str(cueID), ' N: ', num2str(noiselevels(NL,cue)),')'];
                legstr{cue} = ['C: ', num2str(cueID)];
            end
            legend(h,legstr,'location','northwest');
            title( plotstr );
            xlabel('dependent variable');
            ylabel('proportion comparison chosen');
            
            % make the subplot for cue combination conditions
            subplot(nrplotrow,nrplotcol,((cueC-1)*2+1)*nrplotcol+NL); hold on;
            nrConf = length(results.cueCombos{cueC}.conflicts(NL,:,1));
            curve_count = 0;
            
            for CF = 1:nrConf,
                curve_count = curve_count + 1;
                if curve_count > length(colstring),
                    curve_count = 1;
                end
                
                data = results.cueCombos{cueC}.rawdata{NL,CF};
                cp = results.cueCombos{cueC}.params(NL,CF,:);
                cpp = results.cueCombos{cueC}.prediction(NL,CF,:);
                
                for cue = 1:nrCueInCombo,
                    plot(results.cueCombos{cueC}.conflicts(NL,CF,cue)*[1,1],[0,1],'--k');
                end
                plot(plotrange(1:2),[0.5,0.5],'--k');
                plot(data(:,1),data(:,2)./data(:,3),['.',colstring(curve_count)]); hold on;
                plot(IVscale,cp(3)+(1-2*cp(3))*normcdf(IVscale,cp(1),sqrt(nrStim)*cp(2)),['-',colstring(curve_count)]);
                plot(IVscale,normcdf(IVscale,cpp(1),sqrt(nrStim)*cpp(2)),['--',colstring(curve_count)]);
                plot(cpp(1),0.5,['d',colstring(curve_count)]);
                axis(plotrange);

            end
            title('combined');
            xlabel('dependent variable');
            ylabel('proportion comparison chosen');
            
        end
    end
end





% -------------------------------------------------------------------------
%   if plotopts(2) == 1, make summary plot
% -------------------------------------------------------------------------
if plotopts(2) == 1,
    MLE_SummaryPlot(results);
end