function [Model_sim cd_sim parameter] = Threshold_analysis(Model_sim,parameter)

%    This file is part of Threshold_analysis.
%
%    Threshold_analysis 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.
%
%    Threshold_analysis 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 Threshold_analysis.  If not, see
%    <http://www.gnu.org/licenses/>.

%% This function determines the trial categorization for change detection
% tasks
%

if nargin == 0
    
    load model_stuff_sanity_test
    parameter.MCMC = 1;
    
end

index = 0;

search_width = 0.005;
% search_width = .001 found good stuff

total_sims = 100;


checkpoint_trial =100;
variance_check = 50;
variance_threshold = .0001;  %Threshold of variance required to end the simulation
variance_threshold = 1;

jumpyness = 1;   %modulate the tendency to jump:  <1 = more conservative  > 1 = more liberal

stuck_threshold = 10;

variance_jump = 0;
threshold_jump = 1;

show_gui = 0;
no_plots = 0;

good_threshold = 0;
current_threshold = 0;
finish_sim = 0;


          bins = 11;
      Keshvari_binwidth = 2*((parameter.sensory_space/2)/(bins-1))-((parameter.sensory_space/2)/(bins-1));
                                      

for length_slope = parameter.length_threshold_variable
    for length_intercept = parameter.length_threshold
        for deviation_slope = parameter.deviation_threshold_variable
            for deviation_intercept = parameter.deviation_threshold
                while good_threshold == 0
                    
                    count = 0;
                    
                    index = index + 1;
                    
%                     if ispc || ismac
                        if mod(index,(total_sims*.05)) == 0
                           sprintf('%g out of %g possible chains',index,total_sims)
                        end
%                     end
%                     
                    for setsize = parameter.setsize
                        
                        count = count + 1;
                        
                        Model_sim(setsize).miss = 0;
                        
                        length_threshold = length_intercept+(length_slope*setsize);
                        deviation_threshold = deviation_intercept+(deviation_slope*setsize);
                        
              
                        changed = zeros(length(unique(parameter.setsize)),bins);
                        no_changed = zeros(length(unique(parameter.setsize)),bins);
                        
                        reported_change = [];
                        reported_no_change = [];
                        
                        
                        count_hit_stimulus = 0;
                        count_miss_stimulus = 0;
                        count_cr_stimulus = 0;
                        count_fa_stimulus = 0;
                        count_change_trial = 0;
                        count_nochange_trial = 0;
                        
                        for trial = 1:parameter.trialcount
                            
                            hitormiss = zeros(1,parameter.probe_items(setsize==parameter.setsize));
                            crorfa = zeros(1,parameter.probe_items(setsize==parameter.setsize));
                            % This will need to be fixed to do single display
                            % something up here that %50 of the time, selects a
                            % changed item given the number of display size
                            
                            changetrial = 0;
                            
                            % If there are fewer display items than the
                            % setsize, then we have to select which items to
                            % probe.
                            if parameter.probe_items(setsize==parameter.setsize) ~= setsize
                                
                                check = 0;
                                for item = 1:parameter.probe_items(setsize==parameter.setsize)
                                    
                                    % If this is a change trial, we need to
                                    % make sure we pick an item that actually
                                    % changed. Currently, this will just pick
                                    % the first item, which may or may not be a
                                    % bug.
                                    if sum(Model_sim(setsize).changed_items(trial,:)) ~= 0 && check == 0
                                        
                                        new_item = Shuffle(Model_sim(setsize).changed_items(trial,:));
                                        display_items(item) = new_item(1);
                                        check = 1;
                                    elseif check == 1
                                        % If we've already grabbed a new item,
                                        % let's grab an old item!
                                        old_item = 1:setsize;
                                        display_items(item) = old_item(old_item ~= display_item(item));
                                        
                                    else
                                        % If this isn't a change trial, let's
                                        % just grab regular items.
                                        display_items(item) = item;
                                    end
                                    
                                end
                            else
                                display_items = 1:setsize;
                            end
                            
                            for token = display_items
                                try
                                    if max(Model_sim(setsize).changed_items(trial) == token)
                                        % changed stimulus
                                        changetrial = 1;
                                        
                                        if  Model_sim(setsize).deviation(trial,token) > deviation_threshold && Model_sim(setsize).retrieved_length(trial,token) > length_threshold
                                            % hit stimulus
                                            
                                            hitormiss(token) = 1;
                                            
                                            count_hit_stimulus = count_hit_stimulus + 1;
                                            
                                            if parameter.scatter
                                            Model_sim(setsize).hit_stimulus_deviation(count_hit_stimulus) = Model_sim(setsize).deviation(trial,token);
                                            Model_sim(setsize).hit_stimulus_retrieved_length(count_hit_stimulus) = Model_sim(setsize).retrieved_length(trial,token);
                                            Model_sim(setsize).hit_stimulus_old_to_new(count_hit_stimulus) = Model_sim(setsize).old_to_new(trial,token);
                                            end
                                            
                                        else
                                            % miss stimulus
                                            
                                            hitormiss(token) = -1;
                                            
                                            
                                            count_miss_stimulus = count_miss_stimulus + 1;
                                            
                                            if parameter.scatter
                                            Model_sim(setsize).miss_stimulus_deviation(count_miss_stimulus) = Model_sim(setsize).deviation(trial,token);
                                            Model_sim(setsize).miss_stimulus_retrieved_length(count_miss_stimulus) = Model_sim(setsize).retrieved_length(trial,token);
                                            Model_sim(setsize).miss_stimulus_old_to_new(count_miss_stimulus) = Model_sim(setsize).old_to_new(trial,token);
                                            end
                                        end
                                        
                                    else
                                        % non changed stimulus
                                        
                                        if Model_sim(setsize).deviation(trial,token) > deviation_threshold && Model_sim(setsize).retrieved_length(trial,token) > length_threshold
                                            % False alarm stimulus
                                            crorfa(token) = -1;
                                            
                                            
                                            count_fa_stimulus = count_fa_stimulus + 1;
                                            if parameter.scatter
                                            Model_sim(setsize).fa_stimulus_deviation(count_fa_stimulus) = Model_sim(setsize).deviation(trial,token);
                                            Model_sim(setsize).fa_stimulus_retrieved_length(count_fa_stimulus) = Model_sim(setsize).retrieved_length(trial,token);
                                            end
                                        else
                                            % CR stimulus
                                            crorfa(token) = 1;
                                            count_cr_stimulus = count_cr_stimulus + 1;
                                            if parameter.scatter
                                            Model_sim(setsize).cr_stimulus_deviation(count_cr_stimulus) = Model_sim(setsize).deviation(trial,token);
                                            Model_sim(setsize).cr_stimulus_retrieved_length(count_cr_stimulus) = Model_sim(setsize).retrieved_length(trial,token);
                                            end
                                        end
                                    end
                                    
                                catch
                                    keyboard
                                end
                                
                            end
                            
                            
                            if max((hitormiss ==1))    %this is a hit trial
                                false_alarm_trial(trial) = 0;
                                miss_trial(trial) = 0;
                                
                                count_change_trial = count_change_trial +1;
                                reported_change(count_change_trial,index) = Model_sim(setsize).old_to_new(trial,Model_sim(setsize).changed_items(trial));
                                
                            elseif max(crorfa == -1) && changetrial  %his is a hit trial because of a false alarm stimulus
                                false_alarm_trial(trial) = 0;
                                miss_trial(trial) = 0;
                                
                                count_change_trial = count_change_trial +1;
                                reported_change(count_change_trial,index) = Model_sim(setsize).old_to_new(trial,Model_sim(setsize).changed_items(trial));
                            elseif max((hitormiss == -1)) %this is a miss trial
                                miss_trial(trial) = 1;
                                false_alarm_trial(trial) = 0;
                                
                                count_nochange_trial = count_nochange_trial + 1;
                                reported_no_change(count_nochange_trial,index) =  Model_sim(setsize).old_to_new(trial,Model_sim(setsize).changed_items(trial));
                            elseif max(crorfa == -1) && changetrial ==0  %his is a false alarm trial
                                false_alarm_trial(trial) = 1;
                                miss_trial(trial) = 0;
                                
                                count_change_trial = count_change_trial + 1;
                                reported_change(count_change_trial,index) = 0;
                                
                            else  %this is a CR trial
                                false_alarm_trial(trial) = 0;
                                miss_trial(trial) = 0;
                                
                                count_nochange_trial = count_nochange_trial + 1;
                                reported_no_change(count_nochange_trial,index) = 0;
                                
                            end
                            
                        end
                        
                        
                        % Plots scatter plots of length compared to deviation
                        if parameter.scatter
                            
                            figure('Color',[1 1 1])
                            subplot(1,2,1);
                            hold on
                            if count_fa_stimulus ~= 0
                                scatter(Model_sim(setsize).fa_stimulus_deviation,Model_sim(setsize).fa_stimulus_retrieved_length,40,'MarkerEdgeColor','r')
                            end
                            if count_cr_stimulus ~= 0
                                scatter(Model_sim(setsize).cr_stimulus_deviation,Model_sim(setsize).cr_stimulus_retrieved_length,40,'MarkerEdgeColor','g')
                            end
                            
                            title(num2str(setsize),'FontSize',16);
                            title('');
                            xlim([0 1]);
                            ylim([-.01 1]);
                            
                            set(gca,'Color',[1 1 1]);
                            box off
                            xlims = get(gca,'XLim');
                            ylims = get(gca,'YLim');
                            line(xlims,[ylims(2) ylims(2)],'Color','black');
                            line([xlims(2) xlims(2)],ylims,'Color','black');
                            set(gca,'Xtick',[0 .5 1],'FontSize',22,'FontName','Arial');
                            
                            
                            subplot(1,2,2);
                            hold on
                            
                            if count_hit_stimulus ~= 0
                                for(i = 1:length(Model_sim(setsize).hit_stimulus_deviation))
                                    
                                    colorhit = Model_sim(setsize).hit_stimulus_old_to_new(i);
                                    scatter(Model_sim(setsize).hit_stimulus_deviation(i),Model_sim(setsize).hit_stimulus_retrieved_length(i),40,'MarkerEdgeColor','g')
                                    
                                end
                            end
                            
                            if count_miss_stimulus ~= 0
                                for(i = 1:length(Model_sim(setsize).miss_stimulus_deviation))
                                    
                                    colormiss = Model_sim(setsize).miss_stimulus_old_to_new(i);
                                    scatter(Model_sim(setsize).miss_stimulus_deviation(i),Model_sim(setsize).miss_stimulus_retrieved_length(i),40,'MarkerEdgeColor','r')
                                    
                                end
                            end
                            
                            title(num2str(setsize),'FontSize',16)
                            title('');
                            xlim([0 1]);
                            ylim([-.01 1]);
                            
                            set(gca,'Color',[1 1 1]);
                            box off
                            xlims = get(gca,'XLim');
                            ylims = get(gca,'YLim');
                            line(xlims,[ylims(2) ylims(2)],'Color','black');
                            line([xlims(2) xlims(2)],ylims,'Color','black');
                            set(gca,'Xtick',[0 .5 1],'FontSize',22,'FontName','Arial');
                            
                            
                        end
                        
                        m(setsize,index) = sum(miss_trial)/(parameter.trialcount/2);
                        fa(setsize,index) = sum(false_alarm_trial)/(parameter.trialcount/2);
                        
                        
                        
                        % Plots the proportion of responding 'change' given a deviation
                        
                        
                            if range(parameter.changed_items) == 0
                                if max(parameter.changed_items)== 1
                                    
                                    if length(reported_change) > 1
                                        for k = 1:length(reported_change(:,index))
                                            
                                            bin = ceil(reported_change(k,index)/  Keshvari_binwidth)+1;
                                            changed(count,bin) = changed(count,bin) + 1;

                                        end
                                    end
                                    
                                    if length(reported_no_change) > 1
                                        for k = 1:length(reported_no_change(:,index))
                                            
                                            bin = ceil(reported_no_change(k,index)/  Keshvari_binwidth)+1;
                                            no_changed(count,bin) = no_changed(count,bin) + 1;

                                        end
                                    end
                                    
                                    proportions_of_change(count,:,index) = (changed(count,:)./(changed(count,:)+no_changed(count,:)));
                                    
                                else
                                    
                                    proportions_of_change(count,:,index) = 0;
                                end
                            end
                        
                    end
                    
                    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;
                    rmse_count = 0;
                    for i = parameter.setsize
                        if i == 2 || i == 4 || i == 6 || i == 8
                            rmse_count = rmse_count + 1;
                            rmse_hit(rmse_count,index) = (hit(rmse_count) - (1-m(i,index)))^2;
                            rmse_fa(rmse_count,index) = (false_alarm(rmse_count) - fa(i,index))^2;
                            
                            
                            if range(parameter.changed_items) == 0
                                for bins = 1:bins
                                    rmse_prop_change(rmse_count,bins,index) = (prop_change(rmse_count,bins) - (proportions_of_change(rmse_count,bins,index)))^2;
                                end
                            end
                        end
                    end
                    
                    rmse_hit = sqrt(mean(rmse_hit(:,index)));
                    rmse_fa = sqrt(mean(rmse_fa(:,index)));
                    
                    rmse_prop_change = sqrt(mean(mean(rmse_prop_change(:,:,index))));
                    rmse_thresholds(:,index) = [length_intercept length_slope deviation_intercept deviation_slope];
                    
                    if(parameter.fixed_variance)  %ignore Keshvari fit during RMSE because we're not using all the data
                        rmse_mean(index) = mean([rmse_hit rmse_fa ]);
                    else
                        rmse_mean(index) = mean([rmse_hit rmse_fa rmse_prop_change]);
                    end
                    
                    if index > checkpoint_trial
                        rmse_variance = var(rmse_mean((index-variance_check):index));
                    else
                        rmse_variance = 1;
                    end
                    
                    if parameter.MCMC == 0 || index > total_sims || rmse_variance < variance_threshold || finish_sim
                        
                        good_threshold = 1;
                        [best_rmse, best_index] = sort(rmse_mean);
                        
                    elseif parameter.MCMC
                        
                        jumped(index) = 0;
                        prob_of_staying(index) = 0;
                        
                        if index == 2
                            comparison_index = index-1;
                        end
                        
                        
                        if variance_jump
                            if index > 1
                                
                                prob_of_staying(index) = (rmse_mean(index)/(rmse_mean(index)+rmse_mean(comparison_index)))^jumpyness;
                                
                                if prob_of_staying(index) > rand
                                    new_var_value = var_values(1);
                                else
                                    new_var_value = var_values(2);
                                    jumped(index) = 1;
                                    comparison_index = index;
                                end
                            end
                        elseif threshold_jump
                            if index > 1
                                if rmse_mean(index) < rmse_mean(comparison_index)
                                    new_var_value = var_values(2);
                                    jumped(index) = 1;
                                    comparison_index = index;
                                else
                                    new_var_value = var_values(1);
                                end
                            end
                        end
                        
                        if index > 1
                            switch whichvar
                                case 1
                                    length_slope = new_var_value;
                                case 2
                                    length_intercept = new_var_value;
                                case 3
                                    deviation_slope = new_var_value;
                                case 4
                                    deviation_intercept = new_var_value;
                            end
                        end
                        
                        whichvar = randi(4);
                        
                        %                         if index > stuck_threshold
                        %
                        %                             if mod(index,stuck_threshold) ==1
                        %
                        % %                                 jumped_rmse = rmse_mean(jumped == 1);
                        % %                                 size_jumped_rmse = length(jumped_rmse);
                        % %                                 recent_rmse = [jumped_rmse(size_jumped_rmse - 1) jumped_rmse(size_jumped_rmse)];
                        % %
                        %                                 if (sum(jumped((index-stuck_threshold):index)) == 0)
                        %                                     search_width = search_width * 1.5;
                        %                                 else
                        %                                     search_width = search_width *.5;
                        %                                 end
                        %
                        %                             end
                        %                         end
                        
                        if search_width > .5
                            keyboard
                        end
                        
                        
                        switch whichvar
                            case 1
                                var_values(1) = length_slope;
                                var_values(2) = normrnd(var_values(1),search_width);
                                length_slope = var_values(2);
                            case 2
                                var_values(1) = length_intercept;
                                var_values(2) = normrnd(var_values(1),search_width*5);
                                length_intercept = var_values(2);
                            case 3
                                var_values(1) = deviation_slope;
                                var_values(2) = normrnd(var_values(1),search_width);
                                deviation_slope = var_values(2);
                            case 4
                                var_values(1) = deviation_intercept;
                                var_values(2) = normrnd(var_values(1),search_width*5);
                                deviation_intercept = var_values(2);
                        end
                        
                        if no_plots
                            if show_gui
                                output = MCMC_tracking(rmse_mean(index),index,total_sims,jumped(index),search_width,1);
                                
                                if output.abort_sim
                                    keyboard
                                elseif output.finish_sim
                                    finish_sim = 1;
                                end
                            else
                                
                                if index == 1
                                    figure
                                end
                                hold on
                                
                                color = [0 1 0];
                                if jumped(index)
                                    color = [1 0 0];
                                end
                                
                                scatter(index,rmse_mean(index),20,color);
                                
                                pause(.005);
                                
                            end
                        end
                    end
                end
            end
        end
    end
end

index = best_index(1);
count = 0;

for setsize = parameter.setsize
    count = count + 1;
    
    Model_sim(setsize).miss = m(setsize,index);
    Model_sim(setsize).false_alarm = fa(setsize,index);
    %Model_sim(setsize).reported_change = reported_change([2 4 6 8],index);
    %Model_sim(setsize).reported_no_change = reported_no_change(:,index);
    parameter.length_threshold = rmse_thresholds(1,index);
    parameter.length_threshold_variable = rmse_thresholds(2,index);
    parameter.deviation_threshold = rmse_thresholds(3,index);
    parameter.deviation_threshold_variable = rmse_thresholds(4,index);
    
    cd_sim.proportions_of_change(count,:) = proportions_of_change(count,:,index);
    cd_sim.bins = bins;
end

if parameter.MCMC
    RMSE_and_Plot_Model_Sim(Model_sim,parameter,cd_sim);
end

end

