function data = RunExperiment(inputparams,totalsubjects,whichtasks,plot_stuff,trialcount,experiment_count)

%    This file is part of RunExperiment.
%
%    RunExperiment is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    RunExperiment is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with RunExperiment.  If not, see <http://www.gnu.org/licenses/>.


% RunExperiment is used to collapse the data from multiple pseudo-subjects,
% as in an experiment.

addpath Functions

if nargin == 0
    
    % If you don't want to use one of the pre-designed experiments in, say,
    % RunAPP_data, you are welcome to input different parameters here. For
    % a list of parameters, check out SetParameters.
    inputparams{1,1} = 'parameter.trialcount';
    inputparams{1,2} = 500;
    
    inputparams{2,1} = 'parameter.bp_size';
    inputparams{2,2} = 800;
    
    inputparams{3,1} = 'parameter.token_conn';
    inputparams{3,2} = .45;
    
    inputparams{4,1} = 'parameter.type_conn';
    inputparams{4,2} = .45;
    
    inputparams{5,1} = 'parameter.type_2_conn';
    inputparams{5,2} = .45;
    
    inputparams{6,1} = 'parameter.type_1_overlap';
    inputparams{6,2} = .3;

    inputparams{7,1} = 'parameter.token_individuation_threshold';
    inputparams{7,2} = .016;
    
    cr_setsize = [1 2 4 6]; % cr = continuous report
    cd_setsize = [2 4 6 8]; % cd = change detection
    
    totalsubjects = 5;
    whichtasks = 1; % Task 1 = Continuous report, Tasl 2 = Change detection
    
    plot_stuff = 1;
    
    experiment_count = 1;
    
elseif nargin==1 
    
    cr_setsize = [1 2 4 6];
    cd_setsize = [2 4 6 8];
    
    totalsubjects = 5;
    whichtasks = 1;
    
    plot_stuff = 1;
    
    experiment_count = 1;
    
    
else
    
    if isstr(inputparams{1,2})
        cr_setsize = str2num(inputparams{1,2});
        cd_setsize = str2num(inputparams{1,2});
    else
        cr_setsize = inputparams{1,2};
        cd_setsize = inputparams{1,2};
    end
    
    inputparams{40,1} = 'parameter.trialcount';
    if isstr(trialcount)
        inputparams{40,2} = str2num(trialcount);
    else
        inputparams{40,2} = trialcount;
    end
    
end

standard_error = 0;

% RUN MODEL

for subnum = 1:totalsubjects
    
    for task = whichtasks
        
        if task == 1 % Always continuous report
            
            [cr(subnum),parameter] = RunModel(inputparams,1,cr_setsize);
            
        elseif task ==2 % Always change detection
            [cd(subnum),parameter]  = RunModel(inputparams,2,cd_setsize);
            
        elseif task == 0 % neither, but useful for other things
            
            [cr(subnum),parameter] = RunModel(inputparams,2,inputparams{1,2});
            
        end
    end
    
end

% PUTTING-THEM-ALTOGETHER

for task = whichtasks
    
    if task == 1
        setsize = cr_setsize;
    elseif task == 2
        setsize = cd_setsize;
    elseif task == 0
        setsize = inputparams{1,2};
    end
    
    setsize_count = 0;
    
    for load = setsize
        
        setsize_count = setsize_count + 1;
        
        cr_count = 0;
        cd_count = 0;
        all_con_count = 0;
        
        guess_count = 0;
        correct_count = 0;
        incorrect_count = 0;
        
        x = [];
        t = [];
        nt = [];
        
        
        for subnum = 1:totalsubjects
            for trial = 1:parameter.trialcount 
                
                if task == 1
                    cr_count = cr_count + 1;
                    
                    if parameter.fixed_target_item == 0
                        target(trial) = randi(load);
                    else
                        target(trial) = parameter.fixed_target_item;
                    end
 
%                     target_token(cr_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).retrieved_token_representation(target(trial));
%                     other_target_token(cr_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).retrieved_token_representation(mod(target(trial),2)+1);
%        
%                     target_token_nt(cr_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,mod(target(trial),2)+1).retrieved_token_representation(target(trial));
%                     other_target_token_nt(cr_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,mod(target(trial),2)+1).retrieved_token_representation(mod(target(trial),2)+1);
%                     
                    x(cr_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).retrieved_location;
                    t(cr_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).target_radians;
                     
%                     error(cr_count) = x(cr_count)-t(cr_count);
%                     
                    if load ~= 1 
                        nt(:,cr_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).nontarget_radians;
                    end
%                     
%                     if cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).token == cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).retrieved_token
%                         
%                         correct_count = correct_count + 1;
%                          
%                         x_c(correct_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).retrieved_location;
%                         t_c(correct_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).target_radians;
%                          
%                         correct_error(correct_count) = x(correct_count)-t(correct_count);
%                         
%                         if load ~= 1 
%                             nt_c(:,correct_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).nontarget_radians;
%                         end
%                         
%                     elseif cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).retrieved_token == 0
%                         
%                         guess_count = guess_count + 1;
%                         
%                         x_g(guess_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).retrieved_location;
%                         t_g(guess_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).target_radians;
%                         
%                         guess_error(guess_count) = x_g(guess_count)-t_g(guess_count);
%                          
%                         if load ~= 1
%                             nt_g(:,guess_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).target_radians;
%                         end
%                         
%                     else 
%                         
%                         incorrect_count = incorrect_count + 1;
%                         
% %                         x_b(incorrect_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).retrieved_location;
% %                         t_b(incorrect_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).nontarget_radians;
% %                         
% %                         incorrect_error(incorrect_count) = x_b(incorrect_count)-t_b(incorrect_count);
% %                         
% %                         if load ~= 1
% %                             nt_b(:,incorrect_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).target_radians;
% %                         end
% %                         
%                          
%                     end
%                     
%                                         
                    
                   
                    
                    
                    if load == 1
                        for con_load = setsize
                            all_con_count = all_con_count + 1;
                            retrieved_length(all_con_count) = cr(subnum).Model_sim(con_load).trial_datastore.data(trial,target(trial)).retrieved_length;
                        end
                    end
                    
                elseif task == 2
                    
                    cd_count = cd_count + 1;
                    
                    for load_again = 1:load
                        
                        cd_stuff(load).deviation(cd_count,load_again) = cd(subnum).Model_sim(load).deviation(trial,load_again);
                        cd_stuff(load).retrieved_length(cd_count,load_again) = cd(subnum).Model_sim(load).retrieved_length(trial,load_again);
                        cd_stuff(load).changed_items(cd_count) = cd(subnum).Model_sim(load).changed_items(trial);
                        cd_stuff(load).old_to_new(cd_count,load_again) = cd(subnum).Model_sim(load).old_to_new(trial,load_again);
                        
                    end
                    
                end
            end
        end
%         
%          
%          figure
%          hist(x,36);
%          
%          keyboard
%          
%         
        
        
        if task == 1 % Continous report
            
            if parameter.confidence == 1 % confidence
                
                bins = prctile(retrieved_length,[33 66]);
                consplit_tertile1 = bins(1);
                consplit_tertile2 = bins(2);
                
                for confidence_category = 1:3
                    
                    con_count = 0;
                    
                    if confidence_category == 1
                        consplit_min = 0;
                        consplit_max = consplit_tertile1;
                    elseif confidence_category == 2
                        consplit_min = consplit_tertile1;
                        consplit_max = consplit_tertile2;
                    else
                        consplit_min = consplit_tertile2;
                        consplit_max = 10000;
                    end
                    
                    x_c = [];
                    t_c = [];
                    nt_c = [];
                    
                    for subnum = 1:totalsubjects
                        for trial = 1:parameter.trialcount
                            
                            trial_length = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).retrieved_length;
                            
                            if trial_length < consplit_max && trial_length > consplit_min
                                con_count = con_count + 1;
                                
                                x_c(con_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).retrieved_location;
                                t_c(con_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).target_radians;
                                
                                
                                if load > 1
                                    nt_c(:,con_count) = cr(subnum).Model_sim(load).trial_datastore.data(trial,target(trial)).nontarget_radians;
                                end
                                
                            end
                        end
                    end
                    
                    if isempty(x_c) == 0
                        if load > 1
                            output = JV10_fit(x_c',t_c',nt_c');
                        elseif load == 1
                            output = JV10_fit(x_c',t_c');
                        end
                    else
                        output = [0 0 0 0];
                    end
                    
                    c(confidence_category).error(setsize_count) = rad2deg(sqrt(1/output(1)));
                    c(confidence_category).swap(setsize_count) = output(3);
                    c(confidence_category).guess(setsize_count) = output(4);
                    
                end
                
            else % Not confidence Continuous report
                
                
                if load > 1
                    output = JV10_fit(x',t',nt');
                else
                    output = JV10_fit(x',t');
                end
                
                
                
                
                error(setsize_count) = rad2deg(sqrt(1/output(1)));
                swap(setsize_count) = output(3);
                guess(setsize_count) = output(4);
                
                x = [];
                t = [];
                nt = [];
                %    target = [];
                
                
                
                
                
                
                
            end
        end
    end
    
    if task == 2 % Change detection
        
        parameter.trialcount = parameter.trialcount * totalsubjects;
        parameter.setsize = cd_setsize;
        
        [Model_sim cd_sim parameter] = Threshold_analysis(cd_stuff,parameter);
        
        setsize_count = 0;
        
        for load = cd_setsize
            
            setsize_count = setsize_count + 1;
            
            hit(setsize_count) = (1-Model_sim(load).miss);
            false_alarm(setsize_count) = Model_sim(load).false_alarm;
            
        end
        
        prop_change = cd_sim.proportions_of_change;
        
    end
end
% PLOTTING

for task = whichtasks
    
    
    if task == 1 % Continuous Report
        if parameter.confidence == 1
            
            for category = 1:3
                data.continuous.error(category,:) = c(category).error;
                data.continuous.swap(category,:) = c(category).swap;
                data.continuous.guess(category,:) = c(category).guess;
            end
            
            if plot_stuff == 1
                if ispc || ismac
                    
                    figure
                    
                    subplot(1,3,1)
                    
                    plot(cr_setsize,data.continuous.error(1,:),'LineWidth',2,'LineStyle','--','Color','r');
                    ylim([0 30])
                    hold on
                    plot(cr_setsize,data.continuous.error(2,:),'LineWidth',2,'LineStyle','--','Color','b');
                    hold on
                    plot(cr_setsize,data.continuous.error(3,:),'LineWidth',2,'LineStyle','--','Color','g');
                    hold on
                    
                    subplot(1,3,2)
                    
                    plot(cr_setsize,data.continuous.swap(1,:),'LineWidth',2,'LineStyle','--','Color','r');
                    ylim([0 .5])
                    hold on
                    plot(cr_setsize,data.continuous.swap(2,:),'LineWidth',2,'LineStyle','--','Color','b');
                    hold on
                    plot(cr_setsize,data.continuous.swap(3,:),'LineWidth',2,'LineStyle','--','Color','g');
                    hold on
                    
                    subplot(1,3,3)
                    
                    plot(cr_setsize,data.continuous.guess(1,:),'LineWidth',2,'LineStyle','--','Color','r');
                    ylim([0 .5])
                    hold on
                    plot(cr_setsize,data.continuous.guess(2,:),'LineWidth',2,'LineStyle','--','Color','b');
                    hold on
                    plot(cr_setsize,data.continuous.guess(3,:),'LineWidth',2,'LineStyle','--','Color','g');
                    hold on
                    
                    
                end
            end
        else
            count = 0;
            
            data.continuous.error = error;
            data.continuous.guess = guess;
            data.continuous.swap = swap;
            
            
            if plot_stuff == 1
                if ispc||ismac
                    if parameter.confidence == 0
                        
                        figure
                        
                        bays_sd = [14.5 19.2 23.7 24.7];
                        bays_g = [.02 .065 .22 .215];
                        bays_b = [0  .04  .11   .31];
                        
                        subplot(1,3,1)
                        plot(cr_setsize,data.continuous.error,'LineWidth',2,'LineStyle','--');
                        hold on
                        
                        if standard_error
                            errorbar(cr_setsize,data.continuous.error,data.continuous.se_error,'LineWidth',2)
                        end
                        ylim([0 30])
                        
                        plot(cr_setsize,bays_sd,'LineWidth',2,'Color',[0 0 0]);
                        
                        subplot(1,3,2)
                        plot(cr_setsize,data.continuous.swap,'LineWidth',2,'LineStyle','--');
                        hold on
                        if standard_error
                            errorbar(cr_setsize,data.continuous.swap,data.continuous.se_swap,'LineWidth',2)
                        end
                        ylim([0 .5])
                        
                        plot(cr_setsize,bays_b,'LineWidth',2,'Color',[0 0 0]);
                        
                        subplot(1,3,3)
                        plot(cr_setsize,data.continuous.guess,'LineWidth',2,'LineStyle','--');
                        hold on
                        
                        if standard_error
                            errorbar(cr_setsize,data.continuous.guess,data.continuous.se_guess,'LineWidth',2)
                        end
                        ylim([0 .5])
                        
                        plot(cr_setsize,bays_g,'LineWidth',2,'Color',[ 0 0 0]);
                    end
                end
            end
        end
    elseif task == 2 % Change detection
        
        count = 0;
        
        data.change_detection.hit = hit;
        data.change_detection.false_alarm = false_alarm;
        data.change_detection.prop_change = prop_change;
        
        if plot_stuff == 1
            if ispc||ismac
                
                figure
                
                Keshvari_hit =                [.79   .63    .55    .5];
                Keshvari_false_alarm =        [.1    .22    .28     .34];
                
                subplot(1,2,1)
                plot(cd_setsize,data.change_detection.hit,'LineWidth',2,'LineStyle','--');
                hold on
                
                if standard_error
                    errorbar(cd_setsize,data.change_detection.hit,data.change_detection.se_hit,'LineWidth',2)
                end
                
                hold on
                
                plot(cd_setsize,data.change_detection.false_alarm,'LineWidth',2,'LineStyle','--');
                hold on
                
                plot(cd_setsize,Keshvari_hit,'LineWidth',2,'Color',[0 0 0])
                hold on
                plot(cd_setsize,Keshvari_false_alarm,'LineWidth',2,'Color',[0 0 0])
                
                legend('Model','Model','Data','Data');
                
                if standard_error
                    errorbar(cd_setsize,data.change_detection.false_alarm,data.change_detection.se_false_alarm,'LineWidth',2)
                end
                ylim([0 1])
                
                
                subplot(1,2,2)
                for i = 1:length(cd_setsize)
                    hold on
                    
                    for(bin = 2: 10)
                        if(isnan(data.change_detection.prop_change(i,bin)))  %interpolate if there is a NAN
                            
                            interp = data.change_detection.prop_change(i,bin-1)+ data.change_detection.prop_change(i,bin+1);
                            data.change_detection.prop_change(i,bin) = interp /2;
                        end
                        
                    end
                    plot(0:1:10,data.change_detection.prop_change(i,:),'LineWidth',2,'Color',[.2*i .2*i .2*i]);
                    hold on
                    
                    if standard_error
                        errorbar(0:1:10,data.change_detection.prop_change(i,:),data.change_detection.se_prop_change(i,:),'LineWidth',2,'Color',[.2*i .2*i .2*i])
                    end
                    ylim([0 1])
                end
            end
        end
    end
    
    if plot_stuff == 1
        if ispc || ismac
            if experiment_count ==1
                
                figure
                
                y_loc = .8;
                x_loc = .5;
                
                
                %uicontrol('Style','text','Units','normalized','BackgroundColor',[1 1 1],'Position',[.75 .1 .8 .8]);
                uicontrol('Style','text','String','PARAMETERS','Units','normalized','Position',[x_loc y_loc .25 .05]);
                y_loc = y_loc - .05;
                uicontrol('Style','text','String',['trial count = ' num2str(parameter.trialcount)],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                y_loc = y_loc - .05;
                uicontrol('Style','text','String',['binding pool size = ' num2str(parameter.bp_size)],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                
                if parameter.jumps == 0
                    y_loc = y_loc - .05;
                    uicontrol('Style','text','String',['capacity = Uniform ( ' num2str(parameter.capacitybottom) num2str(parameter.capacitytop) ')'],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                else
                    y_loc = y_loc - .05;
                    uicontrol('Style','text','String',['total jumps = ' num2str(parameter.jumps)],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                end
                
                y_loc = y_loc - .05;
                uicontrol('Style','text','String',['token sparsity= ' num2str(parameter.token_conn)],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                y_loc = y_loc - .05;
                uicontrol('Style','text','String',['type 1 sparsity = ' num2str(parameter.type_conn)],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                y_loc = y_loc - .05;
                uicontrol('Style','text','String',['type 2 sparsity = ' num2str(parameter.type_2_conn)],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                y_loc = y_loc - .05;
                uicontrol('Style','text','String',['type 1 overlap = ' num2str(parameter.type_1_overlap*parameter.similar_type_1_connections)],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                y_loc = y_loc - .05;
                uicontrol('Style','text','String',['type 2 overlap = ' num2str(parameter.type_2_overlap*parameter.similar_type_2_connections)],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                y_loc = y_loc - .05;
                uicontrol('Style','text','String',['token individuation threshold = ' num2str(parameter.token_individuation_threshold)],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                y_loc = y_loc - .05;
                
                uicontrol('Style','text','String',['length b = ' num2str(parameter.length_threshold)],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                y_loc = y_loc - .05;
                uicontrol('Style','text','String',['length m = ' num2str(parameter.length_threshold_variable)],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                y_loc = y_loc - .05;
                
                uicontrol('Style','text','String',['deviation b = ' num2str(parameter.deviation_threshold)],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                y_loc = y_loc - .05;
                uicontrol('Style','text','String',['deviation m = ' num2str(parameter.deviation_threshold_variable)],'Units','normalized','Position',[x_loc y_loc .25 .05]);
                y_loc = y_loc - .05;
                
            end
        end
    end
    
end

% GETTING RMSE

for task = whichtasks
    
    if task == 1
        bays_sd = [14.5 19.2 23.7 24.7];
        bays_g = [.02 .065 .22 .215];
        bays_b = [0  .04  .11   .31];
        
        count = 0;
        
        for i = cr_setsize
            
            count = count +1;
            
            cr_rmse_sd(count) = (bays_sd(count) - data.continuous.error(count))^2;
            cr_rmse_g(count) = (bays_g(count) - data.continuous.guess(count))^2;
            cr_rmse_b(count) = (bays_b(count) - data.continuous.swap(count))^2;
            
        end
        
        data.continuous.sd_rmse = sqrt(mean(cr_rmse_sd));
        data.continuous.g_rmse = sqrt(mean(cr_rmse_g));
        data.continuous.b_rmse = sqrt(mean(cr_rmse_b));
        
    elseif task == 2
        
        hit =                [.79   .63    .55    .5];
        false_alarm =        [.1    .22    .28     .34];
        
        prop_change = [.1 .2 .35 .58 .80 .89 .92 .93 .95 .95 .93; ...
            .22 .25 .42 .40 .54 .66 .75 .74 .81 .77 .82; ...
            .27 .27 .30 .40 .42 .52 .56 .64 .69 .72 .71; ...
            .33 .30 .36 .42 .40 .43 .60 .53 .52 .62 .63];
        
        
        
        bins = 11;
        count = 0;
        for i = cd_setsize
            
            count = count + 1;
            
            rmse_hit(count) = (hit(count) - (data.change_detection.hit(count)))^2;
            rmse_fa(count) = (false_alarm(count) - data.change_detection.false_alarm(count))^2;
            
            
            if range(parameter.changed_items) == 0
                for bins = 1:bins
                    rmse_prop_change(count,bins) = (prop_change(count,bins) - (data.change_detection.prop_change(count,bins)))^2;
                end
            end
            
            
        end
        
        data.change_detection.rmse_hit = sqrt(mean(rmse_hit));
        data.change_detection.rmse_fa = sqrt(mean(rmse_fa));
        
        data.change_detection.rmse_prop_change = sqrt(mean(mean(rmse_prop_change)));
    end
end

%
% if ispc|ismac
%
% else
%     exit
% end

data.parameters = parameter;

try
    data.cd = cd;
catch
    data.cd = [];
end

try
    data.cr = cr;
catch
    data.cr = [];
end

end
