function Sim(experimentFlag)

    %###########################################################
    %##################### set  parameters ###########################
    %###########################################################

    rng(1337); % set random seed
    
    resultsFolder = './results';
    savePath = sprintf('%s/%s_%s/', resultsFolder, datestr(now, 'yy-mm-dd_HH-MM'), experimentFlag); % might remain the only variable defined here

    % load/Initialize Model
    loadModel = 0;
    loadModelPath = './Path/to/model/';
    loadWeights = 0;    % whether to use previous values for sparse coder or RL agents
    loadBases = 0;      % load bases from separate file independent of loaded model
    preSample = 1;      % are patches presampled before simulation starts? (decreases runtime for large 'T_train')
    loadSample = 0;     % load presampled patches in advance
    saveSample = 0;     % save presampled patches to file which can be reused
    
    % search in resultsFolder for folder with flags, to load samples from
    if ~isempty(strfind(experimentFlag, 'hc'))        
        sample_Flag = 'hc';   
    elseif ~isempty(strfind(experimentFlag, 'aniso'))        
        sample_Flag = 'aniso';
    elseif ~isempty(strfind(experimentFlag, 'mono'))        
        sample_Flag = 'mono';
    else
       sample_Flag = 'error';
    end
    
    %% create new model & copy source files as backup to results folder
    model = config(savePath); %instantiate model object
    mkdir(savePath);
    copyfile(strcat(mfilename, '.m'), savePath);
    copyfile('config.m', savePath);
    copyfile('Model.m', savePath);
    copyfile('SparseCoding.m', savePath);
    copyfile('ReinforcementLearning.m', savePath);
    copyfile('PatchGenerator.m', savePath);
    copyfile('calculateRightBinocularity.m', savePath);
    copyfile('calculateCorrelation.m', savePath);
    copyfile('laprnd.m', savePath);
    copyfile('Gabor.m', savePath)

    if loadModel  % start with a pretrained model
        modelTrained = load(strcat(loadModelPath, '/model.mat'));
        modelTrained = modelTrained.model;

        % load RLs
        if loadWeights
            model.rlmodel_blur = modelTrained.rlmodel_blur;
            model.rlmodel_disp = modelTrained.rlmodel_disp;
        end

        % load SC Basis
        for s = 1:length(model.scmodel)
            model.scmodel{s} = modelTrained.scmodel{s};
            model.scmodel{s}.Basis_hist(:,:,1) = modelTrained.scmodel{s}.Basis;
        end

        % load parameters for reward normalization
        model.blur_reward_mean = modelTrained.blur_reward_mean;
        model.disp_reward_mean = modelTrained.disp_reward_mean;
        model.blur_reward_variance = modelTrained.blur_reward_variance;
        model.disp_reward_variance = modelTrained.disp_reward_variance;

        clearvars modelTrained  %remove loaded model from workspace to free up RAM
    end

    if loadBases        
        data = load('basis_coarse.mat','basis');
        model.scmodel{1}.Basis = data.basis;
        model.scmodel{1}.Basis_hist(:, :, 1) = data.basis;
        data = load('basis_fine.mat','basis');
        model.scmodel{2}.Basis = data.basis;
        model.scmodel{2}.Basis_hist(:, :, 1) = data.basis;

        copyfile('basis_coarse.mat', savePath);
        copyfile('basis_fine.mat', savePath);
    end

    % prepare parameters for patch generation
    nScales = length(model.scmodel);
    params = {};
    for s = 1:nScales
        params{s} = {model.scmodel{s}.patch_size, model.scmodel{s}.Dsratio, model.blur_scale, model.cataract_l, model.cataract_r};
    end

    %###############################################################
    %################# initialize or load input patches ########################
    %###############################################################

    files = dir([model.texture_directory '*.bmp']);

    if (loadSample && ~exist('sample_Left','var') && ~exist('sample_Right','var'))
        loadingPath = dir(sprintf('%s/*%s-samples',resultsFolder, sample_Flag));
        loaded = 0;
        for directory = 1 : length(loadingPath)
            if loaded
                continue
            end
            try
                display('Load patches_Left.mat');
                load(strcat(resultsFolder, '/' , loadingPath(directory).name, '/patches_Left.mat'));
                load(strcat(resultsFolder, '/' , loadingPath(directory).name, '/rewardsL.mat'));
                display('Load patches_Right.mat');
                load(strcat(resultsFolder, '/' , loadingPath(directory).name, '/patches_Right.mat'));
                load(strcat(resultsFolder, '/' , loadingPath(directory).name, '/rewardsR.mat'));
                display('patches loaded!');
                loaded = 1;
            catch
            end
        end
        if ~loaded
            sprintf('loading failed')
        end
    elseif (~loadSample)
        sample_Left = {};
        sample_Right = {};
        accRewLeft = {};
        accRewRight = {};
        texture_images = {};
        tic;

        for k = 1 : model.n_textures
            img = imread([model.texture_directory files(k).name]);
            if(model.colored_input_images==1)
                img = .2989*img(:,:,1) +.5870*img(:,:,2) +.1140*img(:,:,3); %RGB to grayscale
            end
            img=im2double(img);
            if(size(img,1)~=model.image_size || size(img,2)~=model.image_size)
                img=imcrop(img,[(size(img,1)/2-model.image_size/2) (size(img,2)/2-model.image_size/2) model.image_size-1 model.image_size-1]);
            end
            if(size(img,1)~=model.image_size || size(img,2)~=model.image_size)
                disp('error - incorrect format of input image');
            end
            texture_images{k}=img;
            if preSample
                for blur=-model.maxBlur:1:model.maxBlur
                    for dispa=-model.maxDisp:1:model.maxDisp
                        for s = 1:nScales
                            [patchesLeft, patchesRight, accRewL, accRewR] = PatchGenerator(texture_images{k}, dispa, blur-model.spectacles_l, blur-model.spectacles_r, model, params{s});
                            sample_Left{k, s}{dispa+model.maxDisp+1, blur+model.maxBlur+1} = patchesLeft;
                            sample_Right{k, s}{dispa+model.maxDisp+1, blur+model.maxBlur+1} = patchesRight;
                            accRewLeft{k, s}{dispa+model.maxDisp+1, blur+model.maxBlur+1} = accRewL;
                            accRewRight{k, s}{dispa+model.maxDisp+1, blur+model.maxBlur+1} = accRewR;
                        end
                    end
                end
            end
            display(['texture: ' num2str(k)]); %print rendering progress
        end

        save(strcat(savePath, 'rewardsL.mat'), 'accRewLeft');
        save(strcat(savePath, 'rewardsR.mat'), 'accRewRight');
        toc;
    end

    if (preSample && saveSample)
        display('Save patches_Left.mat');
        save(strcat(savePath, 'patches_Left.mat'),'sample_Left','-v7.3');
        display('Save patches_Right.mat');
        save(strcat(savePath, 'patches_Right.mat'),'sample_Right','-v7.3');
        display('patches saved!');
    end

    
    %###############################################################
    %################# initializing simulation variables ########################
    %###############################################################

    current_texture = randi(model.n_textures);

    %generate object positions uniformly between -3 to 3
    model.distances = randi([-model.maxObjPlane model.maxObjPlane], 1, model.T_train/model.Interval + 1);
    index = 1;
    curr_disp = model.distances(index);
    curr_blur = model.distances(index);
    blurIndexOffset = ceil(model.rlmodel_blur.Action_num / 2);
    dispIndexOffset = ceil(model.rlmodel_disp.Action_num / 2);

    shift = 0;              % initial vergence plane
    lense_blur = 0;         % initial accommodation plane
    monocularity = 0.5;     % initial binocularity index: ranges from -1 (left mon.) to 1 (right mon.)
    contrastLeft = 1;
    contrastRight = 1;
    nFeatVec = 0;           % for initializing the feature vector
    for s = 1:nScales
        nFeatVec = nFeatVec + model.scmodel{s}.Basis_num;
    end
    feature = abs(rand(nFeatVec, 1) - 0.5);     % zero mean init of coeffitients to enable suppression from beginning
    accRewL = {};
    accRewR = {};
    patchesLeft = {};
    patchesRight = {};
    monDetectLeft = 0;
    monDetectRight = 0;
    monDetectLeft = 0;
    monDetectRight = 0;

    runningAvgLC = 0.5;  % initialization for the moving averages
    runningAvgRC = 0.5;
    runningAvgLF = 0.5;
    runningAvgRF = 0.5;

    %% -
    %###########################################################################
    %######################## MAIN SIMULATION ##################################
    %###########################################################################

    tic
    done = 0;                % controle parameter if sim is finished
    t = model.trainedUntil;  % current # of iterations, starts with 0
    while(~done)
        
        % display training process and save model 11 times per simulation
        if(find(t == model.saveAt))
            model.trainedUntil = t;

            disp([num2str(t/model.T_train*100) '% is finished']);
            if model.trainSC
                % save basis
                for s = 1:nScales
                    model.scmodel{s}.saveBasis(find(t==model.saveAt));
                end
            end

            if model.trainRL
                % save weights
                model.rlmodel_blur.saveWeights(find(t==model.saveAt));  %save policy and value net weights
                model.rlmodel_disp.saveWeights(find(t==model.saveAt));  %save policy and value net weights
            end

            % save model
            save(strcat(savePath, 'model.mat'),'model');
        end

        % update variables
        t = t+1;    % iteration counter
        if (t >= model.T_train)
            done = 1;
        end
        
        % change object and its position every 'Interval' iterations
        if ~mod(t,model.Interval)
            index = index + 1;
            current_texture = randi(model.n_textures,1);  % random texture
            curr_disp = model.distances(index);
            curr_blur = model.distances(index);
        end

        if (preSample || loadSample)
            blur = curr_blur - lense_blur + model.maxBlur + 1;   % transform to indices
            dispa = curr_disp - shift + model.maxDisp + 1;

            for s = 1:nScales
                patchesLeft{s} = sample_Left{current_texture, s}{dispa, blur};
                patchesRight{s} = sample_Right{current_texture, s}{dispa, blur};
                accRewL{s} = accRewLeft{current_texture, s}{dispa, blur};
                accRewR{s} = accRewRight{current_texture, s}{dispa, blur};
            end
        else  % online sampling
            d = curr_disp - shift;          % disparity
            b = curr_blur - lense_blur;     % blur
            for s = 1:nScales
                [patchesLeft{s}, patchesRight{s}, accRewL{s}, accRewR{s}] = PatchGenerator(texture_images{current_texture}, d, b-model.spectacles_l, b-model.spectacles_r, model, params{s});
            end
        end

        %#######################################################################
        %#################### Suppression Mechanism ############################
        %#######################################################################

        if model.useSuppression
            for s = 1:nScales
                binocularity = calculateRightBinocularity(model.scmodel{s}.Basis);      % weights to contrast units
                                                                                        % 0: left monocular, 0.5: binocular, 1: right monocular
                if s == 1                                                               
                    activities = feature(1:model.scmodel{s}.Basis_num);                 % feature vector is already normalized since feature vector is normalized
                elseif s == 2
                    activities = feature(model.scmodel{s}.Basis_num+1 : end);
                end

                monDetectLeft = dot(activities, (1 - binocularity));                    % weight binocularities by cortical cell activity
                monDetectRight = dot(activities, binocularity);

                if s == 1
                    runningAvgLC = runningAvgLC + model.eta_supp * (monDetectLeft - runningAvgLC);  % running average of contrast measure
                    runningAvgRC = runningAvgRC + model.eta_supp * (monDetectRight - runningAvgRC);

                    model.monocularity_hist(t, 3) = runningAvgLC;
                    model.monocularity_hist(t, 4) = runningAvgRC;

                    monDetectLeft = max(0, model.slope * (runningAvgLC - model.threshold));  % output of contrast units dependent on slope and threshold
                    monDetectRight = max(0, model.slope * (runningAvgRC- model.threshold));

                    contrastLeft = 1 + model.exci*monDetectLeft - model.suppr*monDetectRight;
                    contrastRight = 1 - model.suppr*monDetectLeft + model.exci*monDetectRight;

                    model.contrast_hist(t, 1) = contrastLeft;  % record contrast values
                    model.contrast_hist(t, 2) = contrastRight;
                    model.monocularity_hist(t, 5) = monDetectLeft;
                    model.monocularity_hist(t, 6) = monDetectRight;

                elseif s == 2
                    runningAvgLF = runningAvgLF + model.eta_supp * (monDetectLeft - runningAvgLF);
                    runningAvgRF = runningAvgRF + model.eta_supp * (monDetectRight - runningAvgRF);

                    model.monocularity_hist(t, 7) = runningAvgLF;
                    model.monocularity_hist(t, 8) = runningAvgRF;

                    monDetectLeft = max(0, model.slope * (runningAvgLF - model.threshold));  % send through linear activation unit
                    monDetectRight = max(0, model.slope * (runningAvgRF - model.threshold));

                    contrastLeft = 1 + model.exci*monDetectLeft - model.suppr*monDetectRight;
                    contrastRight = 1 - model.suppr*monDetectLeft + model.exci*monDetectRight;

                    model.contrast_hist(t, 3) = contrastLeft;
                    model.contrast_hist(t, 4) = contrastRight;
                    model.monocularity_hist(t, 9) = monDetectLeft;
                    model.monocularity_hist(t, 10) = monDetectRight;
                end

                % adjust contrast of input patches patch adaptation
                patchesLeft{s} = patchesLeft{s} .* contrastLeft;
                patchesRight{s} = patchesRight{s} .* contrastRight;
                
                % adjust acc reward according to left/right eye suppression
                accRewL{s} = accRewL{s} * contrastLeft;  
                accRewR{s} = accRewR{s} * contrastRight; 
            end
        end

        if model.trainRL || loadWeights
            [feature, reward, err, coef, monocularity] = model.generateFR(patchesLeft, patchesRight, t, sum([accRewL{:}]), sum([accRewR{:}]));
        end

        %% TRAINING
        
        % if SC AND RL are initialized: update receptive fields
        if model.trainSC && (model.trainRL || loadWeights )
            for s = 1:nScales
                model.scmodel{s}.updateBasis(coef{s},err{s});
            end
            
        %if no RL, encode input and update receptive fields without suppression  
        elseif model.trainSC
            coefs = [];
            for s = 1:nScales
                [err, coef, monocularity] = model.scmodel{s}.stepTrain([patchesLeft; patchesRight]);
                coefs = vertcat(coefs, coef);
            end
            feature = mean(coefs.^2, 2);
            feature = feature./sum(feature); % normalize feature vector
        end

        % get RL actions and then update RL agents
        if model.trainRL
            action_blur = model.rlmodel_blur.stepTrain(feature,reward(1),mod(t-1,model.Interval),0); %index of action
            action_disp = model.rlmodel_disp.stepTrain(feature,reward(2),mod(t-1,model.Interval),0);
        end

        % if no RL model is available, do not move planes and set reward to zero
        if (~model.trainRL && ~loadWeights)
            action_blur = 0 + blurIndexOffset;
            action_disp = 0 + dispIndexOffset;
            reward=[0 0 0];
        end

        %if RL training is off but RL model loaded, choose greedy action, i.e. no exploration
        if (~model.trainRL && loadWeights)
            action_blur = model.rlmodel_blur.Act(feature);  %returns index of action
            action_disp = model.rlmodel_disp.Act(feature);
            reward=[0 0 0];
        end

        % track some variables
        blur_command = model.Act_blur(action_blur);
        shift_command = model.Act_disp(action_disp);
        model.rlmodel_blur.Action_hist(t) = blur_command;
        model.rlmodel_disp.Action_hist(t) = shift_command;
        model.rlmodel_blur.pol_hist(t, :) = model.rlmodel_blur.pol;
        model.rlmodel_disp.pol_hist(t, :) = model.rlmodel_disp.pol;
        model.rlmodel_blur.td_hist(t) = model.rlmodel_blur.td;
        model.rlmodel_disp.td_hist(t) = model.rlmodel_disp.td;
        model.recerr_hist(t, :) = mean(sum([err{:}].^2));  % reconstruction error
        d = curr_disp - shift;          % input disparity
        b = curr_blur - lense_blur;     % defocus blur value
        model.vergerr_hist(t) = d;
        model.accerr_hist(t) = b;
        model.reward_hist(t,:) = [d; b; reward(1); reward(2); reward(3); blur_command; shift_command];
        model.texture_hist(t,1) = current_texture;
        model.monocularity_hist(t, 1:nScales) = monocularity;
        model.objPlane_hist(t) = curr_disp;
        model.vergPlane_hist(t) = shift;
        model.accPlane_hist(t) = lense_blur;

        % apply plane position shifts
        shift = shift + shift_command;
        lense_blur = lense_blur + blur_command;

        if ((abs(shift) > model.maxObjPlane))  % check if eye fixation is out of bound
            shift = model.maxObjPlane * sign(shift);  % set fixation at max/min distance
        end
        if (abs(lense_blur) > model.maxAccDispPlane)  % check if focus reaches near/far-point
            lense_blur = model.maxAccDispPlane * sign(lense_blur);  % set focus at max/min distance
        end
    end
    
    TotT = toc/60;  % total simulation time
    sprintf(' Time [min] = %.2f\non avergage %.1f iters/sec',TotT, model.T_train / (TotT*60))
    model.simulationTime = TotT/60;

    model.trainedUntil = t;  % track when model was last saved

    if model.trainSC
        %save basis
        for s = 1:nScales
            model.scmodel{s}.saveBasis(find(t==model.saveAt));
        end
    end
    if model.trainRL
        %save RL weights
        model.rlmodel_blur.saveWeights(find(t==model.saveAt));  %save policy and value net weights
        model.rlmodel_disp.saveWeights(find(t==model.saveAt));  %save policy and value net weights
    end
    
    %save data
    save(strcat(savePath, 'model.mat'),'model');
    
    % plot results
    model.allPlotSave();
end