function varargout = simulator(varargin)
% SIMULATOR M-file for simulator.fig
%      SIMULATOR, by itself, creates a new SIMULATOR or raises the existing
%      singleton*.
%
%      H = SIMULATOR returns the handle to a new SIMULATOR or the handle to
%      the existing singleton*.
%
%      SIMULATOR('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SIMULATOR.M with the given input arguments.
%
%      SIMULATOR('Property','Value',...) creates a new SIMULATOR or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before simulator_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to simulator_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help simulator

% Last Modified by GUIDE v2.5 17-May-2013 11:10:41

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @simulator_OpeningFcn, ...
                   'gui_OutputFcn',  @simulator_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% Otevírací funkce celého simulátoru
function simulator_OpeningFcn(hObject, ~, handles, varargin)
global arrows
% Choose default command line output for simulator
handles.output = hObject;
% Update handles structure
guidata(hObject, handles);

set(handles.DCM_error_uitable,'Visible','off');
% Zobrazení obrázku s mozkem
axes(handles.estimated_model_axes);
axis off
hold off

axes(handles.specified_model_axes);
imshow(imread('images/background_small.png'));

hold on
arrow1_pure = imread('images/arrow1.png');
arrows.arrow1 = imshow(arrow1_pure);
set(arrows.arrow1,'AlphaData',zeros(size(rgb2gray(arrow1_pure))));

hold on
arrow2_pure = imread('images/arrow2.png');
arrows.arrow2 = imshow(arrow2_pure);
set(arrows.arrow2,'AlphaData',zeros(size(rgb2gray(arrow2_pure))));

hold on
arrow3_pure = imread('images/arrow3.png');
arrows.arrow3 = imshow(arrow3_pure);
set(arrows.arrow3,'AlphaData',zeros(size(rgb2gray(arrow3_pure))));

hold on
arrow4_pure = imread('images/arrow4.png');
arrows.arrow4 = imshow(arrow4_pure);
set(arrows.arrow4,'AlphaData',zeros(size(rgb2gray(arrow4_pure))));

hold on
arrow5_pure = imread('images/arrow5.png');
arrows.arrow5 = imshow(arrow5_pure);
set(arrows.arrow5,'AlphaData',zeros(size(rgb2gray(arrow5_pure))));

hold on
arrow6_pure = imread('images/arrow6.png');
arrows.arrow6 = imshow(arrow6_pure);
set(arrows.arrow6,'AlphaData',zeros(size(rgb2gray(arrow6_pure))));

hold off
axis off


% --- Outputs from this function are returned to the command line.
function varargout = simulator_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

set(handles.from1to1_radiobutton,'Enable','off')
set(handles.from2to2_radiobutton,'Enable','off')
set(handles.from3to3_radiobutton,'Enable','off')

set(handles.run_pushbutton,'Enable','off')
set(handles.save_results_as_matfile_pushbutton,'Enable','off')
set(handles.plot_results_pushbutton,'Enable','off')
set(handles.export_results_to_excel_pushbutton,'Enable','off')

% ---------------- Funkce nastavení vstupních parametrů -------------------
function from1to1_radiobutton_Callback(hObject, eventdata, handles)
handles = check_setting(handles);

function from2to1_radiobutton_Callback(hObject, eventdata, handles)
global arrows

handles = check_setting(handles);
arrow1_pure = imread('images/arrow1.png');
from2to1 = get(handles.from2to1_radiobutton,'Value');
if from2to1 == 1
    set(arrows.arrow1,'AlphaData',rgb2gray(arrow1_pure)>80);
else
    set(arrows.arrow1,'AlphaData',zeros(size(rgb2gray(arrow1_pure))));
end


function from3to1_radiobutton_Callback(hObject, eventdata, handles)
global arrows

handles = check_setting(handles);
arrow2_pure = imread('images/arrow2.png');
from3to1 = get(handles.from3to1_radiobutton,'Value');
if from3to1 == 1
    set(arrows.arrow2,'AlphaData',rgb2gray(arrow2_pure)>80);
else
    set(arrows.arrow2,'AlphaData',zeros(size(rgb2gray(arrow2_pure))));
end


function from1to2_radiobutton_Callback(hObject, eventdata, handles)
global arrows

handles = check_setting(handles);
arrow3_pure = imread('images/arrow3.png');
from1to2 = get(handles.from1to2_radiobutton,'Value');
if from1to2 == 1
    set(arrows.arrow3,'AlphaData',rgb2gray(arrow3_pure)>80);
else
    set(arrows.arrow3,'AlphaData',zeros(size(rgb2gray(arrow3_pure))));
end


function from2to2_radiobutton_Callback(hObject, eventdata, handles)
handles = check_setting(handles);

function from3to2_radiobutton_Callback(hObject, eventdata, handles)
global arrows

handles = check_setting(handles);
arrow4_pure = imread('images/arrow4.png');
from3to2 = get(handles.from3to2_radiobutton,'Value');
if from3to2 == 1
    set(arrows.arrow4,'AlphaData',rgb2gray(arrow4_pure)>80);
else
    set(arrows.arrow4,'AlphaData',zeros(size(rgb2gray(arrow4_pure))));
end


function from1to3_radiobutton_Callback(hObject, eventdata, handles)
global arrows

handles = check_setting(handles);
arrow5_pure = imread('images/arrow5.png');
from1to3 = get(handles.from1to3_radiobutton,'Value');
if from1to3 == 1
    set(arrows.arrow5,'AlphaData',rgb2gray(arrow5_pure)>80);
else
    set(arrows.arrow5,'AlphaData',zeros(size(rgb2gray(arrow5_pure))));
end


function from2to3_radiobutton_Callback(hObject, eventdata, handles)
global arrows

handles = check_setting(handles);
arrow6_pure = imread('images/arrow6.png');
from2to3 = get(handles.from2to3_radiobutton,'Value');
if from2to3 == 1
    set(arrows.arrow6,'AlphaData',rgb2gray(arrow6_pure)>80);
else
    set(arrows.arrow6,'AlphaData',zeros(size(rgb2gray(arrow6_pure))));
end


function from3to3_radiobutton_Callback(hObject, eventdata, handles)
handles = check_setting(handles);


% --- Executes on selection change in signal_choice_listbox.
function signal_choice_listbox_Callback(hObject, eventdata, handles)
% hObject    handle to signal_choice_listbox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns signal_choice_listbox contents as cell array
%        contents{get(hObject,'Value')} returns selected item from signal_choice_listbox


% --- Executes during object creation, after setting all properties.
function signal_choice_listbox_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



function signal_mix_ratio_popupmenu_Callback(hObject, eventdata, handles)


function signal_mix_ratio_popupmenu_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end



function number_of_iterations_edit_Callback(hObject, eventdata, handles)
handles = check_setting(handles);

function number_of_iterations_edit_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%------------------- Kontrola nasatevní vstupních parametrů ---------------
function [handles] = check_setting(handles)

from2to1 = get(handles.from2to1_radiobutton,'Value');
from3to1 = get(handles.from3to1_radiobutton,'Value');
from1to2 = get(handles.from1to2_radiobutton,'Value');
from3to2 = get(handles.from3to2_radiobutton,'Value');
from1to3 = get(handles.from1to3_radiobutton,'Value');
from2to3 = get(handles.from2to3_radiobutton,'Value');

model = [0 from2to1 from3to1;...
         from1to2 0 from3to2;...
         from1to3 from2to3 0];
     
no_of_connections = sum(sum(model));

number_of_iterations = get(handles.number_of_iterations_edit, 'String');

if no_of_connections > 0 && isempty(number_of_iterations) == 0 && ...
        str2num(number_of_iterations) > 0
    set(handles.run_pushbutton,'Enable','on')
else 
    set(handles.run_pushbutton,'Enable','off')
end


% ----------------- Spuštění simulace a analýzy ---------------------------
function run_pushbutton_Callback(hObject, eventdata, handles)
global arrows h_out DCM_error_tab ratio fluct correl physiol str clk

set(handles.run_pushbutton,'Enable','off')
set(handles.save_results_as_matfile_pushbutton,'Enable','off')
set(handles.plot_results_pushbutton,'Enable','off')
set(handles.export_results_to_excel_pushbutton,'Enable','off')
waitbar_fig = waitbar(0,'Vyčkejte prosím, probíhá výpočet.','WindowStyle','modal','Color',[1 1 1],'Name','Please wait...');
set(findobj(waitbar_fig,'type','patch'),'edgecolor',[0 1 0],'facecolor',[0 1 0]);


from2to1 = get(handles.from2to1_radiobutton,'Value');
from3to1 = get(handles.from3to1_radiobutton,'Value');
from1to2 = get(handles.from1to2_radiobutton,'Value');
from3to2 = get(handles.from3to2_radiobutton,'Value');
from1to3 = get(handles.from1to3_radiobutton,'Value');
from2to3 = get(handles.from2to3_radiobutton,'Value');

model = [0 from2to1 from3to1;...
         from1to2 0 from3to2;...
         from1to3 from2to3 0];
no_of_connections = sum(sum(model));

fluct = 0;
correl = 0;
physiol = 0;

signals = get(handles.signal_choice_listbox,'Value');
no_of_signals = length(signals);
for i = 1:no_of_signals
    if signals(i) == 1
        fluct = 1;
    elseif signals(i) == 2
        correl = 1;
    elseif signals(i) == 3
        physiol = 1;
    end
end

h_in.no_of_iterations = str2num(get(handles.number_of_iterations_edit,'String'));
% Nastavení počtu iterací. Podle počtu iterací bude počet výsledných buněk v mat souboru
h_in.SNR=3;%  SNR levels for simulation (- výkonové SNR, podle vzorečku, co je na wikipedii)
h_in.ROI=3;    % Number of ROIs
h_in.TR =2;    % Repeat time [s]
h_in.dt =5;    % Time step fine sampling [ms]
h_in.N  =200;                % Number of time points (BOLD signal after downsampling by TR)
h_in.M(1).dt=h_in.dt/1000;      % Time step fine sampling [s]
h_in.T  =h_in.TR*h_in.N/h_in.M(1).dt; % Number of elements in fine sampling
h_in.M(1).Ud_dt = 0.1; % Time step for stimulation vector for DCM analysis
h_in.no_of_events = 15;% Number of stimuls for stimulation vector(default)
h_in.equal_length = 5; % [s] If selected, all stimulations in vector will last this time
%---------------------------------------------------------------------

h_in.M(1).C = zeros(h_in.ROI,h_in.ROI);
h_in.M(1).C(1,1)=1;
%----------------------------------------------------------------------
% Stimulation vector
h_in.event_duration = 500;    % [ms]

% Positions of every event - scans or secs according to flag in CreateU
h_in.Pos = [1 5.5 17 24.3 30 60 80 100 120 140 150 160 170 180 190];
% Durations of every event - scans or secs according to flag in CreateU
h_in.Dur = [2 10 3 3 20 15 15 15 15 5 5 5 5 5 5 ];
% Amplitude of stimulation vector
h_in.Ampl = [1];
% Signal mix ratio
cells = get(handles.signal_mix_ratio_popupmenu,'String');
ratio = str2num(cells{get(handles.signal_mix_ratio_popupmenu,'Value')});

DCM_error_tab = zeros(h_in.no_of_iterations,1);

% Inicializace promenných pro parfor
parameter_index = cell(1,h_in.no_of_iterations);
h_out = cell(1,h_in.no_of_iterations);
for i = 1:h_in.no_of_iterations
    h_out{i} = h_in;
end

% Changing Parameters - models, connections
for iteration = 1:h_in.no_of_iterations
    waitbar(iteration /(h_in.no_of_iterations+2))
    h_out{iteration}.M(1).d = [];
    h_out{iteration}.M(1).A = [];
    DCM_A  = [];
    DCM_Ap = [];
    DCM_sA = [];
    DCM_error = [];
    
h_out{iteration}.M(1).d = ones(h_out{iteration}.ROI,h_out{iteration}.ROI)*50;% Matice zpoždění, tvar má jako A, ale hodnoty jsou zpoždění signálu mezi ROI v ms
h_out{iteration}.M(1).d = h_out{iteration}.M(1).d.*logical(model);

% ----------------- DCM data generation -----------------------------------
diagonal = 0.5;
if no_of_connections > 3
    diagonal = 1;
end
% h.M(1).A=((rand(h.ROI)*0.15)+0.3);  % Random values of simulated connectivity matrix: 0.3-0.45 (0.45=0.35+0.1)
% h.M(1).A=h.M(1).A.*model;           % Select values on valid positions

h_out{iteration}.M(1).A = ((rand(h_out{iteration}.ROI)*0.15)+0.3);  % parametry nastavené do intervalu 0.3-0.45 (0.45=0.35+0.1)
h_out{iteration}.M(1).A = h_out{iteration}.M(1).A.*model; 
h_out{iteration}.M(1).A = h_out{iteration}.M(1).A-diagonal*eye(h_out{iteration}.ROI);  % Diagonal values -0.5 (or -1)(DCM)

% Creating stimulation vector
h_out{iteration}.xU = CreateU(h_out{iteration}.Pos, h_out{iteration}.Dur, h_out{iteration}.Ampl, h_out{iteration}.N, h_out{iteration}.TR, h_out{iteration}.M(1).dt,'scans', h_out{iteration}.event_duration);
% Specifying the input
U = repmat(h_out{iteration}.xU',h_out{iteration}.ROI,1);
h_out{iteration}.M(1).U = U;
% Creating stimulation vector for DCM estimation
h_out{iteration}.xUd = CreateU(h_out{iteration}.Pos, h_out{iteration}.Dur, h_out{iteration}.Ampl, h_out{iteration}.N, h_out{iteration}.TR, h_out{iteration}.M(1).Ud_dt,'scans', h_out{iteration}.event_duration);
h_out{iteration}.M(1).Ud = repmat(h_out{iteration}.xUd,1,h_out{iteration}.ROI);

h_out{iteration} = create_DCM(h_out{iteration});                    % Generating BOLD data by DCM

h_out{iteration} = signal_mix(h_out{iteration},fluct,correl,physiol,no_of_signals,ratio); 
waitbar(iteration /(h_in.no_of_iterations+1))
h_out{iteration} = analyze_DCM(h_out{iteration});   % DCM analysis
waitbar(iteration /h_in.no_of_iterations)


%% Method based on thresholding to count number of incorrectly estimated
%% connections

% --- Structural similarity estimation ---
model_A = h_out{iteration}.M(1).A;                           % Model matrix
sA = logical(model_A.*(~eye(length(model_A)))); % Diagonal values are not included

DCM_A  = h_out{iteration}.DCM.Ep.A;
DCM_Ap = h_out{iteration}.DCM.Pp.A;
p=0.9;             % Threshold for DCM
DCM_A(DCM_Ap<p) = 0;  % Thresholding
DCM_sA = logical(DCM_A.*(~eye(length(DCM_A)))); % Diagonal values are not included

DCM_error = abs(sA-DCM_sA);
DCM_error = sum(DCM_error(:)); % Number of structural differences between DCM and model

DCM_error_tab(iteration,1) = DCM_error;

end
close(waitbar_fig)
disp([''])
disp(['*******************************************************************************'])
disp(['Simulation complete'])
disp([''])

str = date;
clk = clock;
no_of_rows = h_out{1}.no_of_iterations;
width = 107;
if no_of_rows > 9
    width = 110;
end

set(handles.DCM_error_uitable,'Visible','on');
set(handles.DCM_error_uitable,'Position',[38 560-no_of_rows*18 width 22+18*no_of_rows]);
set(handles.DCM_error_uitable,'Data',DCM_error_tab(1:no_of_rows),'ColumnName','DCM error'); 
set(handles.DCM_error_uitable,'CellSelectionCallback',{@click_on_cell, handles})

set(handles.save_results_as_matfile_pushbutton,'Enable','on')
set(handles.plot_results_pushbutton,'Enable','on')
set(handles.export_results_to_excel_pushbutton,'Enable','on')
set(handles.run_pushbutton,'Enable','on')

function h=create_DCM(h)
% Function is used to create simulated BOLD data according DCM method

%---------------------------------------------------------------------
% Inputs 
h.M(1).n    = h.ROI;       % Number of ROIs
h.M(1).N    = h.N;         % Number of time points (BOLD signal after downsampling with TR)
h.M(1).TR   = h.TR;        % Repeat Time [s]
h.M(1).T    = (h.M(1).N*h.M(1).TR)/h.M(1).dt;   % Number of elements in fine sampling
h.M(1).SNR    = h.SNR;     % Signal to noise ratio of the simulated BOLD signal
h.M(1).G      = [ones(1,5)*0.001]';   % Varince of the additive Wiener process noise to the hemodynamic states
                        % Don't go higher than 0.02 
%---------------------------------------------------------------------
% excitatory connections:
h.M(1).B = zeros(h.ROI,h.ROI,0); % empty so far 
%---------------------------------------------------------------------
% nonlinear term:
h.M(1).D = zeros(h.ROI,h.ROI,0); % empty so far
%---------------------------------------------------------------------
% coding delay matrix
D = h.M(1).d/(h.M(1).dt*1000) + eye(size(h.M(1).d)); % conversion from ms to time steps; 
Db = logical(D);
Ds = [1:h.ROI];
Ds = Ds(ones(h.ROI,1),:);
Ds = Ds.*Db;
D  = D(:); D(D==0) = [];
Ds = Ds(:); Ds(Ds==0) = [];
for i = 1:length(D)
  Dc{i} = zeros(1,D(i));  % delay sequences
end
h.M(1).Db = Db;
h.M(1).Ds = Ds;
h.M(1).Dc = Dc;
%---------------------------------------------
%
% hemodynamic model parameters:
%    (a)   (b)  (c)  (d)  (e)
% --------------------------------------------
parameters = [0.64 0.32 1.00 0.32 0.32]; % Pro každé ROI je jeden řádek
h.M(1).H = repmat(parameters,h.ROI,1);

% (a) Neuronal efficacy
% (b) Signal decay
% (c) Flow-dependent elimination
% (d) Hemodynamic transit time
% (e) Grubb's exponent
%--------------------------------------------------------------------------
h.SIM = DCM_generateDE(h.M);

function h=analyze_DCM(h)

DCM.a=ones(h.ROI,h.ROI);   % Possible endogenous connections
DCM.b=zeros(h.ROI,h.ROI);  % Possible bilinear modulations - zeros
DCM.c=zeros(h.ROI,h.ROI);  % Possible exogenous connections
DCM.c(1,1)=1; 
DCM.U.u=h.M(1).Ud;         % Exogenous inputs - stimulation vector
DCM.U.dt=h.M(1).Ud_dt;     % Time step of stimulation vector

DCM.Y.y=h.SIM.y;           % Downsampled BOLD with noise
DCM.Y.Q = spm_Ce(ones(1,h.ROI)*h.N);
DCM.Y.X0=[];               % Confounds - empty
DCM.Y.dt=h.TR;             % TR sampling of BOLD
DCM.n=h.ROI;               % Number of regions
DCM.v=h.N;                 % Number of scans

DCM.options.nograph=1;     % Turning off SPM interactive screen

[h.DCM] = spm_dcm_estimate(DCM);

function [U] = CreateU (Pos, Dur, Ampl, nscan, TR, dt, unit_flag, event_duration);
% Create fine-sampled activation/stimulation time-series
% function [U] = CreateU (Pos, Dur, Ampl, nscan, TR, dt, unit_flag, event_duration);
%
%   Pos    =  vector of positions (begins of each trial/condition/event in scans or secs)
%   Dur    =  vector of durations  (durations of each trial in scans or secs)
%   Ampl   =  vector of amplitudes (amplitudes of each trial)
%   dt     =  fine resolution step (in seconds);
%   nscan  =  number of scans
%   TR     =  scan repetition time in seconds
%   unit_flag = definition of units for input parameters (Pos and Dur)
%                'secs' or 'seconds' for seconds
%                'scans' for scans
%   event_duration = duration for events defined with zero Dur (in miliseconds)
%                     if zero then duration = TR
%
% *************************************************************************
%  Written by Michal Mikl
%  last update 2011-03-08
%  version 1.0
% *************************************************************************

% main routine
%**************************************************************************

FSfactor = TR/dt;    % fine sampling factor
nf = FSfactor*nscan; % number of time points using fine sampling

nPos = length(Pos);
nAmpl = length(Ampl);
nDur = length(Dur);

if event_duration == 0
    event_duration = tr*1000;
end
evdur = round(event_duration/dt/1000); %event duration for zero-legnth events in fine samples

% ensure the same length of Pos, Dur and Ampl vectors
if nAmpl > nPos
    Ampl=Ampl(1:nPos);
end
if nDur > nPos
    Dur=Dur(1:nPos);
end
if nAmpl < nPos
    Ampl = Ampl*ones(1,nPos);
end;
if nDur < nPos
    Dur = Dur*ones(1,nPos);
end;

% convert seconds to scans
if (strcmp(unit_flag,'seconds')) || (strcmp(unit_flag,'secs'))
    Pos = Pos / TR;
    Dur = Dur / TR;
end


fsPos = ceil(Pos * FSfactor);        % fine sampled positions
fsDur = round(Dur * FSfactor);         % fine sampled durations
fsAmpl = Ampl;                         % fine sampled amplitudes

for i=1:nPos
    if fsDur(i)==0
        fsDur(i)=evdur;
        fsAmpl(i)=1;
    end;
end;

tmpseries=zeros(nf,1);        % raw (zero-filled) fine sampled time series

% fill in timeseries positions with amplitudes = build experimental timecourse
for i=0:nf-1
    k=find(fsPos==i);
    if (k>0)
        for l=1:length(k)
        for j=1:fsDur(k(l))
            tmpseries(i+j,1)=tmpseries(i+j,1)+fsAmpl(k(l));
        end
        end
    end
end

U = tmpseries;

function [y] = spm_gx_fmri2(x,u,P,M)
% Simulated BOLD response to input
% FORMAT [y] = spm_gx_fmri(x,u,P,M)
% y          - BOLD response (%)
% x          - state vector     (see spm_fx_dcm)
% P          - Parameter vector (see spm_fx_dcm)
% M          - model specification structure (see spm_nlsi)
%__________________________________________________________________________
%
% This function implements the BOLD signal model described in: 
%
% Stephan KE, Weiskopf N, Drysdale PM, Robinson PA, Friston KJ (2007)
% Comparing hemodynamic models with DCM. NeuroImage 38: 387-401.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% modified version MH 2011
% Karl Friston & Klaas Enno Stephan
% $Id: spm_gx_fmri.m 3888 2010-05-15 18:49:56Z karl $
 
 
% Biophysical constants for 1.5T
%==========================================================================
 
% time to echo (TE) (default 0.04 sec)
%--------------------------------------------------------------------------
try, TE = M.TE; catch, TE = 0.04; end
 
% resting venous volume (%)
%--------------------------------------------------------------------------
V0  = 4; % check it!!!

% estimated region-specific ratios of intra- to extra-vascular signal 
%--------------------------------------------------------------------------
ep  = 1; % ? check it!!!
 
% slope r0 of intravascular relaxation rate R_iv as a function of oxygen 
% saturation S:  R_iv = r0*[(1 - S)-(1 - S0)] (Hz)
%--------------------------------------------------------------------------
r0  = 25;
 
% frequency offset at the outer surface of magnetized vessels (Hz)
%--------------------------------------------------------------------------
nu0 = 40.3; 
 
% resting oxygen extraction fraction
%--------------------------------------------------------------------------
E0  = 0.4;
 
%-Coefficients in BOLD signal model
%==========================================================================
k1  = 4.3*nu0*E0*TE;
k2  = ep*r0*E0*TE;
k3  = 1 - ep;
 
%-Output equation of BOLD signal model
%==========================================================================
x   = x';
v   = (x(:,4));
q   = (x(:,5));
y   = V0*(k1.*(1 - q) + k2.*(1 - q./v) + k3.*(1 - v));


function [y P] = spm_fx_fmriD(x,u,P)
% state equation for a dynamic [bilinear/nonlinear/Balloon] model of fMRI
% responses
% FORMAT [y] = spm_fx_fmri(x,u,P,M)

% x      - state vector
%   x(:,1) - excitatory neuronal activity             ue
%   x(:,2) - vascular signal                          s
%   x(:,3) - rCBF                                  ln(f)
%   x(:,4) - venous volume                         ln(v)
%   x(:,5) - deoyxHb                               ln(q)
%  [x(:,6) - inhibitory neuronal activity             ui]
%
% y      - dx/dt
%
%___________________________________________________________________________
%
% References for hemodynamic & neuronal state equations:
% 1. Buxton RB, Wong EC & Frank LR. Dynamics of blood flow and oxygenation
%    changes during brain activation: The Balloon model. MRM 39:855-864,
%    1998.
% 2. Friston KJ, Mechelli A, Turner R, Price CJ. Nonlinear responses in
%    fMRI: the Balloon model, Volterra kernels, and other hemodynamics.
%    Neuroimage 12:466-477, 2000.
% 3. Stephan KE, Kasper L, Harrison LM, Daunizeau J, den Ouden HE,
%    Breakspear M, Friston KJ. Nonlinear dynamic causal models for fMRI.
%    Neuroimage 42:649-662, 2008.
% 4. Marreiros AC, Kiebel SJ, Friston KJ. Dynamic causal modelling for
%    fMRI: a two-state model.
%    Neuroimage. 2008 Jan 1;39(1):269-78.
%__________________________________________________________________________

% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Modified version MH 2011
% Karl Friston & Klaas Enno Stephan
% $Id: spm_fx_fmri.m 3888 2010-05-15 18:49:56Z karl $

xlag = zeros(size(x,2));
Dbi  = find(P.Db==1);
for i=1:length(Dbi)
    DD  =[];
    DD = [x(1,P.Ds(i)),P.Dc{i}];
    xlag(Dbi(i)) = DD(end-1);
    P.Dc{i} = DD(1:end-1);
end
% Neuronal motion
%==========================================================================
B   = full(P.B);                       % bi-linear parameters
C   = P.C/2;                           % exogenous parameters
%D   = full(P.D);                       % nonlinear parameters
x     = x';
% excitatory connections
%--------------------------------------------------------------------------
for i = 1:size(P.B,3)
    P.A = P.A + u(i)*B(:,:,i);
end

% and nonlinear (state) terms
%--------------------------------------------------------------------------
% for i = 1:size(P.D,3)
%     P.A = P.A + x(i,1)*D(:,:,i);
% end

% implement differential state equation y = dx/dt (neuronal)
%--------------------------------------------------------------------------
%try, global sigma; if isempty(sigma), sigma = 1; end; end;
y    = x;
if size(x,2) == 5
    
    % one neuronal state per region
    %----------------------------------------------------------------------
    y(:,1) = sum(P.A.*xlag,2) + C*u(:); % 

else

%     % extrinsic (two neuronal states)
%     %----------------------------------------------------------------------
%     A      = exp(P.A)/8;             % enforce positivity
%     IE     = diag(diag(A));          % inhibitory to excitatory
%     EE     = A - IE;                 % excitatory to excitatory
%     EI     = 1;                      % excitatory to inhibitory
%     SE     = 1;                      % self-inhibition (excitatory)
%     SI     = 2;                      % self-inhibition (inhibitory)
% 
%     % motion - excitatory and inhibitory: y = dx/dt
%     %----------------------------------------------------------------------
%     y(:,1) = EE*x(:,1) - SE*x(:,1) - IE*x(:,6) + P.C*u(:);
%     y(:,6) = EI*x(:,1) - SI*x(:,6);

end

% Hemodynamic motion
%==========================================================================

% hemodynamic parameters
%--------------------------------------------------------------------------
%   H(1) - signal decay                                   d(ds/dt)/ds)
%   H(2) - autoregulation                                 d(ds/dt)/df)
%   H(3) - transit time                                   (t0)
%   H(4) - exponent for Fout(v)                           (alpha)
%   H(5) - resting oxygen extraction                      (E0)
%   H(6) - ratio of intra- to extra-vascular components   (epsilon)
%          of the gradient echo signal
%--------------------------------------------------------------------------
%H        = [0.65 0.41 2.00 0.32 0.34];
try 
   H     = P.H;    
catch    
   H     = [0.64 0.32 2.00 0.32 0.32];
end

% signal decay
%--------------------------------------------------------------------------
sd       = H(:,1);

% transit time
%--------------------------------------------------------------------------
tt       = H(:,3);

% Fout = f(v) - outflow
%--------------------------------------------------------------------------
fv       = x(:,4).^(1./H(:,4));

% e = f(f) - oxygen extraction
%--------------------------------------------------------------------------
ff       = (1 - (1 - H(:,5)).^(1./x(:,3)))./H(:,5);


% implement differential state equation y = dx/dt (hemodynamic)
%--------------------------------------------------------------------------
y(:,2)   = x(:,1) - sd.*x(:,2) - H(:,2).*(x(:,3) - 1);
y(:,3)   = x(:,2);
y(:,4)   = (x(:,3) - fv)./(tt);
y(:,5)   = (ff.*x(:,3) - fv.*x(:,5)./x(:,4))./(tt);
y        = y';
y        = y(:);

return

function [SIM] = DCM_generateDE(M)
% Generate BOLD signal according DCM model
% *************************************************************************
%  Written by Michal Mikl and Martin Havlíček
% *************************************************************************

% initialization of hemodynamic states
x = repmat([0 0 1 1 1]',1,M(1).n);              % x(1,:) = neuronal signal
                                                % x(2,:) = vascular signal
                                                % x(3,:) = rCBF 
                                                % x(4,:) = venous volume 
                                                % x(5,:) = deoyxHb   
% Inputs
DCM(1).P.H  = M(1).H;
DCM(1).P.A  = M(1).A;
DCM(1).P.C  = M(1).C;%/sqrt(0.001);
DCM(1).P.B  = M(1).B; 
DCM(1).P.D  = M(1).D;
DCM(1).U    = M(1).U;
DCM(1).P.Ds = M(1).Ds;
DCM(1).P.Db = M(1).Db;
DCM(1).P.Dc = M(1).Dc;

EXPm   = kron(speye(M(1).n),sparse([ones(5),2*ones(5,1);zeros(1,5+1)]));
xt     = repmat([zeros(1,5) 1]',M(1).n,1);
J      = repmat([1 0 0 0 0;
                 1 1 1 0 0;
                 0 1 0 0 0;
                 0 0 1 1 0;
                 0 0 1 1 1;],1,M(1).n);
             
ww    =  zeros(5*M(1).n,1);

X        = zeros(M(1).T,M(1).n*5);
Y        = zeros(M(1).T,M(1).n);

for t = 1:M(1).T,
       if t==400
         disp('k')  
       end   
    
        [f DCM(1).P] = spm_fx_fmriD(x,DCM(1).U(:,t),DCM(1).P); % hemodynamic state equation
        % integration by euler scheme:
        %------------------------------------------------------------------
        dx       = f * M(1).dt;
        x        = x + reshape(dx + ww,5,M(1).n);
        y        = spm_gx_fmri2(x,[],DCM(1).P);    % observation equation:
        X(t,:)   = x(:)';            
        Y(t,:)   = y(:)';
end

% downsample the BOLD signal
Yd      =   Y(1:(1/(M(1).T*M(1).dt))*M(1).T*M(1).TR:end,:);


% save simulated data:
%-----------------------------------------------------------
SIM(1).x  = X;   % simulated hemodynamic states 
SIM(1).y  = Yd; % simulated hemodynamic response - downsampled and with additive noise 
SIM(1).y0 = Y;   % simulated hemodynamic response - high resolution - no noise

function h = signal_mix (h,fluct,correl,physiol,no_of_signals,ratio)

Yd = h.SIM.y;
nsamples = length(Yd);

% Low frequency fluctuations
%--------------------------------------------------------------------------
f = 0.8+0.4.*rand(1); % 0,8 - 1,2 Někam ukládat?
t = 1:200;
fluct_sin = 10*sin(2*pi*f*t);
% fluct_cos = cos(2*pi*f*t);

% fluct_sin = fluct_sin';
% Ydn(:,1) = Ydn(:,1) + 0.5*fluct_sin; % add to the first ROI

fluct_sin = repmat(fluct_sin',1,3);

% Highly correlated signal
%--------------------------------------------------------------------------
% Create U with 200 samples
pos = h.Pos;
npos = length(pos);
dur = h.Dur;
ndur = length(dur);
U_R = zeros(nsamples,1);
for sample = 1:nsamples
    for j = 1:npos
        if sample == round(pos(j))
            U_R(sample:sample+dur(j)-1) = U_R(sample)+1;
        end
    end
end

% Create rectangular signal
% pos = pos(1:2:end);
% for i = 1:2:ndur-1
%     dur(i) = dur(i) + dur(i+1);
% end
% dur = dur(1:2:end);
% rect = zeros(nsamples,1);
% npos = length(pos);
% for sample = 1:nsamples
%     for j = 1:npos
%         if sample == round(pos(j))
%             rect(sample:sample+dur(j)-1) = rect(sample)+1;
%         end
%     end
% end
rect = U_R;
rect(pos(5):pos(5)+dur(5)) = 0; % delete stimules
rect(pos(11):pos(11)+dur(11)) = 0;
hrf = spm_hrf(h.TR);
con = conv(rect,hrf,'same'); % convoles rectangular signal and hrf
% con = con(4:nsamples+3); for correlation
% R = corrcoef(con,U_R) % correlation coefficients
% h.R = R; % save correlation coeffs
con = 10*con; % amplitude as Yd
con = repmat(con,1,3);

% Physiological noise
%--------------------------------------------------------------------------
% 10x pro 0.1 a 0.9, 0.2 a 0.8..., poslat martinovi graf
load('TR166.mat')
pnoise = TR166(9).data_detr;
pnoise = pnoise(1:nsamples);
pnoise = repmat(pnoise,1,3);

% Add to BOLD signal
%--------------------------------------------------------------------------
N = h.M.N;
n = h.M.n;
SNR = h.M.SNR;
z = randn(N,n);
% standard deviation of additive noise, for all areas
r = diag(std(Yd)/SNR);
gauss = z*r;
Ydn = (1-ratio)*Yd + ratio*(fluct*fluct_sin + correl*con + physiol*pnoise)/no_of_signals + gauss; 
% add to all ROIs

% Save simulated data:
%--------------------------------------------------------------------------
h.SIM.y = Ydn; % simulated hemodynamic response mixed with other signals

function DCM=spm_dcm_post_hoc_MG_02(P)
% Modified spm_dcm_post_hoc(P), 23.3.2013, MG
% Function for posthoc evidence approximation.
% Neccesary inputs:
%
% DCM.Ep - posteriorní odhady - střední hodnoty parametrů (matice A,B,C,D, transit, decay a epsilon)
% DCM.Cp - posteriorní kovariance parametrů - matice o rozměru (parametrů x parametrů)
% 
% DCM.M.pC - DCM.M jsou apriorní parametry modelu; DCM.M.pC je kovarianční matice pro apriorní odhad (kovariance mezi parametry)- matice o rozměru (parametrů x parametrů)
% DCM.M.pE - DCM.M.pE jsou apriorní parametry o modelu (matice A,B,C,D, transit, decay a epsilon)
%--------------------------
% DCM.Ce - Error variance for each region
% DCM.y  - Predicted responses (predicted data)
% DCM.R  - Residuals
% DCM.T  - Threshold for inference based on posterior p.d.f (Threshold for Posterior inference)
% Output for direct connectivity matrix:
% DCM.Ep.A
%
%
% Post hoc optimisation of DCMs (under Laplace approximation)
% FORMAT spm_dcm_post_hoc(P)
%
% P         - character/cell array of DCM filenames
%           - or cell array of DCM structures
%
%--------------------------------------------------------------------------
% This routine searches over all possible reduced models of a full model 
% (DCM) and uses post hoc model selection to select the best. Reduced 
% models mean all permutations of free parameters (parameters with a non-
% zero prior covariance), where models are defined in terms of their prior 
% covariance. The full model should be inverted prior to post hoc 
% optimization. If there are more than 16 free-parameters, this routine 
% will implement a greedy search: This entails searching over all 
% permutations of the 8 parameters whose removal (shrinking the prior 
% variance to zero) produces the smallest reduction (greatest increase) 
% in model evidence. This procedure is repeated until all 8 parameters 
% are retained in the best model or there are no more parameters to 
% consider. When several DCMs are optimized together (as in group studies), 
% they are checked to ensure the same free parameters have been specified 
% and the log-evidences are pooled in a fixed effects fashion.
% 
% This application of post hoc optimization assumes the DCMs that are 
% optimized are the same model of different data. Normally, this would be 
% a full model, in the sense of having the maximum number of free 
% parameters, such that the set of reduced models is as large as possible. 
% In contrast spm_dcm_search operates on different DCMs of the same data 
% to identify the best model
% 
% The outputs of this routine are graphics reporting the model reduction 
% (optimisation) and a DCM_opt_??? for every input DCM that contains the 
% reduced conditional parameters estimates (for simplicity, the original 
% kernels and predicted states are retained). The structural and function 
% (spectral embedding) graphs are based on Bayesian parameter averages 
% over multiple DCMs.
%__________________________________________________________________________
% Copyright (C) 2008-2011 Wellcome Trust Centre for Neuroimaging
 
% Karl Friston
% $Id: spm_dcm_post_hoc.m 4261 2011-03-24 16:39:42Z karl $
 
% get filenames

%--------------------------------------------------------------------------
% Dopočítání proměnných
P.n=size(P.y,2);
P.v=size(P.y,1);
%--------------------------------------------------------------------------

if ischar(P),   P = cellstr(P); end
if isstruct(P), P = {P}; end
N = numel(P);
 
%-Check models are compatible in terms of their prior variances
%==========================================================================
for j = 1:N
    
    % get prior covariances
    %----------------------------------------------------------------------
    try, load(P{j}); catch, DCM = P{j}; end
    pC   = diag(DCM.M.pC);
    
    
    % and compare it with the first model
    %----------------------------------------------------------------------
    if j == 1
        C = pC;
    else
        if any(~pC - ~C)
            fprintf('Please check model %i for compatibility',j)
            return
        end
    end
end
 
% find free coupling parameters
%--------------------------------------------------------------------------
k     = spm_fieldindices(DCM.Ep,'A','B','D');
k     = k(~~C(k));
n     = length(C);


% If there are too many find those with the least evidence
%--------------------------------------------------------------------------
if length(k) > 16
    
    %-Loop through DCMs and free parameters and get log-evidences
    %----------------------------------------------------------------------
    for j = 1:N
        
        try, load(P{j}); catch, DCM = P{j}; end
        
        % Get priors and posteriors
        % -----------------------------------------------------------------
        qE    = DCM.Ep;
        qC    = DCM.Cp;
        pE    = DCM.M.pE;
        pC    = DCM.M.pC;
        
        % model search over new prior without the i-th parameter
        % -----------------------------------------------------------------
        for i = 1:length(k)
            R      = speye(n,n) - sparse(k(i),k(i),1,n,n);
            Z(i,j) = spm_log_evidence(qE,qC,pE,pC,pE,R*pC*R);
        end
    end
    
    % find parameters with the least evidence
    %----------------------------------------------------------------------
    Z      = sum(Z,2);
    [Z i]  = sort(-Z);
    k      = k(i(1:8));
    
    % flag a greedy search
    %----------------------------------------------------------------------
    repeat = 1;
    
else
    repeat = 0;
end

% Create model space in terms of free parameter indices
%--------------------------------------------------------------------------
K     = spm_perm_mtx(length(k));
 
%-Loop through DCMs and models and get log-evidences
%==========================================================================
for j = 1:N
    
    try, load(P{j}); catch, DCM = P{j}; end
    fprintf('\nsearching (%i): 00%%',j)
    
    % Get priors and posteriors
    % ---------------------------------------------------------------------
    qE    = DCM.Ep;
    qC    = DCM.Cp;
    pE    = DCM.M.pE;
    pC    = DCM.M.pC;
    
    % model search over new prior (covariance)
    % ---------------------------------------------------------------------
    for i = 1:length(K)
        R      = speye(n,n) - sparse(k,k,K(i,:),n,n);
        G(i,j) = spm_log_evidence(qE,qC,pE,pC,pE,R*pC*R);
    
        fprintf('\b\b\b\b')
        fprintf('%-3.0f%%',i*100/length(K))
    end 
end
fprintf('\n')
 
% Pooled model evidence
% =========================================================================
S     = sum(G,2);
S     = S - S(end);
p     = exp(S - max(S));
p     = p/sum(p);
 
% Get selected model (prior covariance)
%--------------------------------------------------------------------------
[q i] = max(p);
R     = speye(n,n) - sparse(k,k,K(i,:),n,n);
rC    = R*pC*R;
 
% Continue greedy search if any parameters have been eliminated
%--------------------------------------------------------------------------
nelim  = full(sum(K(i,:)));
repeat = repeat & nelim;
if repeat
    fprintf('%i paramters eliminated in this step\n',nelim)
end
 
% % Show results
% % -------------------------------------------------------------------------
% spm_figure('Getwin','Graphics');
%  
% subplot(2,2,1)
% if length(K) > 32, plot(S,'k'), else, bar(S,'c'), end
% title('log-posterior','FontSize',16)
% xlabel('model','FontSize',12)
% ylabel('log-probability','FontSize',12)
% axis square
%  
% subplot(2,2,2)
% if length(K) > 32, plot(p,'k'), else, bar(p,'r'), end
% title('model posterior','FontSize',16)
% xlabel('model','FontSize',12)
% ylabel('probability','FontSize',12)
% axis square
 
% conditional estimates of selected model
% =========================================================================
Eq    = 0;
Pq    = 0;
for j = 1:N
    
    % Get priors and posteriors
    % ---------------------------------------------------------------------
    try, load(P{j}); catch, DCM = P{j}; end
    
    qE   = DCM.Ep;
    qC   = DCM.Cp;
    pE   = DCM.M.pE;
    pC   = DCM.M.pC;
    
    % Get posterior of selected model - rC
    % ---------------------------------------------------------------------
    [F Ep Cp] = spm_log_evidence(qE,qC,pE,pC,pE,rC);
    
    % Bayesian parameter average
    %----------------------------------------------------------------------
    Pp   = spm_inv(Cp);
    Eq   = Eq + Pp*spm_vec(Ep);
    Pq   = Pq + Pp;
 
 
    % Put reduced conditional estimates in DCM
    % =====================================================================
 
    % Bayesian inference and variance
    %----------------------------------------------------------------------
    sw   = warning('off','SPM:negativeVariance');
    Pp   = spm_unvec(1 - spm_Ncdf(DCM.T,abs(spm_vec(Ep)),diag(Cp)),Ep);
    Vp   = spm_unvec(diag(Cp),Ep);
    warning(sw);
    
    
    % Store parameter estimates
    %----------------------------------------------------------------------
    DCM.M.pC = rC;
    DCM.Ep   = Ep;
    DCM.Cp   = Cp;
    DCM.Pp   = Pp;
    DCM.Vp   = Vp;   
    
    % Save approximations to model evidence: negative free energy, AIC, BIC
    %----------------------------------------------------------------------
    evidence = spm_dcm_evidence(DCM);
    DCM.F    = F;
    DCM.AIC  = evidence.aic_overall;
    DCM.BIC  = evidence.bic_overall;
    
    
    %-Save optimised DCM
    %======================================================================
%     try
%         [pth, name] = fileparts(P{j});
%         P{j} = fullfile(pth,['DCM_opt_' name(4:end) '.mat']);
%     catch
%         P{j} = fullfile(pwd,['DCM_opt_' date '.mat']);
%     end
%     if spm_check_version('matlab','7') >= 0
%         save(P{j},'-V6','DCM','F','Ep','Cp');
%     else
%         save(P{j},'DCM','F','Ep','Cp');
%     end
    
end

function help_pushbutton_Callback(hObject, eventdata, handles)
web webhelp.html -browser


function click_on_cell(hObject,eventdata,handles)
global h_out

set(handles.figure1,'CurrentAxes',handles.estimated_model_axes);
imshow(imread('images/background_small.png'));

hold on
arrow1_pure = imread('images/arrow1.png');
arrows_est.arrow1 = imshow(arrow1_pure);
set(arrows_est.arrow1,'AlphaData',zeros(size(rgb2gray(arrow1_pure))));

hold on
arrow2_pure = imread('images/arrow2.png');
arrows_est.arrow2 = imshow(arrow2_pure);
set(arrows_est.arrow2,'AlphaData',zeros(size(rgb2gray(arrow2_pure))));

hold on
arrow3_pure = imread('images/arrow3.png');
arrows_est.arrow3 = imshow(arrow3_pure);
set(arrows_est.arrow3,'AlphaData',zeros(size(rgb2gray(arrow3_pure))));

hold on
arrow4_pure = imread('images/arrow4.png');
arrows_est.arrow4 = imshow(arrow4_pure);
set(arrows_est.arrow4,'AlphaData',zeros(size(rgb2gray(arrow4_pure))));

hold on
arrow5_pure = imread('images/arrow5.png');
arrows_est.arrow5 = imshow(arrow5_pure);
set(arrows_est.arrow5,'AlphaData',zeros(size(rgb2gray(arrow5_pure))));

hold on
arrow6_pure = imread('images/arrow6.png');
arrows_est.arrow6 = imshow(arrow6_pure);
set(arrows_est.arrow6,'AlphaData',zeros(size(rgb2gray(arrow6_pure))));

hold off

coordinates = eventdata;
row = coordinates.Indices(1);

DCM_A  = h_out{row}.DCM.Ep.A;
DCM_Ap = h_out{row}.DCM.Pp.A;
p = 0.9;             % Threshold for DCM
DCM_A(DCM_Ap<p) = 0;  % Thresholding
DCM_sA = logical(DCM_A.*(~eye(length(DCM_A)))); % Estimated model

set(handles.figure1,'CurrentAxes',handles.estimated_model_axes);
if DCM_sA(1,2) == 1
    set(arrows_est.arrow1,'AlphaData',rgb2gray(arrow1_pure)>80);
    hold on
end
if DCM_sA(1,3) == 1
    set(arrows_est.arrow2,'AlphaData',rgb2gray(arrow2_pure)>80);
    hold on
end
if DCM_sA(2,1) == 1
    set(arrows_est.arrow3,'AlphaData',rgb2gray(arrow3_pure)>80);
    hold on
end
if DCM_sA(2,3) == 1
    arrow4_pure = imread('images/arrow4.png');
    set(arrows_est.arrow4,'AlphaData',rgb2gray(arrow4_pure)>80);
    hold on
end
if DCM_sA(3,1) == 1
    set(arrows_est.arrow5,'AlphaData',rgb2gray(arrow5_pure)>80);
    hold on
end
if DCM_sA(3,2) == 1
    set(arrows_est.arrow6,'AlphaData',rgb2gray(arrow6_pure)>80);
    hold off
end



function save_results_as_matfile_pushbutton_Callback(hObject, eventdata, handles)
global h_out DCM_error_tab ratio fluct correl physiol str clk

readme = {'DCM simulation of BOLD in causal model'};

filename = ['DCM_sim_',str,'_',num2str(clk(1,4)),'_',num2str(clk(1,5)),...
    '_',num2str(fluct),'_',num2str(correl),'_',num2str(physiol)];

% Save results on disk.
save(filename,'readme','h_out','DCM_error_tab','ratio')
disp([''])
disp(['*******************************************************************************'])
disp(['Name of created file: ',filename, '.mat'])
msgbox(['Název uloženého souboru: ',filename, '.mat'],'Saved to the current folder','none')

function export_results_to_excel_pushbutton_Callback(hObject, eventdata, handles)
global DCM_error_tab fluct correl physiol str clk

filename = ['DCM_sim_',str,'_',num2str(clk(1,4)),'_',num2str(clk(1,5)),...
    '_',num2str(fluct),'_',num2str(correl),'_',num2str(physiol)];

xlswrite(filename, DCM_error_tab)
disp([''])
disp(['*******************************************************************************'])
disp(['Name of created file: ',filename, '.xls'])
msgbox(['Název uloženého souboru: ',filename, '.xls'],'Saved to the current folder','none')

% --- Executes on button press in plot_results_pushbutton.
function plot_results_pushbutton_Callback(hObject, eventdata, handles)
global DCM_error_tab ratio
data = DCM_error_tab;

figure
distributionPlot(data,'globalNorm',true,'colormap',1-gray(64),'histOpt',0,...
    'divFactor',[0:1:6],'xValues',ratio);

title(['Rozložení počtu chyb DCM odhadu pro ', num2str(ratio*100),...
    ' % šumového signálu a zvolený počet iterací']);

xlabel('Míra šumu');
ylabel('Počet chyb odhadu DCM');
colorbar

function handles = distributionPlot(varargin)
%DISTRIBUTIONPLOT creates violin plots for convenient visualization of multiple distributions
%
% SYNOPSIS: handles = distributionPlot(data,propertyName,propertyValue,...)
%           handles = distributionPlot(ah,...)
%
% INPUT data : m-by-nData array of values, or vector of grouped data (use
%           the 'groups' property to specify the grouping variable), or
%           cell array of length nData.
%           The cell array can either contain vectors with values, or
%           m-by-2 arrays with [bins,counts] if you want to determine the
%           histograms by yourself (m can be different between cell
%           elements). Note that arrays inside cells with any
%           other shape than m-by-2 are reshaped to vector an a warning is
%           thrown (DISTRIBUTIONPLOT:AUTORESHAPE).
%
%       DISTRIBUTIONPLOT accepts the following propertyName/propertyValue
%           pairs (all are optional):
%
%       distWidth :  width of distributions; ideally between 0 and 1.
%           1 means that adjacent distributions might touch. Default: 0.9
%       variableWidth : If true, the width of the distribution changes,
%           reflecting the shape of the histogram of the data. If false,
%           the distribution is only encoded by color levels. Default: true
%       color : uniform coloring of histograms. Supply either a color
%           string ('r'), or a truecolor vector ([1 0 0]). Use a
%           cell array of length nData to specify one color per
%           distribution. Default: 'k'
%           If variableWidth is set to false, a colormap is generated that
%           goes from white to the chose color (or from black, if
%           invert==true).
%           If both 'color', and 'colormap' are specified, 'colormap' takes
%           precedence.
%       colormap : colormap used to describe the distribution (first row
%           corresponds to bins with least data, last row corresponds to
%           bins with most data (invert the grayscale colormap to have
%           black indicate the most data).
%           Supply a cell array of length nData to color distributions
%           individually. Note that using multiple colormaps means that
%           the colorbar doesn't contain much useful information.
%           Default: []
%           Colormap will index into the figure colormap, which will be
%           modified by distributionPlot. This is done to allow editing the
%           distributions in e.g. Adobe Illustrator.
%           If both 'color', and 'colormap' are specified, 'colormap' takes
%           precedence.
%       globalNorm : normalization for bin width (x-direction)
%           0 : every histogram is normalized individually so that the
%               maximum bin width is equal to distWidth. This is best
%               suited to comparing distribution shapes. Default.
%           1 : histograms are normalized such that equal bin width
%               reports equal numbers of counts per bin.
%           2 : histograms are normalized so that the relative areas
%               covered by the histograms reflect the relative total number
%               of data points.
%           3 : histograms areas are normalized so that relative densities
%               are the same across histograms. Thus, if
%               data = {rand(100,1),rand(500,1)},
%               then
%               distributionPlot(data,'globalNorm',2,'histOpt',0,'divFactor',10)
%               shows the left histogram 5x as wide as the right, while
%               distributionPlot(data,'globalNorm',3,'histOpt',0,'divFactor',10)
%               displays both histograms equally wide, since each bin
%               contains ~10% of the data.
%           Options 1 and 2 produce similar results if the bins are spaced
%           equally for the distributions. Options 0 and 3 produce similar
%           results if the data are drawn from the same distributions.
%           Note that colormaps currently always report the number of data
%           points per bin; 'globalNorm' only applies to the distribution
%           shape.
%
%       groups : grouping variable for grouped data. Grouping will be
%                   resolved by calling grp2idx, and unless xNames have
%                   been supplied, group names determine the x-labels.
%                   If the grouping variable is numeric, group labels also
%                   determine x-values, unless the parameter xValues has
%                   been specified.
%       histOpt : histogram type to plot
%                   0 : use hist command (no smoothing, fixed number of
%                       bins)
%                   1 : smoothened histogram using ksdensity with
%                       Normal kernel. Default.
%                   1.1: smoothened histogram using ksdensity where the
%                       kernel is robustly estimated via histogram.m.
%                       Normal kernel.
%                   2 : histogram command (no smoothing, automatic
%                       determination of thickness (y-direction) of bins)
%       divFactor : Parameter dependent on histOpt. If...
%                   histOpt == 0: divFactor = # of bins. Default: 25.
%                       Alternatively, pass a vector which will be
%                       interpreted as bin centers.
%                   histOpt == 1: divFactor decides by how much the default
%                       kernel-width is multiplied in order to avoid an
%                       overly smooth histogram. Default: 1/2
%                   histOpt == 2: divFactor decides by how much the
%                       automatic bin width is multiplied in order to have
%                       more (<1) or less (>1) detail. Default: 1
%       addSpread : if 1, data points are plotted with plotSpread.
%                   distWidth is ideally set to 0.95
%                   This option is not available if the data is supplied as
%                   histograms.
%                   Please download plotSpread.m separately from the File
%                   Exchange using the link in the remarks
%       showMM : if 1, mean and median are shown as red crosses and
%                green squares, respectively. This is the default
%                2: only mean
%                3: only median
%                4: mean +/- standard error of the mean (no median)
%                5: mean +/- standard deviation (no median)
%                6: draw lines at the 25,50,75 percentiles (no mean)
%                0: plot neither mean nor median
%       xValues: x-coordinate where the data should be plotted.
%                If xValues are given, "distWidth" is scaled by the median
%                difference between adjacent (sorted) x-values. Note that
%                this may lead to overlapping distributions. Default:
%                1:nData
%       xNames : cell array of length nData containing x-tick names
%               (instead of the default '1,2,3')
%       xMode  : if 'auto', x-ticks are spaced automatically. If 'manual',
%                there is a tick for each distribution. If xNames is
%                provided as input, xMode is forced to 'manual'. Default:
%                'manual'.
%          NOTE: SPECIFYING XNAMES OR XVALUES OR XMODE WILL ERASE PREVIOUS
%                LABELS IF PLOTTING INTO EXISTING AXES
%       yLabel : string with label for y-axis. Default : ''
%                If empty and data is histograms, ylabel is set to 'counts'
%       invert : if 1, axes color is changed to black, and colormap is
%                   inverted.
%       histOri: Orientation of histogram. Either 'center', 'left', or
%                'right'. With 'left' or 'right', the left or right half of
%                the standard violin plot is shown. Has no effect if
%                variableWidth is false. Default: center
%       xyOri  : orientation of axes. Either 'normal' (=default), or
%                'flipped'. If 'flipped', the x-and y-axes are switched, so
%                that violin plots are horizontal. Consequently,
%                axes-specific properties, such as 'yLabel' are applied to
%                the other axis.
%       widthDiv : 1-by-2 array with [numberOfDivisions,currentDivision]
%                widthDiv allows cutting the stripe dedicated to a single
%                distribution into multible bands, which can be filled with
%                sequential calls to distributionPlot. This is one way
%                to compare two (or more) sequences of distributions. See
%                example below.
%       ah : axes handle to plot the distributions. Default: gca
%
% OUTPUT handles : 1-by-4 cell array with patch-handles for the
%                  distributions, plot handles for mean/median, the
%                  axes handle, and the plotSpread-points handle
%
%
% EXAMPLES
%         %--Distributions contain more information than boxplot can capture
%         r = rand(1000,1);
%         rn = randn(1000,1)*0.38+0.5;
%         rn2 = [randn(500,1)*0.1+0.27;randn(500,1)*0.1+0.73];
%         rn2=min(rn2,1);rn2=max(rn2,0);
%         figure
%         ah(1)=subplot(3,4,1:2);
%         boxplot([r,rn,rn2])
%         ah(2)=subplot(3,4,3:4);
%         distributionPlot([r,rn,rn2],'histOpt',2); % histOpt=2 works better for uniform distributions than the default
%         set(ah,'ylim',[-1 2])
%
%         %--- additional options  
%
%         data = [randn(100,1);randn(50,1)+4;randn(25,1)+8];
%         subplot(3,4,5)
%
%         %--- defaults
%         distributionPlot(data); 
%         subplot(3,4,6)
%
%         %--- show density via custom colormap only, show mean/std,
%         distributionPlot(data,'colormap',copper,'showMM',5,'variableWidth',false) 
%         subplot(3,4,7:8)
%
%         %--- auto-binwidth depends on # of datapoints; for small n, plotting the data is useful
%         % note that this option requires the additional installation
%         % of plotSpread from the File Exchange (link below)
%         distributionPlot({data(1:5:end),repmat(data,2,1)},'addSpread',true,'showMM',false,'histOpt',2) 
%
%         %--- show quantiles
%         subplot(3,4,9),distributionPlot(randn(100,1),'showMM',6)
%
%         %--- horizontal orientation
%         subplot(3,4,10:11),
%         distributionPlot({chi2rnd(3,1000,1),chi2rnd(5,1000,1)},'xyOri','flipped','histOri','right','showMM',0),
%         xlim([-3 13])
%
%         %--- compare distributions side-by-side (see also example below)
%         % plotting into specified axes will throw a warning that you can
%         % turn off using " warning off DISTRIBUTIONPLOT:ERASINGLABELS "
%         ah = subplot(3,4,12);
%         subplot(3,4,12),distributionPlot(chi2rnd(3,1000,1),'histOri','right','color','r','widthDiv',[2 2],'showMM',0)
%         subplot(3,4,12),distributionPlot(chi2rnd(5,1000,1),'histOri','left','color','b','widthDiv',[2 1],'showMM',0)
%
%         %--Use globalNorm to generate meaningful colorbar
%         data = {randn(100,1),randn(500,1)};
%         figure
%         distributionPlot(data,'globalNorm',true,'colormap',1-gray(64),'histOpt',0,'divFactor',[-5:0.5:5])
%         colorbar
%
%         %--Use widthDiv to compare two series of distributions
%         data1 = randn(500,5);
%         data2 = bsxfun(@plus,randn(500,5),0:0.1:0.4);
%         figure
%         distributionPlot(data1,'widthDiv',[2 1],'histOri','left','color','b','showMM',4)
%         distributionPlot(gca,data2,'widthDiv',[2 2],'histOri','right','color','k','showMM',4)
%
%         %--Christmas trees!
%           x=meshgrid(1:10,1:10);
%           xx = tril(x);
%           xx = xx(xx>0);
%           figure
%           hh=distributionPlot({xx,xx,xx},'color','g','addSpread',1,'histOpt',2,'showMM',0);
%           set(hh{4}{1},'color','r','marker','o')
% END
%
% REMARKS To show distributions as clouds of points (~beeswarm plot),
%         and/or to use the option "addSpread", please download the
%         additional function plotSpread.m from the File Exchange
%         http://www.mathworks.com/matlabcentral/fileexchange/37105-plot-spread-points-beeswarm-plot
%         
%         I used to run ksdensity with the Epanechnikov kernel. However,
%         for integer data, the shape of the kernel can produce peaks
%         between the integers, which is not ideal (use histOpt=2 for
%         integer valued data).
%
%         A previous iteration of distributionPlot used the input
%         specifications below. They still work to ensure backward
%         compatibility, but are no longer supported or updated.
%           handles = distributionPlot(data,distWidth,showMM,xNames,histOpt,divFactor,invert,addSpread,globalNorm)
%           where distWidth of 1 means that the maxima
%           of  two adjacent distributions might touch. Negative numbers
%           indicate that the distributions should have constant width, i.e
%           the density is only expressed through greylevels.
%           Values between 1 and 2 are like values between 0 and 1, except
%           that densities are not expressed via graylevels. Default: 1.9
%
%
% SEE ALSO histogram, ksdensity, plotSpread, boxplot, grp2idx
%

% created with MATLAB ver.: 7.6.0.324 (R2008a) on Windows_NT
%
% created by: Jonas Dorn; jonas.dorn@gmail.com
% DATE: 08-Jul-2008
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%====================================
%% TEST INPUT
%====================================

% set defaults
def.xNames = [];
def.showMM = 1;
def.distWidth = 0.9;
def.histOpt = 1;
def.divFactor = [25,2,1];
def.invert = false;
def.colormap = [];
def.color = 'k';
def.addSpread = false;
def.globalNorm = false;
def.variableWidth = true;
def.groups = [];
def.yLabel = '';
def.xValues = '';
def.xMode = 'manual';
def.histOri = 'center';
def.xyOri = 'normal';
def.widthDiv = [1 1];
isHistogram = false; %# this parameter is not set by input


if nargin == 0 || isempty(varargin{1})
    error('not enough input arguments')
end

% check for axes handle
if ~iscell(varargin{1}) && isscalar(varargin{1}) == 1 && ...
        ishandle(varargin{1}) && strcmp(get(varargin{1},'Type'),'axes')
    ah = varargin{1};
    data = varargin{2};
    varargin(1:2) = [];
    newAx = false;
    
    
else
    ah = gca;
    data = varargin{1};
    varargin(1) = [];
    newAx = true;
end

% check for current axes limits. Set NaN if the axes have no children
% yet - we need that in case we're building a complicated set of
% distributions
if ~isempty(get(ah,'children'))
    xAxLim = xlim;
    yAxLim = ylim;
else
    [xAxLim,yAxLim] = deal([NaN NaN]);
end

fh = get(ah,'Parent');

% check data. If not cell, convert
if ~iscell(data)
    [nPoints,nData] = size(data);
    data = mat2cell(data,nPoints,ones(nData,1));
else
    % get nData
    data = data(:);
    nData = length(data);
    % make sure all are vectors
    badCol = ~cellfun(@isvector,data) & ~cellfun(@isempty,data);
    if any(badCol)
        nCols = cellfun(@(x)(size(x,2)),data(badCol));
        if all(nCols==2)
            % bins,counts
            isHistogram = true;
        else
            warning('DISTRIBUTIONPLOT:AUTORESHAPE',...
                'Elements %s of the cell array are not vectors. They will be reshaped automatically',...
                num2str(find(badCol)'));
            data(badCol) = cellfun(@(x)(x(:)),data(badCol),'UniformOutput',false);
        end
    end
end

parserObj = inputParser;
parserObj.FunctionName = 'distributionPlot';
stdWidth = 1; % scaling parameter for variableWidth with uneven x-values
% check whether we're dealing with pN/pV or straight arguments
if ~isempty(varargin) && ~ischar(varargin{1}) && ~isstruct(varargin{1})
    % use old format
    % distWidth,showMM,xNames,histOpt,divFactor,invert,addSpread,globalNorm
    def.distWidth = 1.9;
    parserObj.addOptional('distWidth',def.distWidth);
    parserObj.addOptional('showMM',def.showMM);
    parserObj.addOptional('xNames',def.xNames);
    parserObj.addOptional('histOpt',def.histOpt);
    parserObj.addOptional('divFactor',def.divFactor);
    parserObj.addOptional('invert',def.invert);
    parserObj.addOptional('addSpread',def.addSpread);
    parserObj.addOptional('globalNorm',def.globalNorm);
    parserObj.addOptional('groups',def.groups);
    parserObj.addOptional('yLabel',def.yLabel);
    parserObj.addOptional('color',def.color);
    
    
    parserObj.parse(varargin{:});
    opt = parserObj.Results;
    % fill in defaults that are not supported in the old version of the
    % code
    opt.colormap = [];
    opt.variableWidth = true;
    opt.histOri = 'center';
    opt.xValues = [];
    opt.xMode = 'auto';
    opt.xyOri = 'normal';
    opt.widthDiv = [1 1];
    
    % overwrite empties with defaults - inputParser considers empty to be a
    % valid input.
    fnList = fieldnames(opt);
    for fn = fnList'
        if isempty(opt.(fn{1}))
            opt.(fn{1}) = def.(fn{1});
        end
    end
    
    
    % fix a few parameters
    if opt.distWidth > 1
        opt.distWidth = opt.distWidth - 1;
    else
        opt.colormap = 1-gray(128);
    end
    if opt.distWidth < 0
        opt.variableWidth = false;
        opt.distWidth = abs(opt.distWidth);
    end
    
    if ~isempty(opt.xNames)
        opt.xMode = 'manual';
    end
    
    
else
    defNames = fieldnames(def);
    for dn = defNames(:)'
        parserObj.addParamValue(dn{1},def.(dn{1}));
    end
    
    
    parserObj.parse(varargin{:});
    opt = parserObj.Results;
    
    % if groups: deal with data
    if ~isempty(opt.groups)
        [idx,labels,vals] = grp2idx(opt.groups);
        % convert data to cell array
        data = accumarray(idx,data{1},[],@(x){x});
        nData = length(data);
        % if not otherwise provided, use group labels for xnames
        if isempty(opt.xNames)
            opt.xNames = labels;
            if ~iscell(opt.xNames)
                opt.xNames = num2cell(opt.xNames);
            end
        end
        if isnumeric(vals) && isempty(opt.xValues)
            opt.xValues = vals;
        end
        
    end
    
    if ~ischar(opt.xyOri) || ~any(ismember(opt.xyOri,{'normal','flipped'}))
        error('option xyOri must be either ''normal'' or ''flipped'' (is ''%s'')',opt.xyOri);
    end
    
    
    
    
end
% common checks

% default x-values: 1:n
if isempty(opt.xValues)
    opt.xValues = 1:nData;
elseif length(opt.xValues) ~= nData
    error('please supply as many x-data values as there are data entries')
elseif length(opt.xValues) > 1 % only check for scale if more than 1 value
    % scale width
    stdWidth = median(diff(sort(opt.xValues)));
    opt.distWidth = opt.distWidth * stdWidth;
end

if ~isscalar(opt.divFactor) && length(opt.divFactor) == 3 && all(opt.divFactor==def.divFactor)
    opt.divFactor = opt.divFactor(floor(opt.histOpt)+1);
end
if isHistogram
    opt.histOpt = 99;
    if isempty(opt.yLabel)
        opt.yLabel = 'counts';
    end
end



% check colors/colormaps: do we need to expand colormap?
if ~iscell(opt.colormap)
    opt.colormap = {opt.colormap};
end
if ~iscell(opt.color)
    opt.color = {opt.color};
end
for iColor = 1:length(opt.color)
    if ischar(opt.color{iColor})
        opt.color{iColor} = colorCode2rgb(opt.color{iColor});
    end
end

% expand - if only single colormap specified, we expand only once
if ~opt.variableWidth
    missingColormaps = find(cellfun(@isempty,opt.colormap));
    for iMissing = missingColormaps(:)'
        
        endColor = opt.color{max(iMissing,length(opt.color))};
        % normally, we go from white to color
        cmap = zeros(128,3);
        for rgb = 1:3
            cmap(:,rgb) = linspace(1,endColor(rgb),128);
        end
        opt.colormap{iMissing} = cmap;
        
    end
end

% if we have colormaps, we need to create a master which we add to the
% figure. Invert if necessary, and expand the cell array to nData
colormapLength = cellfun(@(x)size(x,1),opt.colormap);
if any(colormapLength>0)
    
    colormap = cat(1,opt.colormap{:});
    if opt.invert
        colormap = 1-colormap;
    end
    set(fh,'Colormap',colormap)
    if length(opt.colormap) == 1
        opt.colormap = repmat(opt.colormap,nData,1);
        colormapLength = repmat(colormapLength,nData,1);
        colormapOffset = zeros(nData,1);
        singleMap = true;
    else
        colormapOffset = [0;cumsum(colormapLength(1:end-1))];
        singleMap = false;
    end
    
else
    
    colormapLength = zeros(nData,1);
    if length(opt.color) == 1
        opt.color = repmat(opt.color,nData,1);
    end
    if opt.invert
        opt.color = cellfun(@(x)1-x,opt.color,'uniformOutput',false);
    end
end


% set hold on
holdState = get(ah,'NextPlot');
set(ah,'NextPlot','add');

% if new axes: invert
if newAx && opt.invert
    set(ah,'Color','k')
end

%===================================



%===================================
%% PLOT DISTRIBUTIONS
%===================================

% assign output
hh = NaN(nData,1);
[m,md,sem,sd] = deal(nan(nData,1));
if opt.showMM == 6
    md = nan(nData,3,3); % md/q1/q3, third dim is y/xmin/xmax
end

% get base x-array
% widthDiv is a 1-by-2 array with
% #ofDivs, whichDiv
% The full width (distWidth) is split into
% #ofDivs; whichDiv says which "stripe" is active
xWidth = opt.distWidth/opt.widthDiv(1);
xMin = -opt.distWidth/2;
xLow = xMin + xWidth * (opt.widthDiv(2)-1);
xBase = [-xWidth;xWidth;xWidth;-xWidth]/2;
xOffset = xLow + xWidth/2;

% b/c of global norm: loop twice
plotData = cell(nData,2);

% loop through data. Prepare patch input, then draw patch into gca
for iData = 1:nData
    currentData = data{iData};
    % only plot if there is some finite data
    if ~isempty(currentData(:)) && any(isfinite(currentData(:)))
        
        switch floor(opt.histOpt)
            case 0
                % use hist
                [xHist,yHist] = hist(currentData,opt.divFactor);
                
            case 1
                % use ksdensity
                
                if opt.histOpt == 1.1
                    % use histogram to estimate kernel
                    [dummy,x] = histogram(currentData); %#ok<ASGLU>
                    if length(x) == 1
                        % only one value. Make fixed distribution
                        dx = 0.1;
                        yHist = x;
                        xHist = sum(isfinite(currentData));
                    else
                        dx = x(2) - x(1);
                    
                    % make sure we sample frequently enough
                    x = min(x)-dx:dx/3:max(x)+dx;
                    [xHist,yHist] = ksdensity(currentData,x,'kernel','normal','width',dx/(1.5*opt.divFactor));
                    end
                else
                    
                    % x,y are switched relative to normal histogram
                    [xHist,yHist,u] = ksdensity(currentData,'kernel','normal');
                    % take smaller kernel to avoid over-smoothing
                    if opt.divFactor ~= 1
                        [xHist,yHist] = ksdensity(currentData,'kernel','normal','width',u/opt.divFactor);
                    end
                end
                
                % modify histogram such that the sum of bins (not the
                % integral under the curve!) equals the total number of
                % observations, in order to be comparable to hist
                xHist = xHist/sum(xHist)*sum(isfinite(currentData));
                
            case 2
                % use histogram - bar heights are counts as in hist
                [xHist,yHist] = histogram(currentData,opt.divFactor,0);
            case 99
                % bins,counts already supplied
                xHist = currentData(:,2)';
                yHist = currentData(:,1)';
        end
        plotData{iData,1} = xHist;
        plotData{iData,2} = yHist;
    end
end

goodData = find(~cellfun(@isempty,plotData(:,1)));
% get norm
switch opt.globalNorm
    case 3
        % #3 normalizes relative densities
        xNorm(goodData) = cellfun(@(x)min(diff(x)),plotData(goodData,2));
        xNorm(goodData) = xNorm(goodData) .* cellfun(@sum,plotData(goodData,1))';
        maxNorm(goodData) = cellfun(@max,plotData(goodData,1));
        xNorm(goodData) = xNorm(goodData)*max(maxNorm(goodData)./xNorm(goodData));
        
    case 2
        % #2 should normalize so that the integral of the
        % different histograms (i.e. area covered) scale with the
        % respective sum of counts across all bins. Requires evenly spaced
        % histograms at the moment
        xNorm(goodData) = cellfun(@(x)min(diff(x)),plotData(goodData,2));
        maxNorm(goodData) = cellfun(@max,plotData(goodData,1));
        xNorm(goodData) = xNorm(goodData)*max(maxNorm(goodData)./xNorm(goodData));
    case 1
        xNorm(goodData) = max(cat(2,plotData{:,1}));
    case 0
        xNorm(goodData) = cellfun(@max,plotData(goodData,1));
end


for iData = goodData'
    
    % find current data again
    currentData = data{iData};
    
    xHist = plotData{iData,1};
    yHist = plotData{iData,2};
    
    % find y-step
    dy = min(diff(yHist));
    if isempty(dy)
        dy = 0.1;
    end
    
    % create x,y arrays
    nPoints = length(xHist);
    xArray = repmat(xBase,1,nPoints);
    yArray = repmat([-0.5;-0.5;0.5;0.5],1,nPoints);
    
    
    % x is iData +/- almost 0.5, multiplied with the height of the
    % histogram
    if opt.variableWidth
        
        
        tmp = xArray.*repmat(xHist,4,1)./xNorm(iData);
        
        switch opt.histOri
            case 'center'
                % we can simply use xArray
                xArray = tmp;
            case 'right'
                % shift everything to the left
                delta = tmp(1,:) - xArray(1,:);
                xArray = bsxfun(@minus,tmp,delta);
            case 'left'
                % shift everything to the right
                delta = tmp(1,:) - xArray(1,:);
                xArray = bsxfun(@plus,tmp,delta);
        end
        
        xArray = xArray + opt.xValues(iData);
        
    else
        xArray = xArray + iData;
    end
    
    % add offset (in case we have multiple widthDiv)
    xArray = xArray + xOffset;
    
    
    % yData is simply the bin locations
    yArray = repmat(yHist,4,1) + dy*yArray;
    
    % add patch
    vertices = [xArray(:),yArray(:)];
    faces = reshape(1:numel(yArray),4,[])';
    
    if colormapLength(iData) == 0
        colorOpt = {'FaceColor',opt.color{iData}};
    else
        % calculate index into colormap
        if singleMap
            % use scaled mapping so that colorbar is meaningful
            if opt.globalNorm > 0
                colorOpt = {'FaceVertexCData',xHist','CDataMapping','scaled','FaceColor','flat'};
            else
                colorOpt = {'FaceVertexCData',xHist'/xNorm(iData),'CDataMapping','scaled','FaceColor','flat'};
            end
            
        else
            idx = round((xHist/xNorm(iData))*(colormapLength(iData)-1))+1;
            colorOpt = {'FaceVertexCData',idx'+colormapOffset(iData),'CDataMapping','direct','FaceColor','flat'};
        end
    end
    
    
    switch opt.xyOri
        case 'normal'
            hh(iData)= patch('Vertices',vertices,'Faces',faces,'Parent',ah,colorOpt{:},'EdgeColor','none');
        case 'flipped'
            hh(iData)= patch('Vertices',vertices(:,[2,1]),'Faces',faces,'Parent',ah,colorOpt{:},'EdgeColor','none');
    end
    
    if opt.showMM > 0
        if isHistogram
            [m(iData),sem(iData)] = weightedStats(currentData(:,1),currentData(:,2),'w');
            sd(iData) = sem(iData) * sqrt(sum(currentData(:,2)));
            % weighted median: where we're at middle weight
            % may need some tweaking
            goodCurrentData = sortrows(currentData(all(isfinite(currentData),2),:),1);
            weightList = cumsum(goodCurrentData(:,2));
            weightList = weightList / weightList(end);
            md(iData) = goodCurrentData(find(weightList>0.5,1,'first'),1);
        else
            m(iData) = nanmean(currentData);
            md(iData) = nanmedian(currentData);
            sd(iData) = nanstd(currentData);
            sem(iData) = sd(iData)/sqrt(sum(isfinite(currentData)));
        end
        
        if opt.showMM == 6
            % read quantiles - "y"-value, plus x-start-stop
            % re-use md array which allows using a loop below instead of
            % lots of copy-paste
            % md array is md/q1/q3, with third dimension y/xmin/xmax
            
            md(iData,2,1) = prctile(currentData,25);
            md(iData,3,1) = prctile(currentData,75);
            
            for qq = 1:3
                % find corresponding y-bin
                yLoc =  repmat(...
                    any(yArray>md(iData,qq,1),1) & any(yArray<=md(iData,qq,1),1),...
                    [4 1]);
                % look up corresponding x-values. Note that there is a bit
                % of a risk that the line will be exactly between two very
                % different bins - but if we make the line longer, it will
                % be ugly almost all the time
                md(iData,qq,2) = min( xArray( yLoc ) );
                md(iData,qq,3) = max( xArray( yLoc ) );
            end
            
        end
    end
end % loop

sh = [];
if opt.addSpread
    if isHistogram
        disp('Option addSpread is unavailable if data is supplied as histograms. Call plotSpread separately')
    else
        % add spread
        try
            sh = plotSpread(ah,data,'xValues',opt.xValues,'xyOri',opt.xyOri);
            set(sh{1},'color',[0,128,255]/255);
        catch me
            if strcmp(me.identifier,'MATLAB:UndefinedFunction')
                error('plotSpread not found. Please download it from the Matlab File Exchange')
            else
                rethrow(me)
            end
        end
    end
end

mh = [];mdh=[];
if opt.showMM
    % plot mean, median. Mean is filled red circle, median is green square
    % I don't know of a very clever way to flip xy and keep everything
    % readable, thus it'll be copy-paste
    switch opt.xyOri
        case 'normal'
            if any(opt.showMM==[1,2])
                mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12);
            end
            if any(opt.showMM==[1,3])
                mdh = plot(ah,opt.xValues+xOffset,md,'sg','MarkerSize',12);
            end
            if opt.showMM == 4
                mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12);
                mdh = myErrorbar(ah,opt.xValues+xOffset,m,sem);
            end
            if opt.showMM == 5
                mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12);
                mdh = myErrorbar(ah,opt.xValues+xOffset,m,sd);
            end
            if opt.showMM == 6
                mdh(1,:) = plot(ah,squeeze(md(:,1,2:3))',repmat(md(:,1,1)',2,1),'color','r','lineWidth',2);%,'lineStyle','--');
                mdh(2,:) = plot(ah,squeeze(md(:,2,2:3))',repmat(md(:,2,1)',2,1),'color','r','lineWidth',1);%,'lineStyle','--');
                mdh(3,:) = plot(ah,squeeze(md(:,3,2:3))',repmat(md(:,3,1)',2,1),'color','r','lineWidth',1);%,'lineStyle','--');
            end
        case 'flipped'
            if any(opt.showMM==[1,2])
                mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12);
            end
            if any(opt.showMM==[1,3])
                mdh = plot(ah,md,opt.xValues+xOffset,'sg','MarkerSize',12);
            end
            if opt.showMM == 4
                mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12);
                mdh = myErrorbar(ah,m,opt.xValues+xOffset,[sem,NaN(size(sem))]);
            end
            if opt.showMM == 5
                mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12);
                mdh = myErrorbar(ah,m,opt.xValues+xOffset,[sd,NaN(size(sd))]);
            end
            if opt.showMM == 6
                mdh(1,:) = plot(ah,repmat(md(:,1,1)',2,1),squeeze(md(:,1,2:3))','color','r','lineWidth',2);%,'lineStyle','--');
                mdh(2,:) = plot(ah,repmat(md(:,2,1)',2,1),squeeze(md(:,2,2:3))','color','r','lineWidth',1);%,'lineStyle','--');
                mdh(3,:) = plot(ah,repmat(md(:,3,1)',2,1),squeeze(md(:,3,2:3))','color','r','lineWidth',1);%,'lineStyle','--');
            end
    end
end

% find extents of x-axis (or y-axis, if flipped)
minX = min(opt.xValues)-stdWidth;
maxX = max(opt.xValues)+stdWidth;

if ~isnan(xAxLim(1))
    % we have previous limits
    switch opt.xyOri
        case 'normal'
            minX = min(minX,xAxLim(1));
            maxX = max(maxX,xAxLim(2));
        case 'flipped'
            minX = min(minX,yAxLim(1));
            maxX = max(maxX,yAxLim(2));
    end
end


% if ~empty, use xNames
switch opt.xyOri
    case 'normal'
        switch opt.xMode
            case 'manual'
                if newAx == false
                    warning('DISTRIBUTIONPLOT:ERASINGLABELS','Plotting into an existing axes and specifying labels will erase previous labels')
                end
                set(ah,'XTick',opt.xValues);
                if ~isempty(opt.xNames)
                    set(ah,'XTickLabel',opt.xNames)
                end
            case 'auto'
                % no need to do anything
        end
        if ~isempty(opt.yLabel)
            ylabel(ah,opt.yLabel);
        end
        % have plot start/end properly
        xlim([minX,maxX])
    case 'flipped'
        switch opt.xMode
            case 'manual'
                if newAx == false
                    warning('DISTRIBUTIONPLOT:ERASINGLABELS','Plotting into an existing axes and specifying labels will erase previous labels')
                end
                set(ah,'YTick',opt.xValues);
                if ~isempty(opt.xNames)
                    set(ah,'YTickLabel',opt.xNames)
                end
            case 'auto'
                % no need to do anything
        end
        if ~isempty(opt.yLabel)
            xlabel(ah,opt.yLabel);
        end
        % have plot start/end properly
        ylim([minX,maxX])
end


%==========================


%==========================
%% CLEANUP & ASSIGN OUTPUT
%==========================

if nargout > 0
    handles{1} = hh;
    handles{2} = [mh;mdh];
    handles{3} = ah;
    handles{4} = sh;
end

set(ah,'NextPlot',holdState);

function rgbVec = colorCode2rgb(c)
%COLORCODE2RGB converts a color code to an rgb vector
%

% SYNOPSIS rgbVec = colorCode2rgb(c)
%
% INPUT c : color code
%       The following colors are supported:
%        'y' 'yellow'
%        'm' 'magenta'
%        'c' 'cyan'
%        'r' 'red'
%        'g' 'green'
%        'b' 'blue'
%        'w' 'white'
%        'k' 'black'
% 
% OUTPUT rgbVec : vector with the rgb value
%
% EXAMPLE 
%    rgb = colorCode2rgb('r')
%    rgb =
%          [1 0 0]

switch c
case {'y','yellow'}, rgbVec = [1,1,0];
case {'m','magenta'}, rgbVec = [1,0,1];
case {'c','cyan'}, rgbVec = [0,1,1];
case {'r','red'}, rgbVec = [1,0,0];
case {'g','green'}, rgbVec = [0,1,0];
case {'b','blue'}, rgbVec = [0,0,1];
case {'w','white'}, rgbVec = [1,1,1];
case {'k','black'}, rgbVec = [0,0,0];
otherwise, error('unknown color code %s',c)
end;
