function SimOv

%% Init
% goto folder where file located & add to path
[homedir, FileName] = fileparts(mfilename('fullpath'));
% add files and subfolders to matlab path
addpath(genpath(homedir));
cd(homedir);
% close all existing images
close all

%% Parameters
Model = 'sa_z1'; % sensory adapt, no z-scoring
switch Model
    case 'sa_z0' % sensory adapt, no positional sensitivity, no z-scoring
        LhStd  = 10^10; % positional sensitivity: std (not variance): voxel response range / width of
        do_z   = 0;
    case 'sa_z1' % sensory adapt, no pos sens, z-scoring
        LhStd  = 10^10;
        do_z   = 1;
    case 'sa_pos_z1' % sensory adapt, pos sens, no z-scoring
        LhStd  = 10^0;
        do_z   = 1;
    otherwise % do 'sa_z0'
        LhStd  = 10^10;
        do_z   = 0;
end

SNR = 10^1.3;

% Plots
ImgFolder = fullfile(homedir, 'figs');
dircheck(ImgFolder);

%% LOAD experiment data
DataFile = fullfile(homedir, 'Events.mat');
% Ons -- stim onsets in seconds
% Dur -- duration onsets in seconds
[Ons, Dur, Stim] = load_exp_data(DataFile);

% balanced Stim for LDA so num of stims for each class is the same
Stim = kron(ones(1,7), 1:3)';

%% Trim data
% we don't want to simulate all trials but the first N ones so that the
% plots don't become too big
%SubSet = 1:size(Stim,1);
SubSet = 1:20; % subset of trials to use
Stim   = Stim(SubSet);
Ons    = Ons(SubSet);
Dur    = Dur(SubSet);
nStim  = size(Stim,1);

%% SIMULATION
% % population sensitivity to stimuli and prior distribution
% -- voxel's "tuning curves" etc
% can be interpreted as what % of voxels responds most to that stimulus value

% % Simulation Parameters
% stim values -- i.e. possible order positions
StimValues = 1:3;
% number of voxels / dimensions
NumVoxels = 24;
% voxel preferences distribution
PdfType = 'gradient'; % flat | gradient

% create voxel's prefs distribution
switch PdfType
    case 'flat'
        % flat Pdf across order positions
        Pdf = ones(max(StimValues),1);
    case 'gradient'
        % decreasing gradient over order positions --
        % i.e. more voxels respond to items at the beginning of the sequence
        Pdf = [0.5 0.3 0.3];
end
% normalise Pdf to sum = 1
Pdf = Pdf ./ sum(Pdf);

% sample voxel preferences
% i.e. for every voxel we sample the value of the stimulus it is most
% likely to respond to
do_sample = 0;
if do_sample
    % inverse sample voxel preferences according to the Pdf
    PopPref = InvSample(Pdf, StimValues, NumVoxels);
else
    % assign voxel preferences equally ->uniform Pdf
    PopPref = kron(1:3, ones(1, NumVoxels/3));
end

%% generate population responses
% sample actual responses to the presented stimuli
LhDist = 'Normal'; % likelihood distribution
%LhStd  = 10^10;  % std (not variance): voxel response range / width of
% the tuning curve in std; default is 10^0 = 1

% open parallel pool
%par_open
% generate responses
tic
parfor s = 1:size(Stim,1) % for every stimulus presented
    % sample every voxel's response to the presented stimulus
    % based on it's fixed prefernces (PopPref)
    Y(:,s) = arrayfun(@(Mu) pdf(LhDist, Stim(s,1), Mu, LhStd), PopPref);
end
toc
Yy = Y; % save original Y as Yy for later use

%% sensory adaptation for the population
% sensory adaptation vector shows how much each sequence position nodulates
% the signal of ALL the voxels in the ROI
% retain all signal for 1st pos, 0.7 signal for 2nd pos, etc
% plus add some voxel-based SA noise
do_SensAdapt = 1;
if do_SensAdapt
    
    SensAdapt = [1 .7 .4 .2 .1];
    SensAdapt = SensAdapt(Stim);
    %SensAdapt = linspace(0.1,1.1,nStim*NumVoxels);
    %SensAdapt = reshape(SensAdapt',[],NumVoxels)';
    SensAdapt = SensAdapt(ones(NumVoxels,1),:);
    Ys = SensAdapt .* Y;
    
    % add some voxel based noise
    %Ns = random('Normal', 0, mean(Ys(:))/10, NumVoxels, 1);
    %Ns = Ns(:, ones(1,nStim));
    %Ys = Ys + Ns;
    
    Y = Ys;
end

%% population noise
% adding noise to voxel responses depneds on what is the final desired SNR.
% Here SNR (signal-to-noise ratio) is defined as X = S + N, where both
% S and N are random variables. A random variable's % power equals its
% mean-squared value: the signal power thus equals E[S^2].
% If the noise has zero mean (usually), it makes its power equal to its
% variance. Thus, the SNR equals E[S2]/sigma^2_N , where sigma^2_N is the
% variance of noise (N).
% http://www.scholarpedia.org/article/Signal-to-noise_ratio

% calculate the variance of noise given the signal as
% SNR = Signal^2/NoiseSigma
%SNR = 8;
% signal^2 (or, power of the signal)
P = mean(Y(:).^2);
% get required Noise var given the SNR
NoiseSigma = P/SNR;
% for Normal dist matlab takes input as std not var
NoiseStd = sqrt(NoiseSigma);
% same noise Std param for every trial
TrialNoiseStd = ones(nStim,1) .* NoiseStd;

% sample noise for every voxel x trial
NoiseDist = 'Normal';
tic
parfor s = 1:size(Stim,1) % for every trial
    N(:,s) = random(NoiseDist, 0, TrialNoiseStd(s), NumVoxels, 1);
end
toc

% final population response is Ys + N
Yn = Y + N;
% z-score
if do_z
    Yz = zscore(Yn);
else
    Yz = Yn;
end

%% Plots

[~, ReOrdered] = sort(Stim);
yy = {Yy, Yy, Ys, Ys, Yn, Yn, ...
    Yz, Yz, Yz, corr(Yz(:,ReOrdered))};
pf = {@imagesc, @stem, @imagesc, @stem, @imagesc, @stem, ...
    @imagesc, @stem, @lda, @imagesc};
tt = {'Responses', 'Mean', ...
    'Responses+SA', 'Mean+SA', ...
    'Responses+SA+Noise', 'Mean+SA+Noise', ...
    'zscore(Resp+SA+N)', 'Mean zscore(Resp+SA+N)', ...
    'LDA', 'z(SA) similarity -- sorted by trial Pos'};

% subplot properties
numSubPlots = length(yy);
rows = 2;
cols = numSubPlots/rows;
pi = reshape(1:numSubPlots, [],rows);
pi = reshape(pi', 1, []);

% Overview figure

figure('Position', [100 100, 1100, 500])
% subplots
% tight subplots
% h=subtightplot(m,n,p,gap,marg_h,marg_w)
subplot = @(m,n,p) subtightplot(m, n, p, [0.15 0.03], [0.1 0.1], [0.03 0.02]);

for p = 1:numSubPlots
    subplot(rows, cols, pi(p));
    PlotType = func2str(pf{p});
    switch PlotType
        case 'stem'
            for s = 1:3
                stem(find(Stim == s), mean(yy{p}(:,Stim == s)), 'Color', full(sparse(1,s,1,1,3))), hold on
            end
            set(gca, 'YLim', [0 max(mean(Yn))*1.3])
            ylabel('Signal change');
            set(gca, 'XTick', 1:nStim);
            xlim([0 nStim+1]);
            set(gca, 'XTickLabel', Stim);
            xlabel('Stimulus value');
            title(tt{p});
        case 'lda'
            ldaX = Stim';
            ldaY = yy{p};
            LdaOpts.DistClasses = [1 2];
            LdaOpts.Plot = 0;
            Voxels = Lda(ldaX, ldaY, LdaOpts);
            LdaOpts.MarkerSize = 4;
            LdaOpts.FontSize = 10;
            LdaOpts.Voxels = Voxels;
            LdaOpts.Plot = 1;
            LdaOpts.Legend = {'1' '2' '3'};
            Lda(ldaX, ldaY, LdaOpts);
            xlabel(['Voxel ' num2str(Voxels(1))]);
            ylabel(['Voxel ' num2str(Voxels(2))]);
        case 'imagesc'
            pf{p}(yy{p});
            ylabel('Voxel');
            set(gca, 'YTick', 1:2:NumVoxels);
            set(gca, 'XTick', 1:nStim); set(gca, 'XTickLabel', Stim);
            xlabel('Stimulus value');
            title(tt{p});
    end
end % end subplot loop

%export_fig(['fig/' FileName '_sim_' Model], '-png', '-transparent');
%export_fig(fullfile(ImgFolder, [FileName '_sim_' Model]), '-png', '-pdf', '-transparent');



end % end function

function Y = InvSample(Pdf, X, n)
% check that Pdf and X (input values) have the same number of elements
if size(Pdf,2) ~= size(X,2)
    error('Pdf and InputValues not not have the same number of elements');
end
% sample n instances from X according to Pdf
Cdf = cumsum(Pdf); % create cdf
UniRand = rand(n,1); % sample n from Uni(0,1)
% make cdf and unirand matrices compatible
Cdf = Cdf(ones(n,1),:);
UniRand = UniRand(:, ones(1,size(Cdf,2)));
% inverse sample values of X
Y = X( sum(Cdf <= UniRand, 2) + 1);
end

function [Ons, Dur, Stim] = load_exp_data(DataFile)
TR   = 1.75;
load(DataFile, 'E');
L    = E{1,1}.L;
Pre  = strcmp(L(:,3), 'p1') | strcmp(L(:,3), 'p2') | strcmp(L(:,3), 'p3');
Rec  = strcmp(L(:,3), 'r1') | strcmp(L(:,3), 'r2') | strcmp(L(:,3), 'r3');
Stim = cellfun(@(x) x(end), L(:,3));
Ons  = cell2mat(L(:,6));
% filter relevant events
All  = logical(Pre + Rec);
Stim = str2num(Stim(All));
Ons  = Ons(All);
% convert onsets to seconds
Ons  = Ons .* TR;
Dur  = TR * ones(size(Ons,1),1);
end

function dircheck(dirname)

D = 'Dir check: ';
%disp = @(Text) cprintf(rgbconv('3366cc'), Text);
disp = @(Text) fprintf('%s', Text);

if exist(dirname, 'dir') ~= 7;
    mkdir(dirname);
    resp = ['Dir created: ' dirname 10];
else
    resp = ['Dir exists: ' dirname 10];
end

disp(resp);
end
