classdef Model < handle
    properties
        
        scmodel;        % SparseCoding class
        rlmodel_blur;   % ReinforcementLearning class
        rlmodel_disp;   % ReinforcementLearning class
        Act_blur;       % action range for blur agent
        Act_disp;       % action space for diparity agent

        % simulation parameters
        T_train;            % total number of interations in the training
        Interval;           % period to change a new environment for the eye
        configpath;         % configuration file path
        trainedUntil;       % in case the training was disrupted, how long was the model trained?
        savePath;           % results file path
        maxObjPlane;        % maximal distance of the object plane
        maxAccDispPlane;    % maximal distance of the accomodation plane
        maxDisp;            % maximally possible disparity
        maxBlur;            % maximally possible blur
        saveAt;             % time points where the model data is saved
        nscales;            % number of scales, i.e., number of sparse coders
        useSuppression;     % bool. number indicating the use of suppressive mechanisms
        threshold;          % threshold of contrast units
        slope;              % slope of contrast units
        eta_supp;           % update rate of contrast estimates
        suppr;              % suppressive effect of suppression model
        exci;               % excitatory effect of suppression model
        simulationTime;     % number of hours, the model was trained
        window_size;        % size of initial window cutout before disparities are applied
        trainRL;            % training of reinforcement learning
        trainSC;            % training of sparse coder
        n_textures;         % number of input images used during training
        texture_directory;  % where the input images are stored
        image_size;         % size of input images
        colored_input_images;% gray values or rgb

        %% record history of values
        recerr_hist;        % history of rec error
        reward_hist;        % history of reward
        vergerr_hist;       % history of vergence error
        accerr_hist;        % history of accommodation error
        vergPlane_hist;     % history of vergence plane
        accPlane_hist;      % history of accomodation plane
        objPlane_hist;      % history of object position
        texture_hist;       % history of which texture was shown
        feature_hist;       % history of base activation during learning
        monocularity_hist;  % history of monocularity of chosen basis functions
        contrast_hist;      % history of contrast adjustments to right and left eye
        If_avg_right_hist   % history of average freq spectrum right
        If_avg_left_hist    % history of average freq spectrum left
        distances;          %object distances during simulation

        %% model data
        blur_reward_mean;    % moving average of reward before normalization
        blur_reward_variance;% moving average of variance of reward before normalization
        disp_reward_mean;    % moving average of reward before normalization
        disp_reward_variance;% moving average of variance of reward before normalization
        alpha;               % update rate for reward moving averages
        blur_scale;          % set relation between 1 [a.u.] of accommodation error and stdv of defocus blur filter

        % simulate visual impairment (anisometropy + cataract)
        spectacles_l;
        spectacles_r;
        cataract_l;
        cataract_r;
    end

    methods
        function obj = Model(PARAM)
            
            obj.nscales = PARAM{2}{1}{7};
            
            obj.configpath = PARAM{1}{1};
            obj.trainSC = PARAM{1}{17};
            obj.trainRL = PARAM{1}{9};
            
            obj.T_train = PARAM{1}{2};
            obj.Interval = PARAM{1}{3};
            obj.trainedUntil = 0;
            obj.savePath = PARAM{1}{8};
            obj.maxObjPlane = PARAM{1}{22};
            obj.maxAccDispPlane = PARAM{1}{23};
            obj.maxDisp = PARAM{1}{10};
            obj.maxBlur = PARAM{1}{11};
            obj.saveAt = PARAM{1}{12};
            
            obj.n_textures = PARAM{1}{18};
            obj.texture_directory = PARAM{1}{19};
            obj.image_size = PARAM{1}{20};
            obj.window_size = PARAM{1}{14};
            obj.colored_input_images = PARAM{1}{21};
            obj.blur_scale = PARAM{1}{16};
            
            obj.useSuppression = PARAM{1}{13}{1};
            obj.threshold = PARAM{1}{13}{2};
            obj.slope = PARAM{1}{13}{3};
            obj.eta_supp = PARAM{1}{13}{4};
            obj.suppr = PARAM{1}{13}{5};
            obj.exci = PARAM{1}{13}{6};
            
            obj.blur_reward_mean=0;
            obj.disp_reward_mean=0;
            obj.blur_reward_variance=1;
            obj.disp_reward_variance=1;
            obj.alpha = PARAM{1}{15};
            
            obj.recerr_hist = zeros(obj.T_train, 1);          % coarse and fine scale reconstruction error
            obj.reward_hist = zeros(obj.T_train, 7);          % coarse and fine scale reward
            obj.vergerr_hist = zeros(obj.T_train, 1);         % saved vergence values
            obj.accerr_hist = zeros(obj.T_train, 1);          % saved blur values
            obj.vergPlane_hist = zeros(obj.T_train, 1);       % save history of verg. plane positions
            obj.accPlane_hist = zeros(obj.T_train, 1);        % save history of acc. plane positions    
            obj.objPlane_hist = zeros(obj.T_train, 1);        % save history of obj. plane positions
            obj.texture_hist = zeros(obj.T_train, 1);         % texture history
            obj.monocularity_hist = zeros(obj.T_train, 5);    % history of monocularity of chosen basis
            obj.contrast_hist = ones(obj.T_train, 2);         % left and right contrast adjustments
            obj.distances = zeros(1,PARAM{1}{2}/PARAM{1}{3}); % holds object distances

            obj.spectacles_l = PARAM{1}{4};
            obj.spectacles_r = PARAM{1}{5};
            obj.cataract_l = PARAM{1}{6};
            obj.cataract_r = PARAM{1}{7};

            %init sparse coders
            for s = 1:obj.nscales
                obj.scmodel{s} = SparseCoding(PARAM{2}{s}, size(PARAM{1}{12}, 2));
            end

            %init RLs
            obj.rlmodel_blur = ReinforcementLearning(PARAM{3}, PARAM{1}{2}, size(PARAM{1}{12}, 2));
            obj.rlmodel_disp = ReinforcementLearning(PARAM{4}, PARAM{1}{2}, size(PARAM{1}{12}, 2));
            obj.Act_blur = PARAM{3}{11};
            obj.Act_disp = PARAM{4}{11};
        end
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%
        %%% Generate the feature vector, rewards and encoding errors.
        %%% Also returns current monocularity of receptive fields and coef
        %%% of cortical neurons after encoding.
        %%% 'feature' is the concatenation of foveal and peripheral scale
        %%% 'reward' is the sum of rewards for foveal and peripheral scale
        %%% 'error' is the sum of errors of foveal and peripheral scale
        %%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function [feature, reward, err, coef_sparse, monocularity] = generateFR(this, patchesL, patchesR, t, blurRewL, blurRewR)
            
            if length(this.scmodel) == 1
                
                [coef_sparse, err, monocularity] = this.scmodel{1}.sparseEncode([patchesL{1}; patchesR{1}]);
                feature = mean(coef_sparse.^2, 2);
                feature = feature./sum(feature);
                coef_sparse = {coef_sparse};
                rew = {mean(sum(err.^2))};
                
            elseif length(this.scmodel) == 2
                
                % peripheral scale
                [coef_coarse, error_coarse, monocularityC] = this.scmodel{1}.sparseEncode([patchesL{1}; patchesR{1}]);
                featureCoarse = mean(coef_coarse.^2, 2);
                featureCoarse = featureCoarse./sum(featureCoarse);
                rewCoarse = -mean(sum(error_coarse.^2));
                
                % foveal scale
                [coef_fine, error_fine, monocularityF] = this.scmodel{2}.sparseEncode([patchesL{2}; patchesR{2}]);
                featureFine = mean(coef_fine.^2, 2);
                featureFine = featureFine./sum(featureFine);
                rewFine = -mean(sum(error_fine.^2));

                err = {error_coarse; error_fine};
                rew = rewCoarse + rewFine;
                feature = [featureCoarse; featureFine];
                monocularity = [monocularityC; monocularityF];
                coef_sparse = {coef_coarse; coef_fine};
            end
    
            % combine rewards and divide by 100 in order to get small
            % reward at beginning of simulation that grows to unit variance
            % due to online normalization. This prevents instabilities due
            % to random initialization.
            disp_reward = rew/100;                 
            blur_reward = (blurRewL + blurRewR)/100;

            % zero mean unit variance reward via exponentially weighted running average
            r = (blur_reward-this.blur_reward_mean)/(sqrt(this.blur_reward_variance) + eps);
            
            delta=blur_reward-this.blur_reward_mean;
            this.blur_reward_mean = this.blur_reward_mean+this.alpha*delta;
            this.blur_reward_variance = (1-this.alpha)*this.blur_reward_variance+this.alpha*delta^2;

            blur_reward = r;

            r = (disp_reward-this.disp_reward_mean)/(sqrt(this.disp_reward_variance) + eps);
            delta=disp_reward-this.disp_reward_mean;
            this.disp_reward_mean = this.disp_reward_mean+this.alpha*delta;
            this.disp_reward_variance = (1-this.alpha)*this.disp_reward_variance+this.alpha*delta^2;
            
            disp_reward = r;

            reward=[blur_reward disp_reward rew];

        end

        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%
        %%% Plot results 
        %%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %% Show receptive fields
        function h = displayBasis(this, scle, M, N, saveTime, savePlot)

            s=sqrt(this.scmodel{scle}.Basis_size/2);
            n = this.scmodel{scle}.Basis_num;

            if isempty(saveTime)
                block = mat2gray(col2im(this.scmodel{scle}.Basis,[s 2*s],[n*s 2*s],'distinct'));
            else
                block = mat2gray(col2im(this.scmodel{scle}.Basis_hist(:, :, saveTime),[s 2*s],[n*s 2*s],'distinct'));
            end
            k=1;
            pic=[];

            for j=1:N
                row=[];
                for i=1:M
                    base_l=block(1+(k-1)*s:k*s,1:s);
                    base_r=block(1+(k-1)*s:k*s,s+1:2*s);
                    brick=[-ones(1,s); base_l; -ones(1,s); base_r; -ones(1,s)]; %left part on top
                    brick=[-ones(2*s+3,1) brick -ones(2*s+3,1)];
                    row=[row brick];
                    k=k+1;
                end

                pic=[pic;row];
            end
            
            figure;
            h = imshow(pic);
            if savePlot
                if scle == 1
                    tag = '_coarse';
                elseif scle == 2
                    tag = '_fine';
                end
                saveas(gca, strcat(this.savePath, 'basis', tag), 'png');
            end

        end

        %% Display the N x M Basis that where choosen for encoding in the last step
        function h = displaySelectedBasis(this, scle, N, M, saveName)

            s = sqrt(this.scmodel{scle}.Basis_size/2);    % single patch size
            n = this.scmodel{scle}.Basis_num;             % total number of basis
            sBox = 2;                                     % size of black box surrounding the basis

            usedB = zeros(this.scmodel{scle}.Basis_num, 1);
            usedB(this.scmodel{scle}.Basis_selected > 0) = 1; % extract indices of basis that where selected by matching pursuit
            inds = find(usedB);
            if isempty(inds)
                sprintf('model contains no selected basis')
                return;
            end

            b = 1;
            h = figure;
            hold on;
            for i = 1:M
                for j = 1:N
                    subplot(M, N, b)

                    img = zeros(2*s + 2*sBox, s + 2*sBox)';    % create the black borders
                    img(sBox+1 : end-sBox, sBox+1 : end-sBox) = mat2gray(reshape(this.scmodel{scle}.Basis(:,inds(b)), [sqrt(size(this.scmodel{1}.Basis,1)/2), sqrt(size(this.scmodel{1}.Basis,1)/2)*2]));
                    imshow(img)
                    title(sprintf('BI = %0.2f', calculateRightBinocularity(this.scmodel{scle}.Basis(:,inds(b)))));
                    b = b + 1;
                end
            end

            if ~isempty(saveName)
                if scle == 1
                    tag = '_coarse';
                elseif scle == 2
                    tag = '_fine';
                end
                saveas(h, strcat(this.savePath, 'selectedBasis', saveName, tag), 'png');
            end
        end

        %% Display 'nSamples' sample patches with the basis used for encoding
        function h = displayEncoding(this, scle, texture, disparity, blurLeft, blurRight, contrastLeft, contrastRight, nSamples, saveName)
            s = sqrt(this.scmodel{scle}.Basis_size/2);  % single patch size
            sBox = 2;                                   % size of black box surrounding the basis
            suppPostNorm = 0;                           % normalization after patch adaptation?
            blurScale = 0.8;

            params={s, this.scmodel{scle}.Dsratio, blurScale, this.cataract_l, this.cataract_r};
            img = imread(sprintf('./textures/%d.bmp', texture)); %example image

            % generate patches (zero mean, unit norm)
            [patchesLeft, patchesRight] = PatchGenerator(img, disparity, blurLeft-this.spectacles_l, blurRight-this.spectacles_r, this, params);

            % contrast adaptation
            newPatchesLeft = patchesLeft .* contrastLeft;
            newPatchesRight = patchesRight .* contrastRight;

            if suppPostNorm
                newPatchesLeft = newPatchesLeft - repmat(mean(newPatchesLeft),[size(newPatchesLeft,1) 1]);
                newPatchesRight = newPatchesRight - repmat(mean(newPatchesRight),[size(newPatchesRight,1) 1]);
            end
            
            [coefs, err, monocularity] = this.scmodel{scle}.suppressiveEncode([newPatchesLeft; newPatchesRight]);

            h = figure;
            hold on;
            rng(1)
            for p = 1 : nSamples

                pIndex = randi(size(patchesLeft, 2)); % take one random patch
                
                % display original patch
                patch = [patchesLeft(:, pIndex), patchesRight(:, pIndex)];
                subplot(nSamples*2, 10, [1 + (20*(p-1)), 11 + (20*(p-1))]);
                imshow(mat2gray(col2im(patch, [s s], [2*s s], 'distinct')))
                title('orig');

                % display altered patch
                patchAdap = [newPatchesLeft(:, pIndex), newPatchesRight(:, pIndex)];
                subplot(nSamples*2, 10, [2 + (20*(p-1)), 12 + (20*(p-1))]);
                imshow(mat2gray(col2im(patchAdap, [s s], [2*s s], 'distinct')))
                title('adapt');

                % display reconstruction
                encPatch = this.scmodel{scle}.Basis*coefs(:, pIndex);
                recIm = [encPatch(1 : end/2), encPatch(end/2+1 : end)];
                subplot(nSamples*2, 10, [3 + (20*(p-1)), 13 + (20*(p-1))]);
                imshow(mat2gray(col2im(recIm, [s s], [2*s s], 'distinct')))
                title('rec.');

                % display selected basis
                binIndx = calculateRightBinocularity(this.scmodel{scle}.Basis);
                corrs = calculateCorrelation(this.scmodel{scle}.Basis);
                [~, bInds]=sort(abs(coefs(:,pIndex)));
                bInds = flip(bInds(end-this.scmodel{scle}.Basis_num_used+1: end));
                for b = 1 : this.scmodel{scle}.Basis_num_used
                    base = this.scmodel{scle}.Basis(:, bInds(b));
                    baseIm = [base(1 : end/2), base(end/2 + 1 : end)];
                    if b <= (this.scmodel{scle}.Basis_num_used / 2)
                        subplot(nSamples*2, 10, (20*(p-1)) + 3 + b);
                    else
                        subplot(nSamples*2, 10, (20*(p-1)) + 8 + b);
                    end
                    imshow(mat2gray(col2im(baseIm, [s s], [2*s s], 'distinct')))
                    text(1,s+1,sprintf('act:%0.2f\nbin:%0.2f\ncorr:%0.2f', coefs(bInds(b), pIndex), binIndx(1, bInds(b)), corrs(1, bInds(b))), 'color', 'y', 'fontsize', 6);
                end

                % display distribution of binocularity
                contribs = binIndx(coefs(:, pIndex) ~= 0);
                sp = subplot(nSamples*2, 10, [(20*(p-1)) + 9, (20*(p-1)) + 10, (20*(p-1)) + 19, (20*(p-1)) + 20]);
                histogram(contribs, linspace(0,1,10));
                title(sprintf('mean BI: %0.2f', mean(contribs)));
                set(sp, 'fontsize', 7);
            end

            i = regexp(this.savePath, '/');
            expName = this.savePath(i(end-1)+1:end-1);  % set the experiment name and params as title above all
            suptitle(sprintf('Exp: %s\ndisp=%1d, blurL=%d, blurR=%d, contrL=%.1f, contrR=%.1f',...
                strrep(expName, '_', '\_'), disparity, blurLeft-this.spectacles_l, blurRight-this.spectacles_r, contrastLeft, contrastRight));

            if ~isempty(saveName)
                if scle == 1
                    tag = '_coarse';
                elseif scle == 2
                    tag = '_fine';
                end
                saveas(h, sprintf('%s/patchEncoding%s_%s.png', this.savePath, tag, saveName));
            end
        end

        %% Displays the correlation of activated basis-sub-parts over multiple disparities
        function displayCorrelation(this, scle, saveTag)

            disparityRange = [-3:1:3];
            blurLeft = 0;
            blurRight = 0;

            s = sqrt(this.scmodel{scle}.Basis_size/2);    % single patch size
            sBox = 2;                               % size of black box surrounding the basis

            params={s, this.scmodel{scle}.Dsratio, 0.625, this.cataract_l, this.cataract_r};
            img = imread(sprintf('/home/klimmasch/projects/AccomodationLearning_orig/textures_mcgill/%d.bmp', img)); %example image

            nBins = 100;
            histograms = zeros(nBins, length(disparityRange));

            for d = 1 : length(disparityRange)

                [patchesLeft, patchesRight] = PatchGenerator(img, disparityRange(d), blurLeft, blurRight, this, params);
                % imshow(mat2gray(col2im(patchesLeft, [10 10], [110 110], 'distinct')))
                [coefs, err, monocularity] = this.scmodel{scle}.suppressiveEncode([patchesLeft; patchesRight]);
                pIndex = 1;
                binIndx = calculateBinocularity(this.scmodel{scle}.Basis);
                corrs = calculateCorrelation(this.scmodel{scle}.Basis);
                [~, bInds]=sort(abs(coefs(:,pIndex)));
                bInds = flip(bInds(end-this.scmodel{scle}.Basis_num_used+1: end, :)); %for each patch, 10 indices of selected basis
                h = histogram(corrs(bInds(:)), nBins);
                histograms(:, d) = h.Values;
            end

            f = figure;
            hold on;
            surf(histograms);
            set(gca, 'XTickLabel', [-3:1:3]);
            set(gca, 'YTickLabel', [-1:0.4:1]);
            xlabel('Disparity');
            ylabel('Freq. of Correlation');
            view([24, 30]);

            if ~isempty(saveTag)
                if scle == 1
                    tag = '_coarse';
                elseif scle == 2
                    tag = '_fine';
                end
                saveas(h, sprintf('%s/correlation_scale$d_%s', scle, this.savePath, saveTag),'png');
            end
        end

        function displayBinHist(this, scle, nBins, saveName)
            binocs = calculateRightBinocularity(this.scmodel{scle}.Basis);
            h = figure;
            histogram(binocs,linspace(0,1,nBins));
            xlabel('Binocularity Index');
            ylabel('Frequency');
            title(sprintf('mean BI = %0.2f', mean(binocs)));
            if ~isempty(saveName)
                if scle == 1
                    tag = '_coarse';
                elseif scle == 2
                    tag = '_fine';
                end
                saveas(h, strcat(this.savePath, saveName, tag,'.png'));
            end
        end

        %% calculates the reward surface for a given set of basis functions over different blur levels
        % textures specifies an array of images over which the reward is averaged.
        function rewards = rewardSurface(this, scle, basisAt, textures, displ, saveName)
            blurScale = 0.8;

            contrasts = [0, 0.5, 1, 1.5, 2];

            contrastsR = flip(contrasts);
            blurLevels = -3:3;

            rewards = zeros(length(blurLevels), length(contrasts));   % reward coarse, fine
            % rewards = zeros(length(blurLevels), length(textures));  % for this version, show reward surface for each img
                                                                      % over blur levels and contrasts

            for s = 1:length(this.scmodel)
                  for b = 1:length(blurLevels)
                      for texture = 1:length(textures)
                            img = imread(sprintf('%s%d.bmp', this.texture_directory, textures(texture)));

                            for c = 1:length(contrasts)
                                  contrastLeft = contrasts(c);
                                  contrastRight = contrastsR(c);
                                  [~, ~, blurRewLeft, blurRewRight] = PatchGenerator(img, 0, blurLevels(b)+this.spectacles_l, blurLevels(b)+this.spectacles_r, this, {10, this.scmodel{scle}.Dsratio, blurScale, this.cataract_l, this.cataract_r});

                                  blurRewLeft = blurRewLeft * contrastLeft;
                                  blurRewRight = blurRewRight * contrastRight;

                                  rew = blurRewLeft + blurRewRight;
                                  rewards(b, c) = rew + rewards(b, c);
                            end
                      end
                  end
            end

            rewards = rewards ./ length(textures);

            if displ
                h = figure('Position', [500 800 600 600]);
                surf(rewards(:, :));
                axis tight;

                s1 = gca();
                s1.XTick = 1:length(contrasts);
                s1.XTickLabel = contrasts;
                s1.YTick = 1:length(blurLevels);
                s1.YTickLabel = blurLevels;

                xlabel('contrast left');
                ylabel('blur error left');
                title('coarse scale');

                if ~isempty(saveName)
                    if scale == 1
                        tag = '_coarse';
                    elseif scale == 2
                        tag = '_fine';
                    end
                    saveas(gca, sprintf('%srewardSurf_%sScale_%s.png', this.savePath, tag, saveName));
                end
            end
        end

        %% Plots results from one simulation run
        function allPlotSave(this)
            
            % predefining some colors
            orange=[255 153 24]./255;
            blue=[39 48 73]./255;
            green=[146 208 80]./255;
            lightblue=[70 166 199]./255;
            purple=[139 84 166]./255;
            brown=[167 121 29]./255;
            gray=[100 100 100]/255;

            for s = 1:size(this.scmodel,2)
                this.displayBasis(s, 20, 15, [], 1); % display and save basis
                try
                    this.displaySelectedBasis(s, 4, 4, '_');
                catch
                    sprintf('It seems no basis where selected that could be displayed.')
                end

                this.displayEncoding(s, 66, 0, 0, 0, this.contrast_hist(this.trainedUntil, 1), this.contrast_hist(this.trainedUntil, 2), 2, 'example');
            end

            if this.useSuppression
                h = figure;
                windowsize = ceil(this.T_train/100);
                conts = filter(ones(1, windowsize) / windowsize, 1, this.contrast_hist);%(windowsize+1:end));
                plot(conts);
                ylim([0, 2]);
                title('Left and Right Contrasts over Train Time');
                xlabel('Time');
                ylabel('Contrasts');
                legend('left eye, coarse', 'right eye, coarse', 'left eye, fine', 'right eye, fine');
                saveas(gca, strcat(this.savePath, 'conts_hist.png'));
                close(h);
            end


            %% blur level
            windowSize = ceil(this.T_train/500);
            ind = 1:this.T_train;
            i = mod(ind,this.Interval);
            ind = find(i); %to exclude trials where vergence is reset
            ind=ind((this.Interval-1):(this.Interval-1):end);

            blur_r = filter(ones(1,windowSize)/windowSize,1,abs(this.reward_hist(ind,2)+this.spectacles_r));
            blur_r = blur_r(windowSize:end); %discard first values (mean not reliable)
            blur_l = filter(ones(1,windowSize)/windowSize,1,abs(this.reward_hist(ind,2)+this.spectacles_l));
            blur_l = blur_l(windowSize:end); %discard first values (mean not reliable)
            ind = ind(windowSize:end);

            h=figure; hold on;
            plot(ind,blur_r,'-','LineWidth',1,'Color',green); xlabel('Time','FontSize',8); ylabel('Stdv [px]','FontSize',8);
            plot(ind,blur_l,'-','LineWidth',1,'Color',lightblue); xlabel('Time','FontSize',8); ylabel('Stdv [px]','FontSize',8);
            ylim([0, Inf])
            title('Moving Average of Blur Level');
            set(gca,'XTick',0:this.T_train/5:this.T_train,'FontSize',8); grid on

            saveas(gca, strcat(this.savePath, 'blur level'), 'png');
            close(h);

            %% mean VergError
            windowSize = ceil(this.T_train/500);
            ind = 1:this.T_train;
            i = mod(ind,this.Interval);
            ind = find(i); %to exclude trials where vergence is resetted
            ind=ind((this.Interval-1):(this.Interval-1):end);

            vergerr = filter(ones(1,windowSize)/windowSize,1,abs(this.reward_hist(ind,1)));
            vergerr = vergerr(windowSize:end); %discard first values (mean not reliable)
            ind = ind(windowSize:end);
            h = figure; hold on;
            plot(ind,vergerr,'--','LineWidth',1,'Color',blue);
            set(gca,'FontName','LM Roman 8');
            ylim([0 Inf]);
            set(gca,'LineWidth',1);
            set(gca,'XTick',0:this.T_train/2.5:this.T_train,'FontSize',8);
            set(gca,'YTick',0:1:3,'FontSize',8);
            
            title('Moving average of errors during training','FontSize',8,'FontWeight','Normal');

            xlabel('Time','FontSize',8);
            ylabel('Error [px]','FontSize',8);

            saveas(gca, strcat(this.savePath, 'mean VergError'), 'png');
            close(h);

            %% recostruction error
            windowSize = ceil(this.T_train/100);
            recerr = filter(ones(1,windowSize)/windowSize,1,this.recerr_hist);
            h = figure; hold on
            plot(recerr,'k','LineWidth',2);
            xlabel('Time','FontSize',14);
            ylabel('Rec Error [AU]','FontSize',14);
            title('Moving Average of Reconstruction Error');
            set(gca,'XTick',0:this.T_train/10:this.T_train,'FontSize',14); grid on;

            saveas(gca, strcat(this.savePath, 'rec_error'), 'png');
            close(h);

            %% acc and verg error after last step of sequence at the beginning and the end of learning:
            h = figure; hold on;
            if this.T_train >= 1000
                plot(this.vergerr_hist(this.Interval-1+end-1000:this.Interval:end),'b.-','LineWidth',1,'MarkerSize',12);%start
            else
                plot(this.vergerr_hist(this.Interval-1:this.Interval:end),'b.-','LineWidth',1,'MarkerSize',12);
            end
            xlabel('Time','FontSize',14);
            ylabel('VergError [px]','FontSize',14);
            title('VergErr for 250 consecutive iterations','FontSize',10);
            set(gca,'XTick',0:10:250,'FontSize',6); grid on;

            saveas(gca, strcat(this.savePath, 'last 250 last-step-vergErrs'), 'png');
            close(h);

            %% last accerror after last step in sequence
            h = figure; hold on;
            if this.T_train >= 1000
                plot(this.accerr_hist(this.Interval-1+end-1000:this.Interval:end),'b.-','LineWidth',1,'MarkerSize',12);
            else
                plot(this.accerr_hist(this.Interval-1:this.Interval:end),'b.-','LineWidth',1,'MarkerSize',12);
            end
            xlabel('Trial Number','FontSize',14);
            ylabel('AccError [px]','FontSize',14);
            title('AccErr for 250 consecutive iterations','FontSize',10);
            set(gca,'XTick',0:10:250,'FontSize',6); grid on;

            saveas(gca, strcat(this.savePath, 'last 250 last-step-AccErrs'), 'png');
            close(h);

            %% show how acc and verg actions interact
            h = figure; hold on;
            plot(this.accerr_hist(end-100+1:1:end),'.-','LineWidth',1,'MarkerSize',12,'Color',lightblue);
            plot(this.vergerr_hist(end-100+1:1:end),'.-','LineWidth',1,'MarkerSize',12,'Color',blue);
            xlabel('Iteration','FontSize',12);
            ylabel('Error [px]','FontSize',12);
            ylim([-6.5 6.5]);
            title('Vergence and accommodation error in 100 consecutive iterations','FontSize',8);
            legend('acc. error left','verg. error','Location', 'North');grid on;
            set(gca,'XTick',0:10:100,'FontSize',12);
            set(gca,'YGrid','off');

            saveas(gca, strcat(this.savePath, 'verg_and_acc_err'), 'png');
            close(h);

            %% plot accommodation and vergence plane trajectories
            if ~isempty(this.accPlane_hist)
                interv=100;

                h = figure; hold on;
                plot(this.objPlane_hist(this.trainedUntil-interv+1:this.trainedUntil),'--','LineWidth',1,'MarkerSize',12,'Color',gray);
                plot(this.accPlane_hist(this.trainedUntil-interv+1:this.trainedUntil)+this.spectacles_r,'.-','LineWidth',1,'MarkerSize',12,'Color',green);
                plot(this.accPlane_hist(this.trainedUntil-interv+1:this.trainedUntil)+this.spectacles_l,'.-','LineWidth',1,'MarkerSize',12,'Color',lightblue);
                plot(this.vergPlane_hist(this.trainedUntil-interv+1:this.trainedUntil),'.-','LineWidth',1,'MarkerSize',12,'Color',blue);
                xlabel('Iteration','FontSize',12);
                ylabel('Plane Position','FontSize',12);
                % ylim([-this.maxObjPlane-1 this.maxObjPlane+1]);
                title(['Plane positions in ' num2str(interv) ' consecutive iterations - spectacles(l,r):(' num2str(+this.spectacles_l) ',' num2str(+this.spectacles_r) ')'],'FontSize',8);
                legend('Obj. Plane', 'Acc. Plane Right','Acc. Plane Left','Verg. Plane','Location', 'North');grid on;
                set(gca,'XTick',0:10:interv,'FontSize',12);
                set(gca,'YGrid','off');

                saveas(gca, strcat(this.savePath, 'vergErrs-AccErr interaction3'), 'png');
                close(h);
            end

            %% mean errors in each iteration in sequence
            for i=1:this.Interval
                mean_iter_acc_err(i)=mean(abs(this.accerr_hist(this.Interval-1+i:this.Interval:end)));
                mean_iter_verg_err(i)=mean(abs(this.vergerr_hist(this.Interval-1+i:this.Interval:end)));
                sigm_iter_acc_err(i)=std(abs(this.accerr_hist(this.Interval-1+i:this.Interval:end)));
                sigm_iter_verg_err(i)=std(abs(this.vergerr_hist(this.Interval-1+i:this.Interval:end)));
            end

            h = figure; hold on;
            errorbar(1:this.Interval,mean_iter_acc_err,sigm_iter_acc_err);
            errorbar(1:this.Interval,mean_iter_verg_err,sigm_iter_verg_err);
            xlabel('Iteration within Trial','FontSize',14);
            ylabel('Error [px]','FontSize',14);
            legend('Acc. Error','Verg. Error');

            saveas(gca, strcat(this.savePath, 'errors per iteration'), 'png');
            close(h);

            %% error and action path for each step in sequence (also for specific textures)
            err_hi=this.accerr_hist+normrnd(0,0.1,size(this.accerr_hist,1),size(this.accerr_hist,2));
            err_hi=err_hi(1:end-1);

            h = figure; hold on;
            for idx=2:min(1000+2, (size(err_hi,1)+1)/this.Interval-1)
                x_data = 1:this.Interval;
                y_data = err_hi( floor(((size(err_hi,1)+1)/this.Interval-idx) * this.Interval) : floor((((size(err_hi,1)+1)/this.Interval-idx+1)) * this.Interval-1));
                h=plot(x_data,y_data,'Color',blue);
                c = 100*1/min(1000, size(err_hi,1)/this.Interval-2);  % 50% transparent
                if c > 1
                    c = 1;
                elseif c < 0
                    c = 0;
                end

                h.Color(4) = c;
            end

            grid on;
            title('Left Acc. Error Trajectories (scattered)','FontSize',12);
            xlabel('Iteration within Trial','FontSize',12);
            ylabel('Acc. Error Left','FontSize',12);
            xlim([1 this.Interval]);
            ylim([min(err_hi) max(err_hi)]);

            saveas(gca, strcat(this.savePath, 'error and action path for each iteration (also textures)'), 'png');

            %% heatmap of errors during one fixation
            interv = this.trainedUntil/2;
            accerrs = this.objPlane_hist(this.trainedUntil - interv + 1: this.trainedUntil - 1) - (this.accPlane_hist(this.trainedUntil - interv + 1: this.trainedUntil - 1) + this.spectacles_l);
            stepInSeq = repmat(1:this.Interval, 1, interv/10);
            stepInSeq = stepInSeq(2:end);

            blurInd = - min(accerrs);
            maxBlur = max(accerrs) + blurInd;

            pairs = [stepInSeq', accerrs, this.objPlane_hist(this.trainedUntil-interv+1:this.trainedUntil-1), this.accPlane_hist(this.trainedUntil-interv+1:this.trainedUntil-1), this.reward_hist(this.trainedUntil-interv+1:this.trainedUntil-1,3)];
            mapping = zeros(this.Interval, maxBlur + 1);

            for i = 1:length(pairs)
                pair = pairs(i, :);
                mapping(pair(1), pair(2) + blurInd + 1) = mapping(pair(1), pair(2) + blurInd + 1) + 1;
            end
            mapping = mapping';
            
            h = figure;
            imagesc(1:this.Interval, min(accerrs):max(accerrs), mapping);
            colormap('gray');
            title('Left Acc. Error');
            xlabel('Iteration within Trial','FontSize',12);
            ylabel('acc error','FontSize',12);
            set(gca,'XTick',2:2:this.Interval,'FontSize',12);
            set(gca,'YTick',min(accerrs):2:max(accerrs),'FontSize',12);

            saveas(gca, strcat(this.savePath, 'acc errs 2nd half of training'), 'png');
            close(h);

            %% average errors per for each input texture:
            % each dot corresponds to one input image
            ind=this.Interval-1:this.Interval:this.trainedUntil;

            mean_pic_vergerror=accumarray(this.texture_hist(ind),abs(this.vergerr_hist(ind)),[],@(x) mean(x,1));        % accumulates the mean error for each image

            mean_pic_accerror_l=accumarray(this.texture_hist(ind),abs(this.accerr_hist(ind)+this.spectacles_l),[],@(x) mean(x,1));
            mean_pic_accerror_r=accumarray(this.texture_hist(ind),abs(this.accerr_hist(ind)+this.spectacles_r),[],@(x) mean(x,1));

            std_pic_vergerror=accumarray(this.texture_hist(ind),abs(this.vergerr_hist(ind)),[],@(x) std(x,1));

            std_pic_accerror_l=accumarray(this.texture_hist(ind),abs(this.accerr_hist(ind)+this.spectacles_l),[],@(x) std(x,1));
            std_pic_accerror_r=accumarray(this.texture_hist(ind),abs(this.accerr_hist(ind)+this.spectacles_r),[],@(x) std(x,1));

            h = figure; hold on;
            scatter(mean_pic_vergerror,std_pic_vergerror,10,blue);
            scatter(mean_pic_accerror_r,std_pic_accerror_r,10,green);
            scatter(mean_pic_accerror_l,std_pic_accerror_l,10,lightblue);
            title(sprintf('Verg-/Acc-Err for %d input images - LensRight:%dpx', max(this.texture_hist), this.spectacles_r,'FontSize',12));
            ylim([0 2.2]);
            xlim([0 4]);
            set(gca,'XTick',0:1:4,'FontSize',12);
            set(gca,'YTick',0:0.5:2,'FontSize',12);
            xlabel('Mean Error [px]','FontSize',12);
            ylabel('Stdv [px]','FontSize',12);
            legend('Verg. Error','Acc. Error Right','Acc. Error Left','Location', 'NorthEast');
            grid on;

            saveas(gca, strcat(this.savePath, 'error per image'), 'png');
            close(h);

            %% Relation between object distance and left acc error
            objt = this.distances(1:end-1);
            objt = objt';
            acc_err = this.accerr_hist;
            acc_err_seq=acc_err(this.Interval-1:this.Interval:end);  % acc err in last step of sequence

            distIndTransf = - min(objt) + 1;                         % transformation of raw values to indices
            accIndTransf = - min(acc_err_seq) + 1;                   % same for accomodation plane

            %plot heatmap which errors occures with what frequency given the object position
           
            buffer=[objt+distIndTransf acc_err_seq+accIndTransf];
            buffer=accumarray(buffer(:,[1 2]), 1);
            buffer=buffer';
            h = figure;
            colormap('gray');   % set colormap
            imagesc(min(objt):max(objt),min(acc_err_seq)+this.spectacles_l:max(acc_err_seq)+this.spectacles_l,buffer); % draw image and scale colormap to values range
            colorbar;           % show color scale
            title('distribution of acc error in last iteration of sequences for given object distances','FontSize',8);
            xlabel('Object Position','FontSize',12);
            ylabel('Acc. Error Left','FontSize',12);

            saveas(gca, strcat(this.savePath, 'object position vs final AccErr'), 'png');
            close(h);
            
            close all;
        end
    end
end
