function parameter = RMSE_and_Plot_Model_Sim(Model_sim,parameter,cd_sim)

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

%% This function calculates RMSE values and initiates plotting
% for MCMC related stuff.

if parameter.change_detection
    % Calculates RMSE for change detection simulations compared with
    % Keshvari et al (2013)
    
    Keshvari_hit =                [.79   .63    .55    .5];
    Keshvari_false_alarm =        [.1    .22    .28     .34];

    Keshvari_ss = [.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];

    count = 0;
    cr_rmse_cd = 0;
    for i = parameter.setsize
        if i == 2 || i == 4 || i == 6 || i == 8
            count = count +1;
            cr_rmse_hit(count) = (Keshvari_hit(count) - (1-Model_sim(i).miss))^2;
            cr_rmse_fa(count) = (Keshvari_false_alarm(count) - Model_sim(i).false_alarm)^2;

            if max(parameter.changed_items) == 1

                for bins = 1:cd_sim.bins
                    cr_rmse_cd(count,bins) = (Keshvari_ss(count,bins) - (cd_sim.proportions_of_change(count,bins)))^2;
                end
            end
        end
    end

    parameter.hit_rmse = sqrt(mean(cr_rmse_hit));
    parameter.fa_rmse = sqrt(mean(cr_rmse_fa));
    parameter.cd_rmse = sqrt(mean(mean(cr_rmse_cd)));

    if parameter.plot_Model_sim
        Plot_Model_sim(Model_sim,parameter,cd_sim)
    end

elseif parameter.continuous_report
    % Calculates RMSE values from continuous report from Bays et al 2009
    
    if parameter.confidence == 0
        % bays 2009 stimulus duration = 2000ms
        bays_sd =       [13   16.25    21.2    23];
        bays_g =        [.03    .04    .08     .07];
        bays_b =        [0      .02    .11     .28];

        % bays 2009 stimulus duration = 100ms
        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 = parameter.setsize
            if i == 1 || i == 2 || i == 4 || i == 6
                count = count +1;

                cr_rmse_sd(count) = (bays_sd(count) - Model_sim(i).deviation)^2;
                cr_rmse_g(count) = (bays_g(count) - Model_sim(i).prob_guessing)^2;
                cr_rmse_b(count) = (bays_b(count) - Model_sim(i).prob_distractor)^2;

            end
        end

        parameter.sd_rmse = sqrt(mean(cr_rmse_sd));
        parameter.g_rmse = sqrt(mean(cr_rmse_g));
        parameter.b_rmse = sqrt(mean(cr_rmse_b));

    end
    
    if parameter.plot_Model_sim 
        Plot_Model_sim(Model_sim,parameter,cd_sim)
    end

end
