supFunSim.org

#+TITLE: supFunSim.org
#+AUTHOR: Jan Nikadon, Tomasz Piotrowski, Dania Gutiérrez
#+EMAIL: nikadon [at] gmail [dot] com

#+STARTUP: hideall
#+PROPERTY: eval no-export
#+PROPERTY: exports both


* DONE Preliminaries
  CLOSED: [2017-07-19 Wed 14:40]
** DONE [i01] Introduction

   #+BEGIN_QUOTE

SYNOPSIS: =supFunSim= provides a comprehensive simulations framework
for solving the /forward/ and the /inverse/ problems in
electroencephalography (EEG).

   #+END_QUOTE

   =supFunSim= (provided as [[./supFunSim.org][supFunSim.org]] file) contains both ordinary
   human language and /(discrete) computing machine/ code that are
   available in this document in the form of /plain text/.  This text
   tells a story about the reconstruction of neuronal activity
   originating from discrete nodes spread among separate brain
   regions.  This reconstruction is based on the EEG signal.

   First neuronal activity is realistically simulated (we also refer
   to that activity as time series in /source-space/).  The solution
   to a /forward problem/ (FP) is based on (1) the geometry and
   conductivity of head compartments and (2) the position of
   electrodes on the scalp.  This solution allows for generation of
   simulated EEG signals (i.e., time-series in /sensor-space/ that
   result from the activity of neuronal sources).

   Later, sensor-space time-series are used as an input for a variety
   of spatial filtering methods implemented here. This produces
   reconstruction (OR estimation) of the original activity of the
   neuronal sources (i.e., each spatial filter constitutes (OR
   provides) an /inverse problem/ solution).  Next, filtering output
   is compared with the original source-space time-series in order to
   estimate reconstruction errors of a given spatial filter.

   Some of the filters implemented here are known in the literature
   and some are novel generalizations of the aforementioned.  A number
   of parameters is available for manipulation in our simulations
   setup, these are, among others:
   - number of considered sources,
   - signal to noise ratios for a number of noise types (described below).
   We also provide a number of reconstruction error evaluation methods
   (that are relevant for the present neuroscience research).

** DONE [i02] Acknowledgments

   Herein we adapt number of solutions provided in toolboxes developed
   by others, namely:
   - [[http://www.fieldtriptoolbox.org/][FieldTrip]],
   - [[http://www.cs.tut.fi/~gomezher/projects/eeg/mvarica.htm][MVAR-ICA]],
   - [[http://climate-dynamics.org/software/#arfit][ARfit]] and
   - [[https://sccn.ucsd.edu/eeglab/][EEGlab]]

   Also, please note that creation of this document is an ongoing
   process as we intend to expand it (therefore we are not seeking for
   it to be complete in any near future).

** DOIN [i03] Technicalities
*** Introduction to *GNU Emacs* and =org-mode=

    This document was written using [[https://www.gnu.org/software/emacs/][GNU Emacs]] in plain text exercising
    advantages of [[http://www.orgmode.org][org-mode]] and [[http://orgmode.org/worg/org-contrib/babel/][(org-)Babel]].  This Emacs mode allows
    for the application of a *[[https://en.wikipedia.org/wiki/Literate_programming][literate]]* approach to programming.
    Briefly speaking, the /source code blocks/ are interspersed with
    /ordinary human language blocks/ that provide explanations and
    some insights into governing dynamics and intrinsic mechanics of
    the code.

    =org-mode= allows for easy export of whole documents to HTML,
    \LaTeX, PDF and other formats.  Therefore, various versions of
    this document may also be available, but the plain text version is
    *favorable* for studying the simulations setup (as it provides the
    most interactive encounter with the code).  =org-mode= also allows
    for latex formatted equations, e.g.,

    $e^{i \pi} + 1 = 0$

    GNU Emacs (with =org-mode= and ~Babel~ enabled) provides an
    environment in which this document can be read and at the same
    time the code snippets can be executed (evaluated) within the same
    environment and the code execution results can be incorporated in
    this document as a plain text, figures or in another form
    (provided that, apart from GNU Emacs, the given machine can
    execute [[https://www.mathworks.com/products/matlab.html][MATLAB]] or [[https://www.gnu.org/software/octave/][GNU Octave]] scripts).

    In order to execute all of the simulations steps in a single
    run (or multiple runs) please proceed to the section [[*%5Bs07%5D%20BATCH][BATCH]] (or
    [[*%5Bs08%5D%20LOOP][LOOP]] respectively). Those sections extensively exert on the
    *tangling*, a feature of the =org-mode= that allows for source
    code blocks content to be written to separate files that can be
    later executed (e.g., using MATLAB).

    For this simulations files that are tangled from this
    =org-mode= document are stored in the same [[./][directory]] as the
    =supFunSim.org= file and have names that follow the pattern
    =tg_*=.

    Simulations are hosted on [[https://github.com][GitHub.com]]:
    - https://github.com/nikadon/supFunSim
    - https://github.com/nikadon/supFunSim/blob/master/supFunSim.org

*** Specific software versions
**** Basic info

     We have tested =supFunSim= using:
     - Lenovo ThinkPad W540:
       - Intel(R) Core(TM) i7-4910MQ CPU @ 2.90GHz,
       - 32768 MB RAM 1600 MHz DDR3,
     - GNU with Linux operating system:
       - SUSE Linux Enterprise Desktop 12 SP2,
       - x86_64-suse-linux-gnu,
       - kernel release 4.4.21-81-default,
     - software:
       - GNU Emacs 25.1.2 (GTK+ Version 3.20.9),
       - MathWorks MATLAB R2017a (9.2.0.538062).

** DOIN [i04] Licensing (GPL)

   Copyright (c) 2014-2017 Jan Nikadon, Tomasz Piotrowski, Dania Gutiérrez

   =supFunSim= contains both ordinary human language text and discrete
   computing machine code (both) provided in this document.  Ordinary
   human language text is copyrighted by authors.  Digital computing
   machine code is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   =supFunSim=  is provided "AS IS", WITHOUT WARRANTY of ANY KIND,
   EXPRESS or IMPLIED, INCLUDING but NOT LIMITED to the WARRANTIES OF
   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE and
   NONINFRINGEMENT. In NO event shall the authors or copyright holders
   be liable for ANY CLAIM, DAMAGES or OTHER LIABILITY, WHETHER in an
   ACTION of CONTRACT, TORT or OTHERWISE, arising from, out of or in
   connection with the software or the use or other dealings in the
   software.

   See the GNU General Public License for more details.  You should
   have received a copy of the GNU General Public License along with
   =supFunSim=. If not, see http://www.gnu.org/licenses/

* DONE Simulations
** DONE [s00] Prelude
*** DONE Description

    The following code block prepares MATLAB workspace, adds paths to
    the necessary toolboxes and [[./fun]] directory containing some additional
    functions.  It also loads initial data:
    - head compartments geometry (cortex)
      - [[./mat/sel_msh.mat]]
    - mesh containing candidates for the location of deep sources (based on thalami)
      - [[./mat/sel_geo_deep_thalami.mat]]
    - mesh containing candidates for the location of deep sources (based on icosahedron642)
      - [[./mat/sel_geo_deep_icosahedron642.mat]]
    - cortex geometry with (anatomical) ROI parcellation
      - [[./mat/sel_atl.mat]]
    - volume conduction model (head-model)
      - [[./mat/sel_vol.mat]]
    - geometry of electrode positions
      - [[./mat/sel_ele.mat]]

    Additionally we also provide:
    - pre-computed leadfields for the above mentioned source candidates
      - [[./mat/sel_src.mat]]

    Contents of the above mentioned files can be checked with:

    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :results silent :tangle no

whos('-file','./mat/sel_atl.mat')
whos('-file','./mat/sel_ele.mat')
whos('-file','./mat/sel_geo_deep_icosahedron642.mat')
whos('-file','./mat/sel_geo_deep_thalami.mat')
whos('-file','./mat/sel_mri00.mat')
whos('-file','./mat/sel_msh.mat')
whos('-file','./mat/sel_src.mat')
whos('-file','./mat/sel_vol.mat')

    #+END_SRC

    Load each file separately (mostly for testing and didactic purpose).

    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :results silent :tangle no

load('./mat/sel_atl.mat')
load('./mat/sel_ele.mat')
load('./mat/sel_geo_deep_icosahedron642.mat')
load('./mat/sel_geo_deep_thalami.mat')
load('./mat/sel_mri00.mat')
load('./mat/sel_msh.mat')
load('./mat/sel_src.mat')
load('./mat/sel_vol.mat')

    #+END_SRC

    NB: For the sake of *backward compatibility* and to facilitate
    *bug tracking*, it also provides initial settings for the
    simulations.  Pleas remember that /ut supra/ those settings (most
    likely) need *not to be modified*.  During batch execution of
    simulations initial settings are *overwritten* by values defined
    in sections [[*%5Bs07%5D%20BATCH][BATCH]] and [[*%5Bs08%5D%20LOOP][LOOP]] (the very last sections of this
    document.)

    NB: In the code we make extensive use of a number of acronyms
    and abbreviations, these include:
    - SrcActiv :: Activity of interest (biological)
    - IntNoise :: Interference (biological) noise
    - BcgNoise :: Background (biological) noise
    - MesNoise :: Measurement noise

    We assume that this file is located in (or symbolically linked to)
    [[~/supFunSim]] on machines running GNU with Linux (or any similar
    decent operating system). Similarly the necessary toolboxes are
    expected to be found in sub-directories of the [[~/toolboxes]]
    directory (or symlinked there). Otherwise the following code might
    need to be adjusted (e.g., by changing value of the =TMP_TBX_PATH=
    variable etc.)  Some additional necessary functions are located in
    the [[./fun]] directory.

    Simulation parameters are controlled using =SETUP= variable.  The
    geometrical arrangement and number of cortical sources in ROIs can
    be controlled using =SRCS= field (=SETUP.SRCS=), which is a 3
    column ~<int>~ array.
    - *Rows* of =SETUP.SRCS= reppresent consequent ROIs.
    - *Cols* of =SETUP.SRCS= represent signal types, namely:
      - ~SrcActiv~,
      - ~IntNoise~ and
      - ~BcgNoise~,
      respectively.
    - Integers represent a number of sources in the given ROI for the
      given signal type.

    Variable =SETUP= also contains field =DEEP= (=SETUP.DEEP=) which
    defines number of signals in the brain center (around thalami) that
    belong to particular signal type (same as for cortical signal types
    these include: ~SrcActiv~, ~IntNoise~ and ~BcgNoise~).

*** !run Simulations main setup

    #+NAME: ./tg_s00_Prelude___01_Simulations_main_setup.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s00_Prelude___01_Simulations_main_setup.m :results silent :var fPath = (file-name-directory (or load-file-name buffer-file-name))

      % Tidy up workspace and change working directory
      clc; close all;clearvars('-except','fPath');
      if exist('fPath'), cd(fPath); else, try, cd('~/supFunSim/'); catch, warningMessage = 'Problem encoutered while trying to change working directory to ''~/supFunSim/''.'; end; end; disp(['CYBERCRAFT:: pwd is: ',pwd])

      % Add path to additional functions
      addpath([pwd,'/fun']);

      % Add path to toolboxes
      if 0, disp('COMMENT: if having problems with toolboxes consider wise use of:'); restoredefaultpath; end
      clearvars TMP_TOOLB_PATH; if exist('~/toolboxes/','dir') == 7, TMP_TOOLB_PATH = '~/toolboxes/'; else, warningMessage = sprintf('CYBERCRAFT:: Warning:: toolboxes directory was not found (use ''~/toolboxes/'')'); end; if exist('TMP_TOOLB_PATH','var'), disp(['CYBERCRAFT:: looking for toolboxes in: ',TMP_TOOLB_PATH]), end;
      addpath([TMP_TOOLB_PATH,'/arfit']);
      addpath([TMP_TOOLB_PATH,'/mvarica']);
      addpath([TMP_TOOLB_PATH,'spm12/']);
      addpath([TMP_TOOLB_PATH,'spm12/toolbox/aal/']);
      addpath([TMP_TOOLB_PATH,'eeglab/']);
      addpath([TMP_TOOLB_PATH,'fieldtrip/']); ft_defaults;
      if exist([TMP_TOOLB_PATH,'fieldtrip/privatePublic/',filesep]),         addpath([TMP_TOOLB_PATH,'fieldtrip/privatePublic/',filesep]);         else, copyfile([TMP_TOOLB_PATH,'fieldtrip/private/',filesep],        [TMP_TOOLB_PATH,'fieldtrip/privatePublic/',filesep]);         addpath([TMP_TOOLB_PATH,'fieldtrip/privatePublic/',filesep]);         end;
      if exist([TMP_TOOLB_PATH,'fieldtrip/forward/privatePublic/',filesep]), addpath([TMP_TOOLB_PATH,'fieldtrip/forward/privatePublic/',filesep]); else, copyfile([TMP_TOOLB_PATH,'fieldtrip/forward/private/',filesep],[TMP_TOOLB_PATH,'fieldtrip/forward/privatePublic/',filesep]); addpath([TMP_TOOLB_PATH,'fieldtrip/forward/privatePublic/',filesep]); end;
      addpath([TMP_TOOLB_PATH,'fieldtrip/utilities']);
      if exist([TMP_TOOLB_PATH,'fieldtrip/utilities/privatePublic/',filesep]), addpath([TMP_TOOLB_PATH,'fieldtrip/utilities/privatePublic/',filesep]); else, copyfile([TMP_TOOLB_PATH,'fieldtrip/utilities/private/',filesep],[TMP_TOOLB_PATH,'fieldtrip/utilities/privatePublic/',filesep]); addpath([TMP_TOOLB_PATH,'fieldtrip/utilities/privatePublic/',filesep]); end;

      % Initialize FieldTrip
      ft_defaults;

      % Load head geometry, electrode positions and ROIs data.
      load('./mat/sel_msh.mat');                     % head compartments geometry (cortex)
      load('./mat/sel_geo_deep_thalami.mat');        % mesh containing candidates for lacation of deep sources (based on thalami)
      load('./mat/sel_geo_deep_icosahedron642.mat'); % mesh containing candidates for lacation of deep sources (based on icosahedron642)
      load('./mat/sel_atl.mat');                     % cortex geometry with (anatomical) ROI parcellation
      load('./mat/sel_vol.mat');                     % volume conduction model (head-model)
      load('./mat/sel_ele.mat');                     % geometry of electrode positions
      load('./mat/sel_src.mat');                     % all cx leadfields

      clearvars ii jj kk mm nn tmp* SETUP

      % Simulations main setup
      SETUP.rROI   = logical(0);       % random (1) or predefined (0) ROIs
      SETUP.rPNT   = logical(0);       % random (1) or predefined (0) candidate points for source locations: if 0, 
				       % number of sources as in SETUP.SRCS(1,1) will be fixed and in close locations
      SETUP.SRCS   = []; % Cortical sources (avoid placing more than 10 sources in single ROI)
      SETUP.SRCS   = [ SETUP.SRCS;  8  0  3 ];
      SETUP.SRCS   = [ SETUP.SRCS;  2  0  0 ];
      SETUP.SRCS   = [ SETUP.SRCS;  2  1  0 ];
      SETUP.SRCS   = [ SETUP.SRCS;  0  4  0 ];
      SETUP.SRCS   = [ SETUP.SRCS;  0  4  0 ];
      SETUP.SRCS   = [ SETUP.SRCS;  0  4  0 ];
      SETUP.SRCS   = [ SETUP.SRCS;  0  4  3 ];
      SETUP.SRCS   = [ SETUP.SRCS;  0  0  3 ];
      SETUP.SRCS   = [ SETUP.SRCS;  0  0  3 ];
      SETUP.SRCS   = [ SETUP.SRCS;  0  0  3 ];
      SETUP.DEEP   = [              2  1  6 ]; % deep sources
      SETUP.ERPs   = 5;       % Add ERPs (timelocked activity)
      SETUP.rROI   = 0;       % random (1) or predefined (0) ROIs
      SETUP.ELEC   = size(sel_ele.elecpos,1); % number of electrodes
      SETUP.n00    = 500;     % number of time samples per trial
      SETUP.K00    = 1;       % number of independent realizations of signal and noise based on generated MVAR model
			      % note: covariance matrix R of observed signal and noise covariance matrix N are estimated from samples originating 
			      % from all realizations of signal and noise 
      SETUP.P00    = 6;       % order of the MVAR model used to generate time-courses for signal of interest
      SETUP.FRAC   = 0.20;    % proportion of ones to zeros in off-diagonal elements of the MVAR coefficients masking array
      SETUP.STAB   = 0.99;    % MVAR stability limit for MVAR eigenvalues (less than 1.0 results in more stable model producing more stationary signals)
      SETUP.RNG    = [0,2.8]; % range for pseudo-random sampling of eigenvalues for MVAR coefficients range
      SETUP.ITER   = 5e5;     % iterations limit for MVAR pseudo-random sampling and stability verification
      SETUP.PDC_RES = [0:0.01:0.5]; % resolution vector for normalized PDC estimation
      SETUP.TELL   = 1;       % provide additional comments during code execution ("tell me more")
      SETUP.PLOT   = 1;       % plot figures during the intermediate stages
      SETUP.SCRN   = get(0,'MonitorPositions'); % get screens positions
      SETUP.DISP   = SETUP.SCRN(end,:);        % force figures to be displayed on (3dr) screenscreen
      SETUP.FIXED_SEED = 0; % Settings for seed selection
      if SETUP.FIXED_SEED, SETUP.SEED = rng(1964);else,SETUP.SEED = rng(round(1e3*randn()^2*sum(clock)));end
      SETUP.RANK_EIG = sum(SETUP.SRCS(:,1)); % rank of EIG-LCMV filter: set to number of active sources
      SETUP.fltREMOVE = 1; % to keep (0) or remove (1) selected filters
      SETUP.SHOWori = 1; % to show (1) or do not show (0) Original and Dummy signals on Figures
      SETUP.IntLfgRANK = round(0.3*sum(SETUP.SRCS(:,2))); % rank of patch-constrained reduced-rank leadfield
      SETUP.supSwitch = 'loc'; % 'rec': run reconstruction of sources activity, 'loc': find active sources

    #+END_SRC

*** !opt Checkups

    NB: The following code is not executed during simulations, it is
    left here only to facilitate *bug tracking* and for
    *explanatory/illustrative purpose*.

    The above applies to all /Checkup/ sections (usually marked with
    =org-todo-keyword=: ~!opt~).

    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos
SETUP
chkSim___tg_s01_PRE___001(SETUP);
sel_atl
sel_ele
sel_ele.label
sel_msh
{sel_msh.bnd.inf}'
sel_vol
SETUP

    #+END_SRC

*** !fun Checkup functions
**** =chkSim___tg_s01_PRE___000_PARSE_SETUP=

     #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/chkSim___tg_s01_PRE___000_PARSE_SETUP.m :results silent

function chkSim___tg_s01_PRE___000_PARSE_SETUP(SETUP,sel_atl)
% Simulations setup parser

    disp('Checking consistency of ''SETUP'' definitions')

    % check if number of ROIs is smaller than number of cortex parcels
    tmp_MaxROIs = size(sel_atl.Atlas(sel_atl.atl).Scouts,2);
    tmp_ActROIs = size(SETUP.SRCS,1);
    if 0
        disp(SETUP.SRCS)
        tmp_MaxROIs
        tmp_ActROIs
    end
    if tmp_ActROIs > tmp_MaxROIs,
        tmp_msg = ['please avoid putting more than ', num2str(tmp_MaxROIs) ' source ROIs (decrease the number of rows in the ''SETUP.SRCS''.)'];
        error(tmp_msg);
    end

    % check if no roi contains more sources than allowed (allowed number
    % of sources is the number of nodes in the ROI that contains the
    % least nodes)
    tmp_MaxSrc = []; for ii = 1:tmp_MaxROIs, tmp_MaxSrc = [tmp_MaxSrc; size(sel_atl.Atlas(sel_atl.atl).Scouts(ii).Vertices,2)]; end;
    tmp_MaxSrc = min(tmp_MaxSrc);
    [tmp_ActSrcs, tmp_ActColMaxIdx] = max(sum(SETUP.SRCS,2));
    if 0
        SETUP.SRCS
        tmp_MaxSrc
        tmp_ActSrcs
        tmp_ActColMaxIdx
    end
    if tmp_ActSrcs > min(tmp_MaxSrc),
        tmp_msg = ['please avoid putting more than ', num2str(tmp_MaxSrc) ' sources in any of the ROIs (decrease sum of the column number ', num2str(tmp_ActColMaxIdx), ' in the ''SETUP.SRCS''.)'];
        error(tmp_msg);
    end

    % Check if a number of ERPs is not greater than the number of signals
    % of interest.
    tmp_MaxERPs = sum(SETUP.SRCS(:,1));
    tmp_ActERPs = SETUP.ERPs;
    if 0
        SETUP.SRCS
        tmp_MaxERPs
        tmp_ActERPs
    end
    if tmp_ActERPs > tmp_MaxERPs,
        tmp_msg = ['please avoid putting more than ', num2str(tmp_MaxERPs) ' ERPs (this is a number of sources of activity of interest, please change the sum in the 1st column of ''SETUP.SRCS'' or the ''SETUP.ERPs'' value.)'];
        error(tmp_msg);
    end

end

     #+END_SRC

**** =chkSim___tg_s01_PRE___001=

     #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/chkSim___tg_s01_PRE___001.m :results silent

function chkSim___tg_s01_PRE___001(SETUP)
    fprintf('     \n');
    fprintf('CYBERCRAFT:: Number of signal types per CORTEX ROI:\n\n');
    disp(array2table([[1:size(SETUP.SRCS,1)]',SETUP.SRCS,sum(SETUP.SRCS,2)],'VariableNames',{'ROI','SrcActiv','IntNoise','BcgNoise','TOTAL'}))
    fprintf('     \n');
    fprintf('CYBERCRAFT:: TOTAL number of signals for CORTEX:\n\n');
    disp(array2table([size(SETUP.SRCS,1);sum(SETUP.SRCS(:,1));sum(SETUP.SRCS(:,2));sum(SETUP.SRCS(:,3));sum(sum(SETUP.SRCS(:,1:3)))]','VariableNames',{'ROI','SrcActiv','IntNoise','BcgNoise','TOTAL'}));
    fprintf('     \n');
    fprintf('CYBERCRAFT:: Number of signals for DEEP sources:\n\n');
    disp(array2table([size(SETUP.SRCS,1)+1;sum(SETUP.DEEP(1,1));sum(SETUP.DEEP(1,2));sum(SETUP.DEEP(1,3));sum(SETUP.DEEP(1,1:3))]','VariableNames',{'ROI','SrcActiv','IntNoise','BcgNoise','TOTAL'}));
    fprintf('     \n');
    fprintf('CYBERCRAFT:: TOTAL number of signals for CORTEX and DEEP sources:\n\n');
    disp(array2table([NaN;sum(SETUP.SRCS(:,1))+sum(SETUP.DEEP(1,1));sum(SETUP.SRCS(:,2))+sum(SETUP.DEEP(1,2));sum(SETUP.SRCS(:,3))+sum(SETUP.DEEP(1,3));sum(sum(SETUP.SRCS(:,1:3)))+sum(SETUP.DEEP(1,1:3))]','VariableNames',{'ROI','SrcActiv','IntNoise','BcgNoise','TOTAL'}));
    fprintf('     \n');
    disp(['CYBERCRAFT:: MesNoise: ',num2str(SETUP.ELEC)])
end

     #+END_SRC

*** !run SNR adjustment and signal components setup

    This section provides further simulation settings.

    NB: For the sake of *backward compatibility* and to facilitate *bug
    tracking*, it also provides initial settings for the simulations.
    Pleas remember that /ut supra/ those settings (most likely) need
    *not to be modified* because for execution of simulations in batch
    they will be *overwritten* by values defined in sections [[*%5Bs07%5D%20BATCH][BATCH]] and
    [[*%5Bs08%5D%20LOOP][LOOP]], the two sections that can be found very close to the end of
    this document.

    #+NAME: ./tg_s00_Prelude___02_SNR_Adjustments.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s00_Prelude___02_SNR_Adjustments.m :results silent

SETUP.CUBE           = 20;    % perturbation of the leadfields based on the shift of source position within a cube of given edge length (centered at the original leadfields positions), in [mm]
SETUP.CONE           = pi/32; % perturbation of the leadfields based on the rotation of source orientation (azimuth TH, elevation PHI), in [rad]
SETUP.H_Src_pert     = 0;     % use original (0) or perturbed (1) leadfield for signal reconstruction
SETUP.H_Int_pert     = 0;     % use original (0) or perturbed (1) leadfield for nulling constrains
SETUP.SINR           = 5;     % signal to interference noise power ratio expressed in dB (both measured on electrode level)
SETUP.SBNR           = 10;    % signal to biological noise power ratio expressed in dB (both measured on electrode level)
SETUP.SMNR           = 15;    % signal to measurment noise power ratio expressed in dB (both measured on electrode level)
SETUP.WhtNoiseAddFlg = 1;     % white noise admixture in biological noise interference noise (FLAG)
SETUP.WhtNoiseAddSNR = 3;     % SNR of BcgNoise and WhiNo (dB)
SETUP.SigPre = 0;   SETUP.IntPre = 0;   SETUP.BcgPre = 1;   SETUP.MesPre = 1; % final signal components for pre-interval  (use zero or one for signal, interference noise, biological noise, measurement noise)
SETUP.SigPst = 1;   SETUP.IntPst = 1;   SETUP.BcgPst = 1;   SETUP.MesPst = 1; % final signal components for post-interval (as above)

     % For localization, the default is to consider random locations of ROI with some fixed candidate points such that 
     % SETUP.SRCS(1,1) of them are in close locations; additionally, no interfering sources are considered
     % in the localization model. You may comment out the 'if' section below to experiment with other settings.
     if(SETUP.supSwitch == 'loc') 
	 SETUP.rROI   = logical(1);       % random (1) or predefined (0) ROIs
	 SETUP.rPNT   = logical(0);       % random (1) or predefined (0) candidate points for source locations: if 0, 
					  % number of sources as in SETUP.SRCS(1,1) will be fixed and in close locations
	 SETUP.SigPre = 0;   SETUP.IntPre = 0;   SETUP.BcgPre = 1;   SETUP.MesPre = 1; % final signal components for pre-interval  (use zero or one for signal, interference noise, biological noise, measurement noise)
	 SETUP.SigPst = 1;   SETUP.IntPst = 0;   SETUP.BcgPst = 1;   SETUP.MesPst = 1; % final signal components for post-interval (as above)
     end

if SETUP.rPNT
    disp('CC: using random source locations')
else
    disp('CC: using predefined source locations')
end

    #+END_SRC

*** !opt Checkups

    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos
SETUP
chkSim___tg_s01_PRE___001(SETUP)
chkSim___tg_s02_PRE___001(SETUP)

    #+END_SRC

*** !fun Checkup functions
**** =chkSim___tg_s02_PRE___001=

     #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/chkSim___tg_s02_PRE___001.m :results silent

function chkSim___tg_s02_PRE___001(SETUP)
    fprintf('\n');
    fprintf('CYBERCRAFT:: Perturbation cube:\n\n');
    disp(array2table(SETUP.CUBE,'VariableNames',{'mm'},'RowNames',{'shift distance'}))
    fprintf('\n');
    fprintf('CYBERCRAFT:: Perturbation cone:\n\n');
    disp(array2table([SETUP.CONE,rawRad2Deg(SETUP.CONE)],'VariableNames',{'rad','deg'},'RowNames',{'rotation angle'}));
    fprintf('\n');
    fprintf('CYBERCRAFT:: Ratios for signal to:\n\n');
    disp(array2table([SETUP.SINR;SETUP.SBNR;SETUP.SMNR],'VariableNames',{'SxNR'},'RowNames',{'IntNoise','BcgNoise','MesNoise'}));
    fprintf('\n');
    fprintf('CYBERCRAFT:: Signal components:\n\n');
    disp(array2table([SETUP.SigPre,SETUP.IntPre,SETUP.BcgPre,SETUP.MesPre;SETUP.SigPst,SETUP.IntPst,SETUP.BcgPst,SETUP.MesPst]','VariableNames',{'pre','post'},'RowNames',{'SrcActiv','IntNoise','BcgNoise','MesNoise'}));
    fprintf('\n');
end

     #+END_SRC

** DONE [s01] Time series for activity and noise sources
*** DONE Description

    Code in this section produces timeseries for the signals described
    in the table below.

    | Variable name    | Description                                         |
    |------------------+-----------------------------------------------------|
    | sim_sig_SrcActiv | Activity of interest (biological)                   |
    | sim_sig_IntNoise | Interference (biological) noise                     |
    | sim_sig_BcgNoise | Background (biological) noise                       |
    |------------------+-----------------------------------------------------|
    | sim_sig_MesNoise | Measurement noise                                   |
    |------------------+-----------------------------------------------------|
    | sim_sig_AdjSNRs  | All the above signals combined and adjusted for SNR |

*** MVAR modeling for bioelectrical activity sources
**** Info

     In this section multivariate timeseries are generated for:
     =sim_sig_SrcActiv=, =sim_sig_IntNoise=, =sim_sig_BcgNoise=.

     Timeseries are based on Multivariate Autoregressive (MVAR) model.

**** MVAR modeling (theory)
***** MVAR model basics

      See:
      - [[./toolboxes/arfit/arfit.pdf]]
      - [[./toolboxes/arfit/arfit_alg.pdf]]
      - [[./toolboxes/mvarica/arfit/arfit.pdf]]
      - [[./toolboxes/mvarica/arfit/arfit_alg.pdf]]
      - [[~/cc_overkill/pdf/papers/2001 - Neumaier,Schneider - Algorithm 808: ARFIT —A Matlab Package for the Estimation of Parameters and Eigenmodes of Multivariate Autoregressive Models.pdf][2001 - Neumaier,Schneider - Algorithm 808: ARFIT —A Matlab Package for the Estimation of Parameters and Eigenmodes of Multivariate Autoregressive Models.pdf]]
      - https://sccn.ucsd.edu/wiki/Chapter_3._Multivariate_Autoregressive_Modeling

      # {{x}_{t}}=v+\underset{k=1}{\overset{p}{\mathop \sum }}\,{{A}_{k}}{{x}_{t-k}}+{{u}_{t}}

      ${{v}_{\nu}}=w+\underset{l=1}{\overset{p}{\mathop \sum }}\,{{A}_{l}}{{v}_{\nu-l}}+{{\varepsilon}_{\nu}}$, where

      - ${{v}_{\nu}}$ are the $m$-dimensional state vectors of (stationary)
	time series,

      - $p$ is the order of the model,

      - matrices ${{A}_{1}}\ldots{{A}_{p}} \in R^{m\times{}m}$ are the
	coefficient matrices of the AR model,

      - ${{\varepsilon}_{\nu}}$ are the $m$-dimensional uncorrelated
	random vectors with mean zero and covariance matrix $C \in
	R^{m\times{}m}$

***** PDC basics

      - [[~/cc_overkill/pdf/papers/2001 - Baccala, Sameshima, Partial directed coherence: a new concept in neural structure determination.pdf][2001 - Baccala, Sameshima, Partial directed coherence: a new concept in neural structure determination.pdf]]
      - [[~/cc_overkill/pdf/books/2014 - SAMESHIMA,BACCALA - Methods in Brain Connectivity Inference through Multivariate Time Series Analysis.pdf][2014 - SAMESHIMA,BACCALA - Methods in Brain Connectivity Inference through Multivariate Time Series Analysis.pdf]]

****** Normalized frequency

       See
       - [[https://en.wikipedia.org/wiki/Normalized_frequency_(unit)][Wiki: EN: Normalized frequency (unit)]]
       - [[https://en.wikipedia.org/wiki/Normalized_frequency_%2528unit%2529#Alternative_normalizations][Wiki: EN: Normalized frequency (unit): Alternative normalizations]]
         - Some programs (such as MATLAB) that design filters with
           real-valued coefficients use the Nyquist frequency
           ($\textstyle f_s/2$) as the [[https://en.wikipedia.org/wiki/Normalizing_constant][normalization constant]].  The
           resultant normalized frequency has units of
           /half-cycles/sample/ or equivalently /cycles per 2 samples/.
         - Sometimes, the unnormalized frequency is represented in units
           of [[https://en.wikipedia.org/wiki/Radian_per_second][radians/second]] ([[https://en.wikipedia.org/wiki/Angular_frequency][angular frequency]]),
           and denoted by $\textstyle \omega$. When
           $\textstyle \omega$ is normalized by the
           sample-rate (/samples/sec/), the resulting units are
           /radians/sample/.  The normalized Nyquist frequency is
           π /radians/sample/, and the normalized sample-rate is
           2π /radians/sample/.
         - The following table shows examples of normalized frequencies
           for a 1 kHz signal, a sample rate $\textstyle
           f_\mathrm{s}$ = [[https://en.wikipedia.org/wiki/44,100_Hz][44.1 kHz]], and 3 different choices of
           normalized units.  Also shown is the frequency region
           containing one cycle of the [[https://en.wikipedia.org/wiki/Discrete-time_Fourier_transform][discrete-time Fourier transform]],
           which is always a periodic function.

           | Units              | Domain           | *Remarks* | Computation          |   Value |
           |--------------------+------------------+-----------+----------------------+---------|
           | cycles/sample      | [-½, ½] or [0,1] |           | 1000 / 44100         | 0.02268 |
           | half-cycles/sample | [-1,1] or [0,2]  | *MATLAB*  | 1000 / 22050         | 0.04535 |
           | radians/sample     | [-π,π] or [0,2π] |           | 2 ''π'' 1000 / 44100 |  0.1425 |

       - http://www.mathworks.com/help/signal/ug/frequency-response.html#zmw57dd0e1407
       - http://www.mathworks.com/help/signal/ref/freqz.html
       - https://www.mathworks.com/matlabcentral/newsreader/view_thread/286192
       - https://www.mathworks.com/matlabcentral/answers/50575-normalized-frequency-in-analog-filter-design
       - https://www.codementor.io/tips/8102473731/understanding-normalized-frequency-in-matlab
       - http://dsp.stackexchange.com/questions/17490/how-to-normalize-frequency-in-matlab
       - http://dsp.stackexchange.com/questions/15219/convert-normalized-frequency-to-real-frequency-in-ar-model

**** Generation of timeseries for bioelectrical activity of interest
***** !run Code

      #+NAME: ./tg_s01_Timeseries___01_SrcActiv.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s01_Timeseries___01_SrcActiv.m :results output silent

clearvars sim_sig_SrcActiv tmp*
sim_sig_SrcActiv = cccSim___makeSimSig(SETUP,1,'MVAR based signal for sources activity');
clearvars ii jj kk nn tmp*

if SETUP.ERPs > 0

    sim_sig_SrcActiv.sigSRC_pst_orig = sim_sig_SrcActiv.sigSRC_pst;
    sim_sig_SrcActiv.sigSRC_pre_orig = sim_sig_SrcActiv.sigSRC_pre;

    tmp_LB = -5;
    tmp_UB = 7;
    tmp_NN = SETUP.n00;
    tmp_PP = 1;
    for ee = 1:1:SETUP.ERPs
        tmp_PP = mod(ee,3)+1;
        sim_sig_SrcActiv.ERPs(ee,:) = gauswavf(tmp_LB,tmp_UB,tmp_NN,tmp_PP);
    end
    sim_sig_SrcActiv.ERPs_pst = transpose(sim_sig_SrcActiv.ERPs);
    sim_sig_SrcActiv.ERPs_pst(size(sim_sig_SrcActiv.sigSRC_pst,1),size(sim_sig_SrcActiv.sigSRC_pst,2)) = 0;
    sim_sig_SrcActiv.ERPs_pst = cat(3,repmat(sim_sig_SrcActiv.ERPs_pst,[1,1,SETUP.K00]));

    sim_sig_SrcActiv.sigSRC_pst = sim_sig_SrcActiv.sigSRC_pst + 20*sim_sig_SrcActiv.ERPs_pst;

end

      #+END_SRC

      NB.: Fields:
      - =w01=,
      - =A01=,
      - =C01=,
      - =SBC01=,
      - =FPE01= and
      - =th01=
      (in general =*01=) are estimated from the =sigMVAR_pst= signal.

***** !fun MVAR timeseries modeling functions
****** =cccSim___makeSimSig=

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./fun/cccSim___makeSimSig.m :results silent

function sim_sig_SrcActiv = cccSim___makeSimSig(SETUP,set_COL,varargin)
    sim_sig_SrcActiv.inf        = '';
    sim_sig_SrcActiv.S00        = [];
    sim_sig_SrcActiv.P00        = [];
    sim_sig_SrcActiv.M00        = [];
    sim_sig_SrcActiv.A00        = zeros([0,1]);
    sim_sig_SrcActiv.mdlStabH   = [];
    sim_sig_SrcActiv.mdlStabM   = [];
    sim_sig_SrcActiv.w00        = [];
    sim_sig_SrcActiv.C00        = [];
    sim_sig_SrcActiv.n00        = [];
    sim_sig_SrcActiv.K00        = [];
    sim_sig_SrcActiv.sigSRC_pre = permute([],[3,2,1]);
    sim_sig_SrcActiv.sigSRC_pst = permute([],[3,2,1]);
    sim_sig_SrcActiv.w01        = [];
    sim_sig_SrcActiv.A01        = zeros([0,1]);
    sim_sig_SrcActiv.C01        = [];
    sim_sig_SrcActiv.SBC01      = [];
    sim_sig_SrcActiv.FPE01      = [];
    sim_sig_SrcActiv.th01       = [];
    if ~isempty(varargin)
        sim_sig_SrcActiv.inf = varargin{end};
    end
    sim_sig_SrcActiv.S00 = sum(SETUP.SRCS(:,set_COL))+SETUP.DEEP(1,set_COL); % number of signals
    if sim_sig_SrcActiv.S00 > 0
        sim_sig_SrcActiv.P00 = SETUP.P00;                     % order of the MVAR model used to generate time-courses
        sim_sig_SrcActiv.M00 = cccSim___diagonMask(sim_sig_SrcActiv.S00,SETUP.FRAC); % MVAR coefficients mask
        sim_sig_SrcActiv.A00 = cccSim___stableMVAR(sim_sig_SrcActiv.S00,sim_sig_SrcActiv.P00,sim_sig_SrcActiv.M00,SETUP.RNG,SETUP.STAB,SETUP.ITER);
        [sim_sig_SrcActiv.mdlStabH,sim_sig_SrcActiv.mdlStabM] = rawIsStableMVAR(sim_sig_SrcActiv.A00,1);
        sim_sig_SrcActiv.w00 = zeros(sim_sig_SrcActiv.S00,1); % expected value for time-courses
        sim_sig_SrcActiv.C00 = eye(sim_sig_SrcActiv.S00);     % covariance of MVAR model's "driving noise"
        sim_sig_SrcActiv.n00 = SETUP.n00;                     % number of time samples
        sim_sig_SrcActiv.K00 = SETUP.K00;                     % number of independent realizations of the driving AR models
        sim_sig_SrcActiv.sigSRC_pre = 10*arsim(sim_sig_SrcActiv.w00,sim_sig_SrcActiv.A00,sim_sig_SrcActiv.C00,[sim_sig_SrcActiv.n00,2*sim_sig_SrcActiv.K00],1e3); % 10e-6; % Current density (10nA*mm)
        sim_sig_SrcActiv.sigSRC_pst = sim_sig_SrcActiv.sigSRC_pre(:,:,[end/2+1:end]);
        sim_sig_SrcActiv.sigSRC_pre = sim_sig_SrcActiv.sigSRC_pre(:,:,[1:end/2]);
        if SETUP.TELL,disp('Fitting MVAR model to generated signal...');end;
        [sim_sig_SrcActiv.w01,sim_sig_SrcActiv.A01,sim_sig_SrcActiv.C01,sim_sig_SrcActiv.SBC01,sim_sig_SrcActiv.FPE01,sim_sig_SrcActiv.th01] = arfit(sim_sig_SrcActiv.sigSRC_pst,sim_sig_SrcActiv.P00-round(sim_sig_SrcActiv.P00/2),sim_sig_SrcActiv.P00+round(sim_sig_SrcActiv.P00/2));
    end
end

       #+END_SRC

****** =cccSim___stableMVAR=

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./fun/cccSim___stableMVAR.m :results silent

function A00 = cccSim___stableMVAR(S00,P00,M00,set_RNG,set_STAB,set_ITER)
% procedure inspired by function stablemvar that is a part of MVARICA Toolbox.
    if 0
        A00 = cccSim___stableMVAR(S00,P00,M00,set_RNG,set_ITER,set_STAB)
    end
    tmp_iterNow = 0;
    tmp_lambda = Inf;
    while any(abs(tmp_lambda)>set_STAB) && tmp_iterNow < set_ITER
        tmp_V = orth(rand(S00*P00,S00*P00));
        tmp_U = orth(rand(S00*P00,S00*P00));
        lambdatmp = set_RNG(1)+(set_RNG(2)-set_RNG(1))*rand(S00*P00,1);
        tmp_A00 = tmp_V*diag(lambdatmp)*tmp_U';
        A00 = tmp_A00(1:S00,:);
        A00 = A00.*repmat(M00,1,P00); % nulling of some coefficients based on mask
        tmp_lambda = eig([A00; eye((P00-1)*S00) zeros((P00-1)*S00,S00)]);
        tmp_iterNow = tmp_iterNow + 1;
    end
    if tmp_iterNow >=set_ITER
        A0=[];
        error('Could not generate stable MVAR model in given number of iterations (set_ITER)');
    end
    if 0
        S00 = sim_sig00.S00
        P00 = sim_sig00.P00
        w00 = sim_sig00.w00
        C00 = sim_sig00.C00
        n00 = sim_sig00.n00
        K00 = sim_sig00.K00
    end
end

       #+END_SRC

****** =cccSim___diagonMask=

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./fun/cccSim___diagonMask.m :results silent

function M00 = cccSim___diagonMask(S00,M00_frc)
% cccSim___diagonMask produces masking square array M (SxS) whose all diagonal and some off-diagonal elements are ones. All
% othere elements are zeros. The off-diagonal elements are sampled pseudo-randomly and their number is declared by
% second argument of this function which should be a number between zero and one that describes the proportion of ones
% to zeros among all off-diagonal elements of the resulting array.
%
% Usage:
%
%    MSK = cccSim___diagonMask(S,frac)

% Copyright (C) 2000-2015, Jan (CyberCraft) Nikadon (nikadon-AT-SIGN-gmail.com)
%
%    This file is free software and a part of CyberCraft Toolkit.  You can redistribute it and/or modify it under the
%    terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the
%    License, or (at your option) any later version.
%
%    CyberCraft is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with CyberCraft. If not, see <http://www.gnu.org/licenses/>.

    if 0
        % Example useage:
        S00     = 10;
        M00_frc = 0.1;
        M00     = cccSim___diagonMask(S00,M00_frc)
        figure(1);clf;
        imagesc(M00);colorbar;
        (numel(find(M00))-S00) / (numel(M00)-S00)
        imagesc(sim_sig00.M00);colorbar;
    end
    M00     = eye(S00);
    M00_idx = find(~M00);
    M00_smp = datasample(M00_idx,round(M00_frc*numel(M00_idx)),'Replace',false);
    M00(M00_smp) = deal(1);
end

       #+END_SRC

***** !opt Checkups
****** !opt Gauss wave

       #+BEGIN_SRC emacs-lisp :eval no :tangle no

size(sim_sig_SrcActiv.sigSRC_pst_orig)
size(sim_sig_SrcActiv.sigSRC_pst)
size(sim_sig_SrcActiv.ERPs_pst)
size(sim_sig_SrcActiv.ERPs)

figure, imagesc(sim_sig_SrcActiv.sigSRC_pst_orig(:,:,1))
figure, imagesc(sim_sig_SrcActiv.sigSRC_pst(:,:,1))
figure, imagesc(sim_sig_SrcActiv.ERPs_pst(:,:,1))

ee = 0
ee = ee + 1

tmp_LB = -5;
tmp_UB = 7;
tmp_NN = SETUP.n00;
tmp_PP = mod(ee,3)+1;

tmp_PP

[PSI,X] = gauswavf(tmp_LB,tmp_UB,tmp_NN,tmp_PP);

close all
plot(PSI)

size(cat(2,20*repmat(transpose(sim_sig_SrcActiv.ERPs),[1,1,SETUP.K00]),zeros(SETUP.n00,sum(SETUP.SRCS(:,1))-SETUP.ERPs,SETUP.K00)))

close all
figure
plot(sim_sig_SrcActiv.ERPs(3,:))

figure
plot(sim_sig_SrcActiv.sigSRC_pst(:,3,1))

figure
plot(mean(sim_sig_SrcActiv.sigSRC_pst(:,3,:),3))


       #+END_SRC

****** !opt Basics

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos
chkSim___tg_s01_PRE___001(SETUP)
sim_sig_SrcActiv

tmp_arrC = [rawSize(sim_sig_SrcActiv.sigSRC_pre,[1,2,3]); rawSize(sim_sig_SrcActiv.sigSRC_pst,[1,2,3])];
tmp_varN = {'Timesamples','Signals','Realizations'};
tmp_rowN = {'sim_sig_SrcActiv.sigSRC_pre','sim_sig_SrcActiv.sigSRC_pst'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN))

tmp_arrC = [size(sim_sig_SrcActiv.A00),size(sim_sig_SrcActiv.A00,2)/size(sim_sig_SrcActiv.A00,1);size(sim_sig_SrcActiv.A01),size(sim_sig_SrcActiv.A01,2)/size(sim_sig_SrcActiv.A01,1)];
tmp_varN = {'Signals','Signals_cdot_ModOrd','ModOrd'};
tmp_rowN = {'sim_sig_SrcActiv.A00','sim_sig_SrcActiv.A01'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN))

~isequal(sim_sig_SrcActiv.sigSRC_pre,sim_sig_SrcActiv.sigSRC_pst)
~isequal(sim_sig_SrcActiv.A00,sim_sig_SrcActiv.A01)

       #+END_SRC

****** !fig MVAR zeroing matrix (=M00=)

       Matrix produced by =cccSim___diagonMask= (mask for nulling of
       some MVAR coefficients).

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

figure(101);clf;set(gcf,'Position',SETUP.DISP);imagesc(sim_sig_SrcActiv.M00);colormap(rawHotColdColorMap(10));colorbar;
set(gcf,'color','w');

       #+END_SRC

****** !fig MVAR model coefficient matrix (=A00=)

       Plotted are:
       - =A00= the original coefficient matrix used for timeseries generation and
       - =A01= coefficient matrix obtained by fitting to the generated timeseries.

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

figure(105);clf;set(gcf,'Position',SETUP.DISP);subplot(2,1,1);cccSim___rawDispA00(sim_sig_SrcActiv.A00);title('A00');subplot(2,1,2);cccSim___rawDispA00(sim_sig_SrcActiv.A01);title('A01');
set(gcf,'color','w');

       #+END_SRC

****** !fig PDC

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

figure(120);clf;set(gcf,'Position',SETUP.DISP);cccSim___rawPlotA00(sim_sig_SrcActiv.A00);
set(gcf,'color','w');

figure(121);clf;set(gcf,'Position',SETUP.DISP);cccSim___rawPlotA00(sim_sig_SrcActiv.A01);
set(gcf,'color','w');

       #+END_SRC

***** !fun Checkup functions
****** =cccSim___rawDispA00=

       #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/cccSim___rawDispA00.m :results silent

function cccSim___rawDispA00(A00)
    rawImgSC(A00,8);
    tmp_base = 255;
    tmp_colorMapMat = rawHotColdColorMap(tmp_base);
    caxis([-max(abs(min(min(A00))),abs(max(max(A00)))) max(abs(min(min(A00))),abs(max(max(A00))))]);
    colormap(tmp_colorMapMat);
    colorbar;
    hold on;
    tmp_stem = 0.5+size(A00,1):size(A00,1):size(A00,2)-0.5;
    tmp_vals = size(A00,1)+0.5*ones(size(tmp_stem));
    stem(tmp_stem,tmp_vals,'Color','k','LineWidth',0.5,'Marker', 'none');
    set(gca,'XTick',[0:size(A00,1):size(A00,2)])
end

       #+END_SRC

****** =cccSim___rawPlotA00=

       #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/cccSim___rawPlotA00.m :results silent

function cccSim___rawPlotA00(A00)
    rawPlotPDC(abs(PDC(A00,SETUP.PDC_RES)),SETUP.PDC_RES);
end

       #+END_SRC

**** Generation of timeseries for bioelectrical interference noise

     In general this signal is intended to be highly correlated (or
     anti-correlated) with the signal of interest. Also some white
     noise can be added.

***** !run Code

      #+NAME: ./tg_s01_Timeseries___02_IntNoise.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s01_Timeseries___02_IntNoise.m :results output silent

clearvars sim_sig_IntNoise tmp*
sim_sig_IntNoise = sim_sig_SrcActiv;

sim_sig_IntNoise.inf        = 'MVAR based signal for interference (biological) noise activity sources';
tmp_IntNoiseCol = 2;

if SETUP.TELL,disp('Getting base for IntNoise...');end;
sim_sig_IntNoise.sigSRC_pre = -sim_sig_IntNoise.sigSRC_pre;
sim_sig_IntNoise.sigSRC_pst = -sim_sig_IntNoise.sigSRC_pst;
[sim_sig_IntNoise.w01,sim_sig_IntNoise.A01,sim_sig_IntNoise.C01,sim_sig_IntNoise.SBC01,sim_sig_IntNoise.FPE01,sim_sig_IntNoise.th01] = deal([]);

if SETUP.TELL,disp('Populating signals for IntNoise...');end;
sim_sig_IntNoise.sigSRC_pre = repmat(sim_sig_IntNoise.sigSRC_pre,1,ceil( (sum(SETUP.SRCS(:,tmp_IntNoiseCol))+SETUP.DEEP(1,tmp_IntNoiseCol))/size(sim_sig_IntNoise.sigSRC_pre,2)),1);
sim_sig_IntNoise.sigSRC_pst = repmat(sim_sig_IntNoise.sigSRC_pst,1,ceil( (sum(SETUP.SRCS(:,tmp_IntNoiseCol))+SETUP.DEEP(1,tmp_IntNoiseCol))/size(sim_sig_IntNoise.sigSRC_pst,2)),1);

if SETUP.TELL,disp('Reducing dimensionality of IntNoise...');end;
sim_sig_IntNoise.sigSRC_pre = sim_sig_IntNoise.sigSRC_pre(:,1:(sum(SETUP.SRCS(:,tmp_IntNoiseCol))+SETUP.DEEP(1,tmp_IntNoiseCol)),:);
sim_sig_IntNoise.sigSRC_pst = sim_sig_IntNoise.sigSRC_pst(:,1:(sum(SETUP.SRCS(:,tmp_IntNoiseCol))+SETUP.DEEP(1,tmp_IntNoiseCol)),:);

if SETUP.TELL,disp('Backing-up IntNoise before white noise admixture...');end;
sim_sig_IntNoise.sigSRC_pre_b4admix = sim_sig_IntNoise.sigSRC_pre;
sim_sig_IntNoise.sigSRC_pst_b4admix = sim_sig_IntNoise.sigSRC_pst;

if SETUP.WhtNoiseAddFlg
    if SETUP.TELL,disp('Adding some white noise to the biological noise activity...');end;
    sim_sig_IntNoise.noiseAdmix_pre = randn(size(sim_sig_IntNoise.sigSRC_pre));
    sim_sig_IntNoise.noiseAdmix_pst = randn(size(sim_sig_IntNoise.sigSRC_pst));
    for kk = 1:SETUP.K00
        sim_sig_IntNoise.noiseAdmix_pre(:,:,kk) = rawAdjTotSNRdB(sim_sig_IntNoise.sigSRC_pre(:,:,kk),sim_sig_IntNoise.noiseAdmix_pre(:,:,kk),SETUP.WhtNoiseAddSNR);
        sim_sig_IntNoise.noiseAdmix_pst(:,:,kk) = rawAdjTotSNRdB(sim_sig_IntNoise.sigSRC_pst(:,:,kk),sim_sig_IntNoise.noiseAdmix_pst(:,:,kk),SETUP.WhtNoiseAddSNR);
    end
    sim_sig_IntNoise.sigSRC_pre = sim_sig_IntNoise.sigSRC_pre+sim_sig_IntNoise.noiseAdmix_pre;
    sim_sig_IntNoise.sigSRC_pst = sim_sig_IntNoise.sigSRC_pst+sim_sig_IntNoise.noiseAdmix_pst;
end

if SETUP.TELL,disp('Fitting MVAR model to generated signal...');end;
[sim_sig_IntNoise.w01,sim_sig_IntNoise.A01,sim_sig_IntNoise.C01,sim_sig_IntNoise.SBC01,sim_sig_IntNoise.FPE01,sim_sig_IntNoise.th01] = arfit(sim_sig_IntNoise.sigSRC_pst,sim_sig_IntNoise.P00-round(sim_sig_IntNoise.P00/2),sim_sig_IntNoise.P00+round(sim_sig_IntNoise.P00/2));

if 0
    sim_sig_IntNoise = rmfield(sim_sig_IntNoise,{'sigSRC_pre_b4admix','sigSRC_pst_b4admix','noiseAdmix_pre','noiseAdmix_pst'})
end

clearvars ii jj kk nn tmp*

      #+END_SRC

***** !opt Checkups
****** !opt Basics 1

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos
sim_sig_SrcActiv
sim_sig_IntNoise

chkSim___tg_s01_PRE___001(SETUP)

kk = 1;
jj = 1:min(size(sim_sig_SrcActiv.sigSRC_pre,2),size(sim_sig_IntNoise.sigSRC_pre,2));

tmp_arrC = [rawSNRdB(sim_sig_IntNoise.sigSRC_pre_b4admix(:,:,kk),sim_sig_IntNoise.noiseAdmix_pre(:,:,kk));rawSNRdB(sim_sig_IntNoise.sigSRC_pst_b4admix(:,:,kk),sim_sig_IntNoise.noiseAdmix_pst(:,:,kk))];
tmp_varN = {'SNR'};
tmp_rowN = {'sim_sig_IntNoise.sigSRC_pre_b4admix','sim_sig_IntNoise.sigSRC_pst_b4admix'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN))

tmp_arrC = corrcoef(sim_sig_SrcActiv.sigSRC_pre(:,jj,kk),sim_sig_IntNoise.sigSRC_pre(:,jj,kk));
tmp_varN = {'sim_sig_SrcActiv_a','sim_sig_IntNoise_a'};
tmp_rowN = tmp_varN;
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN))

tmp_arrC = corrcoef(sim_sig_SrcActiv.sigSRC_pst(:,jj,kk),sim_sig_IntNoise.sigSRC_pst(:,jj,kk));
tmp_varN = {'sim_sig_SrcActiv_b','sim_sig_IntNoise_b'};
tmp_rowN = tmp_varN;
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN))

tmp_arrC = corrcoef(sim_sig_SrcActiv.sigSRC_pre(:,jj,kk),sim_sig_IntNoise.sigSRC_pst(:,jj,kk));
tmp_varN = {'sim_sig_SrcActiv_a','sim_sig_IntNoise_b'};
tmp_rowN = tmp_varN;
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));disp('CYBERCRAFT:: This should be negligible')

tmp_arrC = corrcoef(sim_sig_SrcActiv.sigSRC_pst(:,jj,kk),sim_sig_IntNoise.sigSRC_pre(:,jj,kk));
tmp_varN = {'sim_sig_SrcActiv_b','sim_sig_IntNoise_a'};
tmp_rowN = tmp_varN;
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));disp('CYBERCRAFT:: This should be negligible')

clearvars ii jj kk nn tmp*

       #+END_SRC

****** !opt Basics 2

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos
sim_sig_SrcActiv
sim_sig_IntNoise

disp(~isequal(sim_sig_SrcActiv.sigSRC_pre,sim_sig_SrcActiv.sigSRC_pst))
disp(~isequal(sim_sig_SrcActiv.sigSRC_pre,sim_sig_IntNoise.sigSRC_pre))
disp(~isequal(sim_sig_SrcActiv.sigSRC_pre,sim_sig_IntNoise.sigSRC_pst))

tmp_arrC = [rawSize(sim_sig_SrcActiv.sigSRC_pre,[1,2,3]); rawSize(sim_sig_SrcActiv.sigSRC_pst,[1,2,3]); rawSize(sim_sig_IntNoise.sigSRC_pre,[1,2,3]); rawSize(sim_sig_IntNoise.sigSRC_pst,[1,2,3])];
tmp_varN = {'Timesamples','Signals','Realizations'};
tmp_rowN = {'sim_sig_SrcActiv.sigSRC_pre','sim_sig_SrcActiv.sigSRC_pst','sim_sig_IntNoise.sigSRC_pre','sim_sig_IntNoise.sigSRC_pst'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));

tmp_arrC = [size(sim_sig_SrcActiv.A00),size(sim_sig_SrcActiv.A00,2)/size(sim_sig_SrcActiv.A00,1);size(sim_sig_SrcActiv.A01),size(sim_sig_SrcActiv.A01,2)/size(sim_sig_SrcActiv.A01,1);size(sim_sig_IntNoise.A00),size(sim_sig_IntNoise.A00,2)/size(sim_sig_IntNoise.A00,1);size(sim_sig_IntNoise.A01),size(sim_sig_IntNoise.A01,2)/size(sim_sig_IntNoise.A01,1)];
tmp_varN = {'Signals','Signals_cdot_ModOrd','ModOrd'};
tmp_rowN = {'sim_sig_SrcActiv.A00','sim_sig_SrcActiv.A01','sim_sig_IntNoise.A00','sim_sig_IntNoise.A01'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));

clearvars ii jj kk nn tmp*

       #+END_SRC

****** !fig MVAR model coefficient matrix (=A00=)

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

close all

figure(55);clf;set(gcf, 'Position', SETUP.DISP);
subplot(4,1,1);cccSim___rawDispA00(sim_sig_SrcActiv.A00);title('sim\_sig00.A00');
subplot(4,1,2);cccSim___rawDispA00(sim_sig_SrcActiv.A01);title('sim\_sig00.A01');
subplot(4,1,3);cccSim___rawDispA00(sim_sig_IntNoise.A00);title('sim\_sig30.A00');
subplot(4,1,4);cccSim___rawDispA00(sim_sig_IntNoise.A01);title('sim\_sig30.A01');
set(gcf,'color','w');

       #+END_SRC

****** !fig PDC

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

close all

figure(100);clf;set(gcf, 'Position', SETUP.DISP);cccSim___rawPlotA00(sim_sig_SrcActiv.A00);
set(gcf,'color','w');
figure(101);clf;set(gcf, 'Position', SETUP.DISP);cccSim___rawPlotA00(sim_sig_SrcActiv.A01);
set(gcf,'color','w');
figure(130);clf;set(gcf, 'Position', SETUP.DISP);cccSim___rawPlotA00(sim_sig_IntNoise.A00);
set(gcf,'color','w');
figure(131);clf;set(gcf, 'Position', SETUP.DISP);cccSim___rawPlotA00(sim_sig_IntNoise.A01);
set(gcf,'color','w');

       #+END_SRC

**** Generation of timeseries for bio-electrical background noise
***** !run Code

      #+NAME: ./tg_s01_Timeseries___03_BcgNoise.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s01_Timeseries___03_BcgNoise.m :results output silent

clearvars sim_sig_BcgNoise tmp*
sim_sig_BcgNoise = cccSim___makeSimSig(SETUP,3,'MVAR based signal for background (biological) noise activity sources');
clearvars ii jj kk nn tmp*

      #+END_SRC

***** !opt Checkups
****** !opt Basics 1

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos
sim_sig_SrcActiv
sim_sig_IntNoise
sim_sig_BcgNoise

chkSim___tg_s01_PRE___001(SETUP)

% Check size of signals generated
tmp_arrC = [rawSize(sim_sig_SrcActiv.sigSRC_pre,[1,2,3]); rawSize(sim_sig_SrcActiv.sigSRC_pst,[1,2,3]); rawSize(sim_sig_IntNoise.sigSRC_pre,[1,2,3]); rawSize(sim_sig_IntNoise.sigSRC_pst,[1,2,3]); rawSize(sim_sig_BcgNoise.sigSRC_pre,[1,2,3]); rawSize(sim_sig_BcgNoise.sigSRC_pst,[1,2,3])];
tmp_varN = {'Timesamples','Signals','Realizations'};
tmp_rowN = {'sim_sig_SrcActiv.sigSRC_pre','sim_sig_SrcActiv.sigSRC_pst','sim_sig_IntNoise.sigSRC_pre','sim_sig_IntNoise.sigSRC_pst','sim_sig_BcgNoise.sigSRC_pre','sim_sig_BcgNoise.sigSRC_pst'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));

       #+END_SRC

****** !opt Basics 2

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

jj = 1
kk = 1

tmp_arrC = corrcoef(sim_sig_SrcActiv.sigSRC_pre(:,jj,kk),sim_sig_BcgNoise.sigSRC_pre(:,jj,kk));
tmp_varN = {'sim_sig_SrcActiv_a','sim_sig_BcgNoise_a'};
tmp_rowN = tmp_varN;
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));disp('CYBERCRAFT:: This should be negligible')

       #+END_SRC

****** !fig MVAR model coefficient matrix (=A00=)

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

close all

figure(105);clf;set(gcf,'Position',SETUP.DISP);subplot(2,1,1);cccSim___rawDispA00(sim_sig_BcgNoise.A00);title('A00');subplot(2,1,2);cccSim___rawDispA00(sim_sig_BcgNoise.A01);title('A01');
set(gcf,'color','w');
figure(100);clf;set(gcf,'Position',SETUP.DISP);cccSim___rawPlotA00(sim_sig_SrcActiv.A01);
set(gcf,'color','w');
figure(130);clf;set(gcf,'Position',SETUP.DISP);cccSim___rawPlotA00(sim_sig_IntNoise.A01);
set(gcf,'color','w');
figure(140);clf;set(gcf,'Position',SETUP.DISP);cccSim___rawPlotA00(sim_sig_BcgNoise.A01);
set(gcf,'color','w');

       #+END_SRC

****** !fig PDC

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

close all

figure(55);clf;set(gcf, 'Position', SETUP.DISP);

subplot(6,1,1);cccSim___rawDispA00(sim_sig_SrcActiv.A00);title('sim\_sig\_SrcActiv.A00')
subplot(6,1,2);cccSim___rawDispA00(sim_sig_SrcActiv.A01);title('sim\_sig\_SrcActiv.A01')
subplot(6,1,3);cccSim___rawDispA00(sim_sig_IntNoise.A00);title('sim\_sig\_IntNoise.A00')
subplot(6,1,4);cccSim___rawDispA00(sim_sig_IntNoise.A01);title('sim\_sig\_IntNoise.A01')
subplot(6,1,5);cccSim___rawDispA00(sim_sig_BcgNoise.A00);title('sim\_sig\_BcgNoise.A00')
subplot(6,1,6);cccSim___rawDispA00(sim_sig_BcgNoise.A01);title('sim\_sig\_BcgNoise.A01')

       #+END_SRC

*** Measurement noise generation
**** !run Code

     #+NAME: ./tg_s01_Timeseries___04_MesNoise.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s01_Timeseries___04_MesNoise.m :results silent

clearvars sim_sig_MesNoise tmp*
sim_sig_MesNoise.inf = 'RANDN based signal for measurment noise';
sim_sig_MesNoise.sigSNS_pre = randn([SETUP.n00,size(sel_ele.chanpos,1),SETUP.K00]);
sim_sig_MesNoise.sigSNS_pst = randn([SETUP.n00,size(sel_ele.chanpos,1),SETUP.K00]);

     #+END_SRC

**** !opt Checkups

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos
sim_sig_SrcActiv
sim_sig_IntNoise
sim_sig_BcgNoise
sim_sig_MesNoise
tmp_arrC = [rawSize(sim_sig_SrcActiv.sigSRC_pre,[1,2,3]); rawSize(sim_sig_SrcActiv.sigSRC_pst,[1,2,3]); rawSize(sim_sig_IntNoise.sigSRC_pre,[1,2,3]); rawSize(sim_sig_IntNoise.sigSRC_pst,[1,2,3]); rawSize(sim_sig_BcgNoise.sigSRC_pre,[1,2,3]); rawSize(sim_sig_BcgNoise.sigSRC_pst,[1,2,3]); rawSize(sim_sig_MesNoise.sigSNS_pre,[1,2,3]); rawSize(sim_sig_MesNoise.sigSNS_pst,[1,2,3])];
tmp_varN = {'Timesamples','Signals','Realizations'};
tmp_rowN = {'sim_sig_SrcActiv.sigSRC_pre','sim_sig_SrcActiv.sigSRC_pst','sim_sig_IntNoise.sigSRC_pre','sim_sig_IntNoise.sigSRC_pst','sim_sig_BcgNoise.sigSRC_pre','sim_sig_BcgNoise.sigSRC_pst','sim_sig_BcgNoise.sigSNS_pre','sim_sig_BcgNoise.sigSNS_pst'};
disp(' ');disp(array2table(tmp_arrC,'VariableNames',tmp_varN,'RowNames',tmp_rowN));disp('NB: measurement noise should be already in electrode space not in the source space.')

chkSim___tg_s01_PRE___001(SETUP)
chkSim___tg_s02_PRE___001(SETUP)

     #+END_SRC

** DONE [s02] Geometry of head compartments, ROIs, source locations and leadfields
*** !opt Initial checkups
**** !fig Deep sources as icosahedron642

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

close all
figure('Color',[1,1,1]);clf;
sel_geo_deep_icosahedron642
trisurf(                   ...
    sel_geo_deep_icosahedron642.tri,      ...
    sel_geo_deep_icosahedron642.pnt(:,1), ...
    sel_geo_deep_icosahedron642.pnt(:,2), ...
    sel_geo_deep_icosahedron642.pnt(:,3), ...
    'facealpha',0.2,       ...
    'facecolor','m',       ...
    'edgecolor','m',       ...
    'edgealpha',0.4);
hold on
ccrender([-150,150],'finish','matte')

     #+END_SRC

**** !fig Deep sources as thalami (L+R)

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

close all
figure('Color',[1,1,1]);clf;
sel_geo_deep_thalami
trisurf(                   ...
    sel_geo_deep_thalami.tri,      ...
    sel_geo_deep_thalami.pnt(:,1), ...
    sel_geo_deep_thalami.pnt(:,2), ...
    sel_geo_deep_thalami.pnt(:,3), ...
    'facealpha',0.05,       ...
    'facecolor','m',       ...
    'edgecolor','m',       ...
    'edgealpha',0.1);
hold on
ccrender([-50,50],'finish','matte')

     #+END_SRC

**** !fig Cortex mesh

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

ft_plot_mesh(                  ...
    sel_atl,                   ...
    'facecolor',[0.9 0.9 0.9], ...
    'facealpha',0.0,           ...
    'edgecolor',[0.7 0.7 0.7], ...
    'edgealpha',0.2);
hold on
ccrender([-160,160],'finish','matte')
view(0,90)
view(90,0)
view(180,0)

     #+END_SRC

**** !fig Brain outer mesh

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

close all
figure('Color',[1,1,1]);clf;
ft_plot_mesh(        ...
    sel_msh.bnd(1),  ...
    'facecolor','m', ...
    'facealpha',0.1, ...
    'edgecolor','m', ...
    'edgealpha',0.05);
ccrender([-160,160],'finish','matte')

     #+END_SRC

**** !fig Skull outer mesh

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

ft_plot_mesh(        ...
    sel_msh.bnd(2),  ...
    'facecolor','c', ...
    'facealpha',0.1, ...
    'edgecolor','c', ...
    'edgealpha',0.1);
ccrender([-160,160],'finish','matte')

     #+END_SRC

**** !fig Scalp outer mesh

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

ft_plot_mesh(           ...
    sel_msh.bnd(3),     ...
    'facecolor','skin', ...
    'facealpha',0.1,    ...
    'edgecolor','red',  ...
    'edgealpha',0.2);
ccrender([-160,160],'finish','matte')

     #+END_SRC

**** !fig Electrode positioning

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

stem3( sel_ele.elecpos(:,1), ...
       sel_ele.elecpos(:,2), ...
       sel_ele.elecpos(:,3), ...
       '.m','filled','LineStyle','none','MarkerSize',12);

     #+END_SRC

**** !fig Electrode labels

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

% for ii=1:size(tmp_elecpos,1),sel_elec00.label{ii,1} = sprintf('Ch%02d',ii);end

text(  sel_ele.elecpos(:,1) * 1.2, ...
       sel_ele.elecpos(:,2) * 1.2, ...
       sel_ele.elecpos(:,3) * 1.2, ...
       sel_ele.label,'Color',[0.5 0.3 0.5],'FontSize',10);

     #+END_SRC

**** !fig View control

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

view(0,90)
view(90,0)
view(180,0)
view(225,60)

     #+END_SRC

*** !run ROI random sampling
**** !run Code

     #+NAME: ./tg_s02_Geometry___01_RandSamp.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s02_Geometry___01_RandSamp.m :results silent

clearvars sim_geo_cort ii jj kk nn tmp*
sim_geo_cort.numROIs = size(sel_atl.Atlas(sel_atl.atl).Scouts,2);

if SETUP.rROI
    disp('CC: using random ROI locations')
    sim_geo_cort.lstROIs = randsample([1:sim_geo_cort.numROIs],size(SETUP.SRCS,1));
else
    disp('CC: using predefined ROI locations')
    sim_geo_cort.lstROIs = [31 30 34 29 56 51 84 53 83 58 57];
    sim_geo_cort.lstROIs = sim_geo_cort.lstROIs(1:size(SETUP.SRCS,1));
end


if SETUP.rPNT
    disp('CC: using random source locations')
else
    disp('CC: using predefined source locations')
    sim_geo_cort.lstROIs(1) = 31;
end

     #+END_SRC

**** !opt Checkups

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

SETUP.rROI
SETUP.rPNT

sim_geo_cort
sim_geo_cort.lstROIs

chkSim___tg_s01_PRE___001(SETUP)
chkSim___tg_s02_Geometry___01_RandSamp(sim_geo_cort,sel_atl)

     #+END_SRC

**** !fun Checkup functions
***** =chkSim___tg_s02_Geometry___01_RandSamp=

      #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/chkSim___tg_s02_Geometry___01_RandSamp.m :results silent

function chkSim___tg_s02_Geometry___01_RandSamp(sim_geo_cort,sel_atl)
    fprintf('\n');
    disp(cell2table([num2cell(sim_geo_cort.lstROIs)',{sel_atl.Atlas(sel_atl.atl).Scouts(sim_geo_cort.lstROIs).Label}'],'VariableNames',{'ROI_number','ROI_name'}));
    fprintf('\n');
end

      #+END_SRC

*** !run Indices for ROI vertices and triangles ON CORTEX
**** !run Code

     Select all *vertices* and all *triangles* for each ROI.

     #+NAME: ./tg_s02_Geometry___02_Indices.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s02_Geometry___02_Indices.m :results silent

sim_geo_cort.indROIs = [];
for ii = 1:length(sim_geo_cort.lstROIs)
    tmp_roi = sim_geo_cort.lstROIs(ii);
    sim_geo_cort.indROIs(ii).pntNum00 = sel_atl.Atlas(sel_atl.atl).Scouts(tmp_roi).Vertices';
    sim_geo_cort.indROIs(ii).triNum00 = find(all(ismember(sel_atl.tri,sim_geo_cort.indROIs(ii).pntNum00),2));
end
clearvars ii jj kk nn tmp*

     #+END_SRC

**** !opt Checkups

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

sim_geo_cort
sim_geo_cort.indROIs
ii = 1
sim_geo_cort.indROIs(ii)
sim_geo_cort.indROIs(ii).pntNum00
sim_geo_cort.indROIs(ii).triNum00

     #+END_SRC

**** !opt ROI visualization functions
***** [num-000] Whole cortex

      #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/chkSim___tg_s02_Geometry___02_Indices___Plot_000.m :results silent

function chkSim___tg_s02_Geometry___02_Indices___Plot_000(sel_atl)
    ft_plot_mesh(                   ...
        sel_atl,                  ...
        'facecolor',[0.95 0.95 0.95],  ...
        'facealpha',0.1,            ...
        'edgecolor',[0.85 0.85 0.85],  ...
        'edgealpha',0.2);
    hold on
    ccrender([-160,160],'finish','matte')
end

      #+END_SRC

***** [num-001] ROIs

      #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/chkSim___tg_s02_Geometry___02_Indices___Plot_001.m :results silent

function chkSim___tg_s02_Geometry___02_Indices___Plot_001(sel_atl,sim_geo_cort)
    hold on
    tmp_colorBar = lines(length(sim_geo_cort.lstROIs));
    for ii = 1:length(sim_geo_cort.lstROIs)
        trisurf(                      ...
            sel_atl.tri(sim_geo_cort.indROIs(ii).triNum00,:),...
            sel_atl.pnt(:,1)-0.1*sel_atl.vn1(:,1),  ...
            sel_atl.pnt(:,2)-0.1*sel_atl.vn1(:,2),  ...
            sel_atl.pnt(:,3)-0.1*sel_atl.vn1(:,3),  ...
            'facealpha',0.2,                            ...
            'facecolor',tmp_colorBar(ii,:),             ...
            'edgecolor',tmp_colorBar(ii,:),             ...
            'edgealpha',0.4);
        hold on
    end
    ccrender([-160,160],'finish','matte')

%    caxis([1 length(sim_geo_cort.lstROIs)]);
%    colormap(tmp_colorBar)
%    colorbar
end

      #+END_SRC

***** [num-002] Cortex sources

      #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/chkSim___tg_s02_Geometry___02_Indices___Plot_002.m :results silent

function chkSim___tg_s02_Geometry___02_Indices___Plot_002(sim_geo_cort)
    hold on
    ii = 1
    qh = quiver3(...
        sim_geo_cort.pos_orig{ii}(:,1),sim_geo_cort.pos_orig{ii}(:,2),sim_geo_cort.pos_orig{ii}(:,3),...
        sim_geo_cort.ori_orig{ii}(:,1),sim_geo_cort.ori_orig{ii}(:,2),sim_geo_cort.ori_orig{ii}(:,3),...
        1,...
        'color','k');
    set(qh,'linewidth',1);
    qh = quiver3(...
        sim_geo_cort.pos_pert{ii}(:,1),sim_geo_cort.pos_pert{ii}(:,2),sim_geo_cort.pos_pert{ii}(:,3),...
        sim_geo_cort.ori_pert{ii}(:,1),sim_geo_cort.ori_pert{ii}(:,2),sim_geo_cort.ori_pert{ii}(:,3),...
        1,':',...
        'color','k');
    set(qh,'linewidth',1);
    ii = 2
    qh = quiver3(...
        sim_geo_cort.pos_orig{ii}(:,1),sim_geo_cort.pos_orig{ii}(:,2),sim_geo_cort.pos_orig{ii}(:,3),...
        sim_geo_cort.ori_orig{ii}(:,1),sim_geo_cort.ori_orig{ii}(:,2),sim_geo_cort.ori_orig{ii}(:,3),...
        1,...
        'color','r');
    set(qh,'linewidth',1);
    qh = quiver3(...
        sim_geo_cort.pos_pert{ii}(:,1),sim_geo_cort.pos_pert{ii}(:,2),sim_geo_cort.pos_pert{ii}(:,3),...
        sim_geo_cort.ori_pert{ii}(:,1),sim_geo_cort.ori_pert{ii}(:,2),sim_geo_cort.ori_pert{ii}(:,3),...
        1,':',...
        'color','r');
    set(qh,'linewidth',1);
    ii = 3
    qh = quiver3(...
        sim_geo_cort.pos_orig{ii}(:,1),sim_geo_cort.pos_orig{ii}(:,2),sim_geo_cort.pos_orig{ii}(:,3),...
        sim_geo_cort.ori_orig{ii}(:,1),sim_geo_cort.ori_orig{ii}(:,2),sim_geo_cort.ori_orig{ii}(:,3),...
        1,...
        'color','b');
    set(qh,'linewidth',1);
    qh = quiver3(...
        sim_geo_cort.pos_pert{ii}(:,1),sim_geo_cort.pos_pert{ii}(:,2),sim_geo_cort.pos_pert{ii}(:,3),...
        sim_geo_cort.ori_pert{ii}(:,1),sim_geo_cort.ori_pert{ii}(:,2),sim_geo_cort.ori_pert{ii}(:,3),...
        1,':',...
        'color','b');
    set(qh,'linewidth',1);
    ccrender([-160,160],'finish','matte')
end

      #+END_SRC

**** !opt ROI visualization

     Figure giving *ultimate situation view*.

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

close all
chkSim___tg_s01_PRE___001(SETUP)
chkSim___tg_s02_Geometry___01_RandSamp(sim_geo_cort,sel_atl)

chkSim___tg_s02_Geometry___02_Indices___Plot_000(sel_atl)
chkSim___tg_s02_Geometry___02_Indices___Plot_001(sel_atl,sim_geo_cort)

view(-115,55)
view(-115,-55)
view(0,60)
view(45,60)
view(90,60)
view(135,60)
view(180,60)
view(225,60)
view(270,60)
view(315,60)
view(315,-20)
view(135,45)

     #+END_SRC

*** !run Coordinates for source position and orientation
**** !run Code

     #+NAME: ./tg_s02_Geometry___03_Coordinates.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s02_Geometry___03_Coordinates.m :results silent

sim_geo_cort.distrROIs    = SETUP.SRCS;
sim_geo_cort.distrROIsSum = sum(SETUP.SRCS,2);
sim_geo_cort.bulkSRC = {};
for ii = 1:length(sim_geo_cort.distrROIsSum)
    sim_geo_cort.bulkSRC{ii,1} = randsample(sim_geo_cort.indROIs(ii).pntNum00,sim_geo_cort.distrROIsSum(ii));
end
sim_geo_cort.splitSRC = {};
for ii = 1:length(sim_geo_cort.distrROIsSum)
    sim_geo_cort.splitSRC(ii,:) = mat2cell(sim_geo_cort.bulkSRC{ii},sim_geo_cort.distrROIs(ii,:),1)';
end

if 0
    sim_geo_cort.splitSRC(1,1)
    sim_geo_cort.splitSRC{1,1}
end

if SETUP.rPNT
    disp('CC: using random source locations')
else
    disp('CC: using predefined source locations')
    sim_geo_cort.splitSRC{1,1} = [5441;5506;5481;5283;5144;5063;5823;5987;6154;6065;6166;6346;6087;6446;6367;6726;6859;6829;6613;6698;6966;7162;7227;6949;6996;7365;7247;7480];
    sim_geo_cort.splitSRC{1,1} = sim_geo_cort.splitSRC{1,1}(1:SETUP.SRCS(1,1));
    if 0
        sim_geo_cort.splitSRC{1,1}
    end
end

sim_geo_cort.mergeSRC = {};
for ii = 1:3
    sim_geo_cort.mergeSRC{ii} = cat(1,sim_geo_cort.splitSRC{:,ii});
end
sim_geo_cort.pos_orig = {};
sim_geo_cort.ori_orig = {};
sim_geo_cort.pos_pert = {};
sim_geo_cort.ori_pert = {};
for ii = 1:3
    sim_geo_cort.pos_orig{ii} = sel_atl.pnt(sim_geo_cort.mergeSRC{ii},:);
    sim_geo_cort.ori_orig{ii} = sel_atl.vn1(sim_geo_cort.mergeSRC{ii},:);
    sim_geo_cort.pos_pert{ii} = sim_geo_cort.pos_orig{ii};
    sim_geo_cort.ori_pert{ii} = sim_geo_cort.ori_orig{ii};
end
clearvars ii jj kk nn tmp*

     #+END_SRC

**** !opt Checkups

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

SETUP.SRCS
SETUP.DEEP

sim_geo_cort
sim_geo_cort.distrROIs
sim_geo_cort.distrROIsSum
sim_geo_cort.bulkSRC
sim_geo_cort.splitSRC
sim_geo_cort.distrROIs
sim_geo_cort.mergeSRC
sim_geo_cort.mergeSRC{1}

chkSim___tg_s02_Geometry___01_RandSamp(sim_geo_cort,sel_atl)
chkSim___tg_s01_PRE___001(SETUP)
sim_geo_cort.mergeSRC

sim_geo_cort.pos_orig
sim_geo_cort.ori_orig
sim_geo_cort.pos_pert
sim_geo_cort.ori_pert

     #+END_SRC

*** !run Deep sources
**** !run Code

     #+NAME: ./tg_s02_Geometry___04_DeepSrc_01.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s02_Geometry___04_DeepSrc.m :results silent

clearvars sel_geo_deep ii jj kk ll nn tmp*

sim_geo_deep = sel_geo_deep_icosahedron642;
sim_geo_deep = sel_geo_deep_thalami;

sim_geo_deep.bulkSRC =  randsample(1:size(sim_geo_deep.pnt,1),sum(SETUP.DEEP))';


[sim_geo_deep.mergeSRC{1},sim_geo_deep.mergeSRC{2},sim_geo_deep.mergeSRC{3}] = deal([]);


sim_geo_deep.mergeSRC{1,1} = sim_geo_deep.bulkSRC(1:SETUP.DEEP(1));
sim_geo_deep.mergeSRC{1,2} = sim_geo_deep.bulkSRC(1+SETUP.DEEP(1):SETUP.DEEP(1)+SETUP.DEEP(2));
sim_geo_deep.mergeSRC{1,3} = sim_geo_deep.bulkSRC(1+SETUP.DEEP(1)+SETUP.DEEP(2):SETUP.DEEP(1)+SETUP.DEEP(2)+SETUP.DEEP(3));



sim_geo_deep.pos_orig = {};
sim_geo_deep.ori_orig = {};
sim_geo_deep.pos_pert = {};
sim_geo_deep.ori_pert = {};

for ii = 1:3
    sim_geo_deep.pos_orig{ii} = sim_geo_deep.pnt(sim_geo_deep.mergeSRC{ii},:);
    sim_geo_deep.ori_orig{ii} = sim_geo_deep.vn1(sim_geo_deep.mergeSRC{ii},:);
    sim_geo_deep.pos_pert{ii} = sim_geo_deep.pos_orig{ii};
    sim_geo_deep.ori_pert{ii} = sim_geo_deep.ori_orig{ii};
end

clearvars ii jj kk nn tmp*

     #+END_SRC

**** !opt Checkups

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

sim_geo_deep
trisurf(                   ...
    sim_geo_deep.tri,      ...
    sim_geo_deep.pnt(:,1), ...
    sim_geo_deep.pnt(:,2), ...
    sim_geo_deep.pnt(:,3), ...
    'facealpha',0.2,       ...
    'facecolor','m',       ...
    'edgecolor','m',       ...
    'edgealpha',0.4);
ccrender([-150,150],'finish','matte')

     #+END_SRC

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos sel_*
whos sim_*
whos sim_geo*
chkSim___tg_s01_PRE___001(SETUP)
disp(SETUP.DEEP)
[sim_geo_deep.bulkSRC,[sim_geo_deep.mergeSRC{1}; sim_geo_deep.mergeSRC{2};sim_geo_deep.mergeSRC{3}]]
sim_geo_deep
sim_geo_cort

     #+END_SRC

*** !run Merge cortex and deep sources
**** !run Code

     #+NAME: ./tg_s02_Geometry___04_DeepSrc_02.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s02_Geometry___04_DeepSrc.m :results silent

sim_geo = sim_geo_cort;
for ii = 1:3,
    sim_geo.pos_orig{ii} = [sim_geo.pos_orig{ii};sim_geo_deep.pos_orig{ii}];
    sim_geo.ori_orig{ii} = [sim_geo.ori_orig{ii};sim_geo_deep.ori_orig{ii}];
    sim_geo.pos_pert{ii} = [sim_geo.pos_pert{ii};sim_geo_deep.pos_pert{ii}];
    sim_geo.ori_pert{ii} = [sim_geo.ori_pert{ii};sim_geo_deep.ori_pert{ii}];
end;

     #+END_SRC

**** !opt Checkups

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos sim* sel*
sim_geo
sim_geo_deep
sim_geo_cort

     #+END_SRC

*** !run Perturbation for source position and orientation
**** !opt Modify perturbation parameters here

     To facilitate *bug tracking* and for *explanatory/illustrative
     purpose* it might be plausible to modify the perturbation
     parameters here.

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

SETUP.CUBE
SETUP.CUBE = 30
SETUP.CUBE
SETUP.CONE
SETUP.CONE = pi/16
SETUP.CONE

     #+END_SRC

**** !run Code

     #+NAME: ./tg_s02_Geometry___05_Perturbation.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s02_Geometry___05_Perturbation.m :results silent

for ii = 1:3
    sim_geo.ori_pert{ii} = normr(rawSphToCart(rawCartToSph(sim_geo.ori_orig{ii}) + [SETUP.CONE*0.5*(2*rand(size(sim_geo.ori_pert{ii},1),2)-1),zeros(size(sim_geo.ori_pert{ii},1),1)]));
end
tmp_count = 0;
for ii = 1:3
    for jj = 1:size(sim_geo.pos_pert{ii},1)
        tmp_srcInside = false;
        while ~tmp_srcInside
            tmp_count = tmp_count+1;
            sim_geo.pos_pert{ii}(jj,:) = sim_geo.pos_orig{ii}(jj,:)+SETUP.CUBE*0.5*(2*rand([1,3])-1);
            tmp_srcInside = bounding_mesh(sim_geo.pos_pert{ii}(jj,:),sel_msh.bnd(1).pnt,sel_msh.bnd(1).tri);
        end
    end
end
if SETUP.TELL
    disp(['CYBERCRAFT:: Perturbation took ',tmp_count,' iterations'])
end
clearvars ii jj kk nn tmp*

     #+END_SRC

     NOTE 1: The following procedure is applied to the orientation vectors of activity sources:
     in order to achieve their random perturbation:
     1. Original orientations are vertex vormal vectors of the triangulation mesh;
     2. Coordinates are transformed to spherical;
     3. =SETUP.CONE= is the max perturbation (angle of max rotation) considered
     4. Multiplied here by random number \in [-1,1];
     5. Azimuth [TH] and elevation [PHI] are modified by adding the above random numbers;
     6. Radius is not changed (the added matrix has last column padded with zeros);
     7. Coordinates are transformed back to the cartesian system.

**** !opt Checkups

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

ii = 1

SETUP.CUBE
min(sim_geo.pos_orig{ii}-sim_geo.pos_pert{ii})
max(sim_geo.pos_orig{ii}-sim_geo.pos_pert{ii})

SETUP.CONE
min(rawCartToSph(sim_geo.ori_orig{ii})-rawCartToSph(sim_geo.ori_pert{ii}))
max(rawCartToSph(sim_geo.ori_orig{ii})-rawCartToSph(sim_geo.ori_pert{ii}))

sim_geo
whos

     #+END_SRC

**** !opt Source visualization

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

close all
% figure(1)
figure('Color',[1 1 1]);
clf

chkSim___tg_s01_PRE___001(SETUP)
chkSim___tg_s02_Geometry___01_RandSamp(sim_geo,sel_atl)

% mesh for ROIs on cortex
chkSim___tg_s02_Geometry___02_Indices___Plot_001(sel_atl,sim_geo)

% mesh for deep sources ROI
chkSim___tg_s02_Geometry___02_Indices___Plot_003(sim_geo_deep)

% sources
chkSim___tg_s02_Geometry___02_Indices___Plot_002(sim_geo)
legend('SrcActiv\_orig','SrcActiv\_pert','IntNoise\_orig','IntNoise\_pert','BcgNoise\_orig','BcgNoise\_pert')

% chkSim___tg_s02_Geometry___02_Indices___Plot_000(sel_atl)

view(0,90)
view(90,0)
view(180,0)

view(-115,55)
view(-115,-55)
view(0,60)
view(45,60)
view(90,60)
view(135,60)
view(180,60)
view(225,60)
view(270,60)
view(315,60)
view(315,-20)
view(135,45)

     #+END_SRC

**** !fun Source visualization functions
***** [num-003] Deep sources

      This one works for deep sources!!!

      #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/chkSim___tg_s02_Geometry___02_Indices___Plot_003.m :results silent

function chkSim___tg_s02_Geometry___02_Indices___Plot_003(sim_geo_deep)
    hold on
    trisurf(                   ...
        sim_geo_deep.tri,      ...
        sim_geo_deep.pnt(:,1), ...
        sim_geo_deep.pnt(:,2), ...
        sim_geo_deep.pnt(:,3), ...
        'facealpha',0.1,       ...
        'facecolor','m',       ...
        'edgecolor','m',       ...
        'edgealpha',0.2);
    ccrender([-150,150],'finish','matte')
end

      #+END_SRC

*** !run Leadfields computation

    Please note that if meshes for =brain=, =skull= and =scalp= change
    the =headmodel= also needs to be recalculated

    Leadfields are named with reference to signals.

    | Variable name         | Description                                 |
    |-----------------------+---------------------------------------------|
    | sim_lfg_SrcActiv_orig | Original activity of interest (biological)  |
    | sim_lfg_IntNoise_orig | Original interference (biological) noise    |
    | sim_lfg_BcgNoise_orig | Original background (biological) noise      |
    | sim_lfg_SrcActiv_pert | Perturbed activity of interest (biological) |
    | sim_lfg_IntNoise_pert | Perturbed interference (biological) noise   |
    | sim_lfg_BcgNoise_pert | Perturbed background (biological) noise     |


    #+NAME: ./tg_s02_Geometry___06_lfg.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s02_Geometry___06_lfg.m :results silent

% WARNING:: suppress warning about MATLAB version and CFG tracking
if 0
    TMP_S = warning('off','all');
end

tmp_cfg            = [];
tmp_cfg.grid.unit  = 'mm';
tmp_cfg.reducerank = 'no';
tmp_cfg.normalize  = 'no';
tmp_cfg.elec       = sel_ele;
tmp_cfg.vol        = sel_vol;

ii = 1;
tmp_cfg.grid.pos          = sim_geo.pos_orig{ii};
tmp_cfg.grid.mom          = transpose(sim_geo.ori_orig{ii});
sim_lfg_SrcActiv_orig.lfg = ft_prepare_leadfield(tmp_cfg);
sim_lfg_SrcActiv_orig.LFG = cat(2,sim_lfg_SrcActiv_orig.lfg.leadfield{:});

tmp_cfg.grid.pos          = sim_geo.pos_pert{ii};
tmp_cfg.grid.mom          = transpose(sim_geo.ori_pert{ii});
sim_lfg_SrcActiv_pert.lfg = ft_prepare_leadfield(tmp_cfg);
sim_lfg_SrcActiv_pert.LFG = cat(2,sim_lfg_SrcActiv_pert.lfg.leadfield{:});

ii = 2;
tmp_cfg.grid.pos          = sim_geo.pos_orig{ii};
tmp_cfg.grid.mom          = transpose(sim_geo.pos_orig{ii});
sim_lfg_IntNoise_orig.lfg = ft_prepare_leadfield(tmp_cfg);
sim_lfg_IntNoise_orig.LFG = cat(2,sim_lfg_IntNoise_orig.lfg.leadfield{:});

tmp_cfg.grid.pos          = sim_geo.pos_pert{ii};
tmp_cfg.grid.mom          = transpose(sim_geo.pos_pert{ii});
sim_lfg_IntNoise_pert.lfg = ft_prepare_leadfield(tmp_cfg);
sim_lfg_IntNoise_pert.LFG = cat(2,sim_lfg_IntNoise_pert.lfg.leadfield{:});

ii = 3;
tmp_cfg.grid.pos          = sim_geo.pos_orig{ii};
tmp_cfg.grid.mom          = transpose(sim_geo.pos_orig{ii});
sim_lfg_BcgNoise_orig.lfg = ft_prepare_leadfield(tmp_cfg);
sim_lfg_BcgNoise_orig.LFG = cat(2,sim_lfg_BcgNoise_orig.lfg.leadfield{:});

tmp_cfg.grid.pos          = sim_geo.pos_pert{ii};
tmp_cfg.grid.mom          = transpose(sim_geo.pos_pert{ii});
sim_lfg_BcgNoise_pert.lfg = ft_prepare_leadfield(tmp_cfg);
sim_lfg_BcgNoise_pert.LFG = cat(2,sim_lfg_BcgNoise_pert.lfg.leadfield{:});

clearvars ii jj kk nn tmp*

    #+END_SRC

** DONE [s03] Forward modeling
*** !run Code

    #+NAME: ./tg_s03_Forward_modeling.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s03_Forward_modeling.m :results silent

clearvars ii jj kk nn tmp*
for kk = 1:SETUP.K00
    sim_sig_SrcActiv.sigSNS_pre(:,:,kk) = (sim_lfg_SrcActiv_orig.LFG*sim_sig_SrcActiv.sigSRC_pre(:,:,kk)')';
    sim_sig_SrcActiv.sigSNS_pst(:,:,kk) = (sim_lfg_SrcActiv_orig.LFG*sim_sig_SrcActiv.sigSRC_pst(:,:,kk)')';
    sim_sig_IntNoise.sigSNS_pre(:,:,kk) = (sim_lfg_IntNoise_orig.LFG*sim_sig_IntNoise.sigSRC_pre(:,:,kk)')';
    sim_sig_IntNoise.sigSNS_pst(:,:,kk) = (sim_lfg_IntNoise_orig.LFG*sim_sig_IntNoise.sigSRC_pst(:,:,kk)')';
    sim_sig_BcgNoise.sigSNS_pre(:,:,kk) = (sim_lfg_BcgNoise_orig.LFG*sim_sig_BcgNoise.sigSRC_pre(:,:,kk)')';
    sim_sig_BcgNoise.sigSNS_pst(:,:,kk) = (sim_lfg_BcgNoise_orig.LFG*sim_sig_BcgNoise.sigSRC_pst(:,:,kk)')';
end
clearvars ii jj kk nn tmp*

    #+END_SRC
** DONE [s04] Preparations for both reconstruction and localization
*** DONE SNRs adjustments
**** !opt Comments

     We introduced desired SNR by adjusting power for:
     - biological interference,
     - background biological activity,
     - measurement noise.
       In order to allow for evaluation of:
       - reconstruction accuracy, and
       - functional dependencies estimation (PDC, DTF, dDTF etc.)
       we did not change the power of the sources signal!

**** !run Code

     #+NAME: ./tg_s04_Preparation___00_SNRs_Adjustment.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s04_Preparation___00_SNRs_Adjustment.m :results silent

clearvars ii jj kk nn tmp*

sim_sig_AdjSNRs.SrcActivPre = sim_sig_SrcActiv.sigSNS_pre;
sim_sig_AdjSNRs.SrcActivPst = sim_sig_SrcActiv.sigSNS_pst;
sim_sig_AdjSNRs.IntNoisePre = sim_sig_IntNoise.sigSNS_pre;
sim_sig_AdjSNRs.IntNoisePst = sim_sig_IntNoise.sigSNS_pst;
sim_sig_AdjSNRs.BcgNoisePre = sim_sig_BcgNoise.sigSNS_pre;
sim_sig_AdjSNRs.BcgNoisePst = sim_sig_BcgNoise.sigSNS_pst;
sim_sig_AdjSNRs.MesNoisePre = sim_sig_MesNoise.sigSNS_pre;
sim_sig_AdjSNRs.MesNoisePst = sim_sig_MesNoise.sigSNS_pst;
for kk = 1:SETUP.K00
    sim_sig_AdjSNRs.IntNoisePre(:,:,kk) = rawAdjTotSNRdB(sim_sig_AdjSNRs.SrcActivPre(:,:,kk),sim_sig_AdjSNRs.IntNoisePre(:,:,kk),SETUP.SINR);
    sim_sig_AdjSNRs.IntNoisePst(:,:,kk) = rawAdjTotSNRdB(sim_sig_AdjSNRs.SrcActivPst(:,:,kk),sim_sig_AdjSNRs.IntNoisePst(:,:,kk),SETUP.SINR);
    sim_sig_AdjSNRs.BcgNoisePre(:,:,kk) = rawAdjTotSNRdB(sim_sig_AdjSNRs.SrcActivPre(:,:,kk),sim_sig_AdjSNRs.BcgNoisePre(:,:,kk),SETUP.SBNR);
    sim_sig_AdjSNRs.BcgNoisePst(:,:,kk) = rawAdjTotSNRdB(sim_sig_AdjSNRs.SrcActivPst(:,:,kk),sim_sig_AdjSNRs.BcgNoisePst(:,:,kk),SETUP.SBNR);
    sim_sig_AdjSNRs.MesNoisePre(:,:,kk) = rawAdjTotSNRdB(sim_sig_AdjSNRs.SrcActivPre(:,:,kk),sim_sig_AdjSNRs.MesNoisePre(:,:,kk),SETUP.SMNR);
    sim_sig_AdjSNRs.MesNoisePst(:,:,kk) = rawAdjTotSNRdB(sim_sig_AdjSNRs.SrcActivPst(:,:,kk),sim_sig_AdjSNRs.MesNoisePst(:,:,kk),SETUP.SMNR);
end

clearvars ii jj kk nn tmp*

     #+END_SRC

**** !opt Checkups

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos

SETUP

sim_sig_SrcActiv % signal
sim_sig_IntNoise % interference
sim_sig_BcgNoise % background
sim_sig_MesNoise % measurment
sim_sig_AdjSNRs % merged signals (bulk) with SNR adjusted
kk = 1
2*pow2db(rawNrm(sim_sig_AdjSNRs.SrcActivPre(:,:,kk)))
2*pow2db(rawNrm(sim_sig_AdjSNRs.SrcActivPst(:,:,kk)))
rawSNRdB(sim_sig_AdjSNRs.SrcActivPre(:,:,kk),sim_sig_AdjSNRs.IntNoisePre(:,:,kk))
rawSNRdB(sim_sig_AdjSNRs.SrcActivPst(:,:,kk),sim_sig_AdjSNRs.IntNoisePst(:,:,kk))
rawSNRdB(sim_sig_AdjSNRs.SrcActivPre(:,:,kk),sim_sig_AdjSNRs.BcgNoisePre(:,:,kk))
rawSNRdB(sim_sig_AdjSNRs.SrcActivPst(:,:,kk),sim_sig_AdjSNRs.BcgNoisePst(:,:,kk))
rawSNRdB(sim_sig_AdjSNRs.SrcActivPre(:,:,kk),sim_sig_AdjSNRs.MesNoisePre(:,:,kk))
rawSNRdB(sim_sig_AdjSNRs.SrcActivPst(:,:,kk),sim_sig_AdjSNRs.MesNoisePst(:,:,kk))

% sim_sig_SrcActiv.sigSNS_pre should remain unchanged after SNR adjustments
isequal(sim_sig_AdjSNRs.SrcActivPre,sim_sig_SrcActiv.sigSNS_pre)
% sim_sig_SrcActiv.sigSNS_pre should remain unchanged after SNR adjustments
isequal(sim_sig_AdjSNRs.SrcActivPst,sim_sig_SrcActiv.sigSNS_pst)

     #+END_SRC

*** DONE Measured signal and covariance matrices
**** !opt Initial checkups

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos

     #+END_SRC

**** !run Code
***** Initialization
****** !run Code

       #+NAME: ./tg_s04_Preparation___01_Signal_Covariances.m
       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s04_Preparation___01_Signal_Covariances.m :results silent

clearvars ii jj kk nn tmp*

       #+END_SRC

***** Computation
****** !run Code

       #+NAME: ./tg_s04_Preparation___01_Signal_Covariances.m
       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s04_Preparation___01_Signal_Covariances.m :results silent

% clearvars ii jj kk nn tmp*

% Note: in this and subsequent subsections we aim at following closely notation from the paper.

if SETUP.H_Src_pert
    H_Src = sim_lfg_SrcActiv_pert.LFG;
else
    H_Src = sim_lfg_SrcActiv_orig.LFG;
end

if SETUP.H_Int_pert
    H_Int = sim_lfg_IntNoise_pert.LFG;
else
    H_Int = sim_lfg_IntNoise_orig.LFG;
end

% reduced-rank leadfield for patch constraints
[U_H_Int Si_H_Int V_H_Int] = svd(H_Int);

for ii = SETUP.IntLfgRANK+1:size(H_Int,2)
    Si_H_Int(ii,ii)=0;
end

H_Int = U_H_Int*Si_H_Int*V_H_Int';

y_Pre =     SETUP.SigPre*sim_sig_AdjSNRs.SrcActivPre ...
    + SETUP.IntPre*sim_sig_AdjSNRs.IntNoisePre ...
    + SETUP.BcgPre*sim_sig_AdjSNRs.BcgNoisePre ...
    + SETUP.MesPre*sim_sig_AdjSNRs.MesNoisePre;

y_Pst =     SETUP.SigPst*sim_sig_AdjSNRs.SrcActivPst ...
    + SETUP.IntPst*sim_sig_AdjSNRs.IntNoisePst ...
    + SETUP.BcgPst*sim_sig_AdjSNRs.BcgNoisePst ...
    + SETUP.MesPst*sim_sig_AdjSNRs.MesNoisePst;

% background activity (as in Pre segment) in Pst segment
y_PstNoise =     SETUP.SigPre*sim_sig_AdjSNRs.SrcActivPst ...
    + SETUP.IntPre*sim_sig_AdjSNRs.IntNoisePst ...
    + SETUP.BcgPre*sim_sig_AdjSNRs.BcgNoisePst ...
    + SETUP.MesPre*sim_sig_AdjSNRs.MesNoisePst;

N = cov(reshape(permute(y_Pre,[1 3 2]),[],size(y_Pre,2),1));
% N = cov(reshape(permute(y_PstNoise,[1 3 2]),[],size(y_PstNoise,2),1));
R = cov(reshape(permute(y_Pst,[1 3 2]),[],size(y_Pst,2),1));

G = H_Src'*pinv(N)*H_Src;
S = H_Src'*pinv(R)*H_Src;

G_SrcInt = [H_Src H_Int]'*pinv(N)*[H_Src H_Int];
S_SrcInt = [H_Src H_Int]'*pinv(R)*[H_Src H_Int];
C_SrcInt = pinv(S_SrcInt) - pinv(G_SrcInt);

% Estimated C for both regular and nulling filters
C_NL = C_SrcInt(1:size(H_Src,2),1:size(H_Src,2));
C = pinv(S)-pinv(G);

clear U_H_Int Si_H_Int V_H_Int G_SrcInt S_SrcInt;

       #+END_SRC

****** !opt Checkups 1

       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos
size(H_Src)
size(H_Int)
size(y_Pre)
size(y_Pst)
size(N)
size(R)
size(G)
size(S)
size(C)

% mean values of signals
y_Pst_RSH = reshape(permute(y_Pst,[1 3 2]),[],size(y_Pst,2),1);
max(abs(mean(y_Pst_RSH)))

clear y_Pst_RSH;

SrcPst = reshape(permute(sim_sig_AdjSNRs.SrcActivPst,[1 3 2]),[],size(sim_sig_AdjSNRs.SrcActivPst,2),1);
max(abs(mean(SrcPst)))

clear SrcPst;

IntPst = reshape(permute(sim_sig_AdjSNRs.IntNoisePst,[1 3 2]),[],size(sim_sig_AdjSNRs.IntNoisePst,2),1);
max(abs(mean(IntPst)))

clear IntPst;

BcgPst = reshape(permute(sim_sig_AdjSNRs.BcgNoisePst,[1 3 2]),[],size(sim_sig_AdjSNRs.BcgNoisePst,2),1);
max(abs(mean(BcgPst)))

clear BcgPst;

BcgPre = reshape(permute(sim_sig_AdjSNRs.BcgNoisePre,[1 3 2]),[],size(sim_sig_AdjSNRs.BcgNoisePre,2),1);
% note: max and min values vector across electrodes
max(abs(BcgPre))'
min(abs(BcgPre))'

clear BcgPre;

MesPst = reshape(permute(sim_sig_AdjSNRs.MesNoisePst,[1 3 2]),[],size(sim_sig_AdjSNRs.MesNoisePst,2),1);
max(abs(mean(MesPst)))

clear MesPst;

% correlation between signals check and # of highly correlated electrodes for Src, Bcg, Mes signals
SrcPst = reshape(permute(sim_sig_AdjSNRs.SrcActivPst,[1 3 2]),[],size(sim_sig_AdjSNRs.SrcActivPst,2),1);
BcgPst = reshape(permute(sim_sig_AdjSNRs.BcgNoisePst,[1 3 2]),[],size(sim_sig_AdjSNRs.BcgNoisePst,2),1);
Corr_SrcPst = corrcoef(SrcPst);

for ii = 1:length(Corr_SrcPst)
    [a k]=find(Corr_SrcPst(:,ii)>.9);
    x(ii)=sum(k);
end

x
x=[];
Corr_BcgPst = corrcoef(BcgPst);

for ii = 1:length(Corr_BcgPst)
    [a k]=find(Corr_BcgPst(:,ii)>.9);
    x(ii)=sum(k);
end

x
x=[];
Corr_SrcBcgPst = corrcoef([SrcPst BcgPst]);
Cross_SrcBcgPst = Corr_SrcBcgPst(length(Corr_SrcPst)+1:end,1:length(Corr_SrcPst));
max([max(Cross_SrcBcgPst) abs(min(Cross_SrcBcgPst))])

clear Corr_SrcPst Corr_BcgPst Corr_SrcBcgPst Cross_SrcBcgPst x;

MesPst = reshape(permute(sim_sig_AdjSNRs.MesNoisePst,[1 3 2]),[],size(sim_sig_AdjSNRs.MesNoisePst,2),1);
Corr_SrcPst = corrcoef(SrcPst);

for ii = 1:length(Corr_SrcPst)
    [a k]=find(Corr_SrcPst(:,ii)>.9);
    x(ii)=sum(k);
end

x
x=[];
Corr_MesPst = corrcoef(MesPst);

for ii = 1:length(Corr_MesPst)
    [a k]=find(Corr_MesPst(:,ii)>.9);
    x(ii)=sum(k);
end

x
x=[];
Corr_SrcMesPst = corrcoef([SrcPst MesPst]);
Cross_SrcMesPst = Corr_SrcMesPst(length(Corr_SrcPst)+1:end,1:length(Corr_SrcPst));
max([max(Cross_SrcMesPst) abs(min(Cross_SrcMesPst))])

clear Corr_MesPst Corr_SrcMesPst Cross_SrcMesPst x;

Corr_BcgMesPst = corrcoef([BcgPst MesPst]);
Cross_BcgMesPst = Corr_BcgMesPst(length(Corr_SrcPst)+1:end,1:length(Corr_SrcPst));
max([max(Cross_BcgMesPst) abs(min(Cross_BcgMesPst))])

clear BcgPst MesPst Corr_BcgMesPst Cross_BcgMesPst;

IntPst = reshape(permute(sim_sig_AdjSNRs.IntNoisePst,[1 3 2]),[],size(sim_sig_AdjSNRs.IntNoisePst,2),1);
Corr_SrcIntPst = corrcoef([SrcPst IntPst]);
Cross_SrcIntPst = Corr_SrcIntPst(length(Corr_SrcPst)+1:end,1:length(Corr_SrcPst));
max([max(Cross_SrcIntPst) abs(min(Cross_SrcIntPst))])

clear ii a k SrcPst Corr_SrcPst IntPst Corr_SrcIntPst Cross_SrcIntPst;

       #+END_SRC

**** !opt Checkups

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

rec_opt
whos

     #+END_SRC

** DONE [s05] Spatial filtering
*** DONE Filter preparations
***** NL

      #+NAME: ./tg_s05_Spatial_filtering___01_Preparation.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___01_Preparation.m :results silent

%  NL filter as proposed in Hui & Leahy 2010
H_SrcInt = [H_Src H_Int];
P_NL = [eye(size(H_Src,2)) zeros(size(H_Src,2),size(H_Int,2))];

      #+END_SRC
      
***** EIG-LCMV

      #+NAME: ./tg_s05_Spatial_filtering___01_Preparation.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___01_Preparation.m :results silent

%  EIG-LCMV filter
[UR, SiR, ~] = svd(R);
P_EIG = UR(:,1:SETUP.RANK_EIG)*UR(:,1:SETUP.RANK_EIG)';

clear UR SiR;

      #+END_SRC

***** sMVP
****** Initialization
       #+NAME: ./tg_s05_Spatial_filtering___01_Preparation.m
       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___01_Preparation.m :results silent

%  sMVP filters
H_Src_R = pinv(sqrtm(R))*H_Src;
H_Src_N = pinv(sqrtm(N))*H_Src;

       #+END_SRC

****** sMVP_MSE

       #+NAME: ./tg_s05_Spatial_filtering___01_Preparation.m
       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___01_Preparation.m :results silent

%  sMVP_MSE filter
% note: S and G use only H_Src, as we want to recover only activity of interest, not interference
% note: pinv(H_Src_R)*pinv(H_Src_R)' = pinv(S)
K_MSE = pinv(S)-2*C;
[U_K_MSE, W_K_MSE] = eig(K_MSE);
[~, p_K_MSE] = sort(diag(W_K_MSE));
% note: for selecting optimal rank, Theobald's Th. can be used without evaluation of the cost function as is done here
U_K_MSE = U_K_MSE(:,p_K_MSE);
sMVP_MSE_ranks = zeros(1,size(H_Src,2));

for jj = 1:size(H_Src,2)
    sMVP_MSE_ranks(jj) = trace(U_K_MSE(:,1:jj)*U_K_MSE(:,1:jj)'*K_MSE);
end

[~, sMVP_MSE_rank_opt] = min(sMVP_MSE_ranks);
P_sMVP_MSE_opt = U_K_MSE(:,1:sMVP_MSE_rank_opt)*U_K_MSE(:,1:sMVP_MSE_rank_opt)';

rec_opt.ranks.sMVP_MSE = sMVP_MSE_rank_opt;
clear K_MSE U_K_MSE W_K_MSE p_K_MSE sMVP_MSE_ranks sMVP_MSE_rank_opt;

       #+END_SRC

****** sMVP_R

       #+NAME: ./tg_s05_Spatial_filtering___01_Preparation.m
       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___01_Preparation.m :results silent

%  sMVP_R filter
% note: S and G use only H_Src, as we want to recover only activity of interest, not interference
% note: pinv(H_Src_R)*pinv(H_Src_R)' = pinv(S)
K_MSE = pinv(S)-2*C;

K_R = K_MSE+2*C;
[U_K_R, W_K_R] = eig(K_R);
[~, p_K_R] = sort(diag(W_K_R));
U_K_R = U_K_R(:,p_K_R);
sMVP_R_ranks = zeros(1,size(H_Src,2));

for jj = 1:size(H_Src,2)
    sMVP_R_ranks(jj) = trace(U_K_R(:,1:jj)*U_K_R(:,1:jj)'*K_MSE);
end

[~, sMVP_R_rank_opt] = min(sMVP_R_ranks);
P_sMVP_R_opt = U_K_R(:,1:sMVP_R_rank_opt)*U_K_R(:,1:sMVP_R_rank_opt)';

rec_opt.ranks.sMVP_R = sMVP_R_rank_opt;
clear K_MSE K_R U_K_R W_K_R p_K_R sMVP_R_ranks;

       #+END_SRC

****** sMVP_N

       #+NAME: ./tg_s05_Spatial_filtering___01_Preparation.m
       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___01_Preparation.m :results silent

%  sMVP_N filter
% note: S and G use only H_Src, as we want to recover only activity of interest, not interference
% note: pinv(H_Src_N)*pinv(H_Src_N)' = pinv(G)
K_MSE_N = pinv(G)-C;

K_N = K_MSE_N+C;
[U_K_N, W_K_N] = eig(K_N);
[~, p_K_N] = sort(diag(W_K_N));
U_K_N = U_K_N(:,p_K_N);
sMVP_N_ranks = zeros(1,size(H_Src,2));

for jj = 1:size(H_Src,2)
    sMVP_N_ranks(jj) = trace(U_K_N(:,1:jj)*U_K_N(:,1:jj)'*K_MSE_N);
end

[~, sMVP_N_rank_opt] = min(sMVP_N_ranks);
P_sMVP_N_opt = U_K_N(:,1:sMVP_N_rank_opt)*U_K_N(:,1:sMVP_N_rank_opt)';

rec_opt.ranks.sMVP_N = sMVP_N_rank_opt;
clear K_MSE_N K_N U_K_N W_K_N p_K_N sMVP_N_ranks;

       #+END_SRC

****** Cleanup
       #+NAME: ./tg_s05_Spatial_filtering___01_Preparation.m
       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___01_Preparation.m :results silent

%  sMVP filters
clear H_Src_R H_Src_N;

       #+END_SRC

***** sMVP_NL
****** Initialization
       #+NAME: ./tg_s05_Spatial_filtering___01_Preparation.m
       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___01_Preparation.m :results silent

%  sMVP-NL filters
H_Src_R = pinv(sqrtm(R))*H_Src;
H_Src_N = pinv(sqrtm(N))*H_Src;

H_Int_R = pinv(sqrtm(R))*H_Int;
H_Int_N = pinv(sqrtm(N))*H_Int;

P_sMVP_NL_R = eye(size(H_Int_R,1))-H_Int_R*pinv(H_Int_R);
P_sMVP_NL_N = eye(size(H_Int_N,1))-H_Int_N*pinv(H_Int_N);

       #+END_SRC

****** sMVP_NL_MSE

       #+NAME: ./tg_s05_Spatial_filtering___01_Preparation.m
       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___01_Preparation.m :results silent

%  sMVP_NL_MSE filter
K_NL_MSE = pinv(P_sMVP_NL_R*H_Src_R)*P_sMVP_NL_R*(pinv(P_sMVP_NL_R*H_Src_R))'-2*C_NL;
[U_K_NL_MSE, W_K_NL_MSE] = eig(K_NL_MSE);
[~, p_K_NL_MSE] = sort(diag(W_K_NL_MSE));
% note: for selecting optimal rank, Theobald's Th. can be used without evaluation of the cost function as is done here
U_K_NL_MSE = U_K_NL_MSE(:,p_K_NL_MSE);
sMVP_NL_MSE_ranks = zeros(1,size(H_Src,2));

for jj = 1:size(H_Src,2)
    sMVP_NL_MSE_ranks(jj) = trace(U_K_NL_MSE(:,1:jj)*U_K_NL_MSE(:,1:jj)'*K_NL_MSE);
end

[~, sMVP_NL_MSE_rank_opt] = min(sMVP_NL_MSE_ranks);
P_sMVP_NL_MSE_opt = U_K_NL_MSE(:,1:sMVP_NL_MSE_rank_opt)*U_K_NL_MSE(:,1:sMVP_NL_MSE_rank_opt)';

rec_opt.ranks.sMVP_NL_MSE = sMVP_NL_MSE_rank_opt;
clear K_NL_MSE U_K_NL_MSE W_K_NL_MSE p_K_NL_MSE sMVP_NL_MSE_ranks sMVP_NL_MSE_rank_opt;

       #+END_SRC

****** sMVP_NL_R

       #+NAME: ./tg_s05_Spatial_filtering___01_Preparation.m
       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___01_Preparation.m :results silent

%  sMVP_NL_R filter
K_NL_MSE = pinv(P_sMVP_NL_R*H_Src_R)*P_sMVP_NL_R*(pinv(P_sMVP_NL_R*H_Src_R))'-2*C_NL;

%K_NL_R = pinv(P_sMVP_NL_R*H_Src_R)*P_sMVP_NL_R*(pinv(P_sMVP_NL_R*H_Src_R))';
K_NL_R = K_NL_MSE+2*C_NL;
[U_K_NL_R, W_K_NL_R] = eig(K_NL_R);
[~, p_K_NL_R] = sort(diag(W_K_NL_R));
U_K_NL_R = U_K_NL_R(:,p_K_NL_R);
sMVP_NL_R_ranks = zeros(1,size(H_Src,2));

for jj = 1:size(H_Src,2)
    sMVP_NL_R_ranks(jj) = trace(U_K_NL_R(:,1:jj)*U_K_NL_R(:,1:jj)'*K_NL_MSE);
end

[~, sMVP_NL_R_rank_opt] = min(sMVP_NL_R_ranks);
P_sMVP_NL_R_opt = U_K_NL_R(:,1:sMVP_NL_R_rank_opt)*U_K_NL_R(:,1:sMVP_NL_R_rank_opt)';

rec_opt.ranks.sMVP_NL_R = sMVP_NL_R_rank_opt;
clear K_NL_MSE K_NL_R U_K_NL_R W_K_NL_R p_K_NL_R sMVP_NL_R_ranks sMVP_NL_R_rank_opt;

       #+END_SRC

****** sMVP_NL_N

       #+NAME: ./tg_s05_Spatial_filtering___01_Preparation.m
       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___01_Preparation.m :results silent

%  sMVP_NL_N filter
K_NL_MSE_N = pinv(P_sMVP_NL_N*H_Src_N)*P_sMVP_NL_N*(pinv(P_sMVP_NL_N*H_Src_N))'-C_NL;

%K_NL_N = pinv(P_sMVP_NL_N*H_Src_N)*P_sMVP_NL_N*(pinv(P_sMVP_NL_N*H_Src_N))';
K_NL_N = K_NL_MSE_N+C_NL;
[U_K_NL_N, W_K_NL_N] = eig(K_NL_N);
[~, p_K_NL_N] = sort(diag(W_K_NL_N));
U_K_NL_N = U_K_NL_N(:,p_K_NL_N);
sMVP_NL_N_ranks = zeros(1,size(H_Src,2));

for jj = 1:size(H_Src,2)
    sMVP_NL_N_ranks(jj) = trace(U_K_NL_N(:,1:jj)*U_K_NL_N(:,1:jj)'*K_NL_MSE_N);
end

[~, sMVP_NL_N_rank_opt] = min(sMVP_NL_N_ranks);
% note the simplified notation for P_N_opt
P_sMVP_NL_N_opt = U_K_NL_N(:,1:sMVP_NL_N_rank_opt)*U_K_NL_N(:,1:sMVP_NL_N_rank_opt)';

rec_opt.ranks.sMVP_NL_N = sMVP_NL_N_rank_opt;
clear K_NL_MSE_N K_NL_N U_K_NL_N W_K_NL_N p_K_NL_N sMVP_NL_N_ranks sMVP_NL_N_rank_opt;

       #+END_SRC

****** Cleanup
       #+NAME: ./tg_s05_Spatial_filtering___01_Preparation.m
       #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___01_Preparation.m :results silent

%  sMVP-NL filters
clear H_Int_R H_Int_N;

       #+END_SRC

*** DONE Filter definitions
**** !run Code
***** LCMV

      #+NAME: ./tg_s05_Spatial_filtering___02_Definitions.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___02_Definitions.m :results silent
% note: S and G use only H_Src, as we want to recover only activity of interest, not interference 
rec_flt.LCMV_R = pinv(S)*H_Src'*pinv(R);
rec_flt.LCMV_N = pinv(G)*H_Src'*pinv(N);

      #+END_SRC

***** NL

      #+NAME: ./tg_s05_Spatial_filtering___02_Definitions.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___02_Definitions.m :results silent

% As proposed in Hui & Leahy 2010
rec_flt.NL = P_NL*pinv(H_SrcInt'*pinv(R)*H_SrcInt)*H_SrcInt'*pinv(R);

      #+END_SRC

***** MMSE

      #+NAME: ./tg_s05_Spatial_filtering___02_Definitions.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___02_Definitions.m :results silent

rec_flt.MMSE = C*H_Src'*pinv(R);

      #+END_SRC

***** MMSE_INT

     #+NAME: ./tg_s05_Spatial_filtering___02_Definitions.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___02_Definitions.m :results silent

rec_flt.MMSE_INT = C_SrcInt(1:size(H_Src,2),:)*H_SrcInt'*pinv(R);

     #+END_SRC

***** ZF

      #+NAME: ./tg_s05_Spatial_filtering___02_Definitions.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___02_Definitions.m :results silent

rec_flt.ZF = pinv(H_Src);

      #+END_SRC
      
***** RANDN

      #+NAME: ./tg_s05_Spatial_filtering___02_Definitions.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___02_Definitions.m :results silent

rec_flt.RANDN = randn(size(H_Src'));

      #+END_SRC

***** ZEROS

      #+NAME: ./tg_s05_Spatial_filtering___02_Definitions.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___02_Definitions.m :results silent

rec_flt.ZEROS = zeros(size(H_Src'));

      #+END_SRC

***** EIG-LCMV

      #+NAME: ./tg_s05_Spatial_filtering___02_Definitions.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___02_Definitions.m :results silent

rec_flt.EIG_LCMV_R = rec_flt.LCMV_R*P_EIG;
rec_flt.EIG_LCMV_N = rec_flt.LCMV_N*P_EIG;

      #+END_SRC

***** PREWHITENED-LCMV

      #+NAME: ./tg_s05_Spatial_filtering___02_Definitions.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___02_Definitions.m :results silent

% not applicable in our case as it assumes access to covariance matrix of interference + background noise + measurement noise, whereas in our case N does not contain interference

      #+END_SRC

***** sMVP

      #+NAME: ./tg_s05_Spatial_filtering___02_Definitions.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___02_Definitions.m :results silent

rec_flt.sMVP_MSE = P_sMVP_MSE_opt*rec_flt.LCMV_R;
rec_flt.sMVP_R = P_sMVP_R_opt*rec_flt.LCMV_R;
rec_flt.sMVP_N = P_sMVP_N_opt*rec_flt.LCMV_N;

      #+END_SRC

***** sMVP_NL

      #+NAME: ./tg_s05_Spatial_filtering___02_Definitions.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___02_Definitions.m :results silent

rec_flt.sMVP_NL_MSE = P_sMVP_NL_MSE_opt*pinv(P_sMVP_NL_R*H_Src_R)*P_sMVP_NL_R*pinv(sqrtm(R));
rec_flt.sMVP_NL_R = P_sMVP_NL_R_opt*pinv(P_sMVP_NL_R*H_Src_R)*P_sMVP_NL_R*pinv(sqrtm(R));
rec_flt.sMVP_NL_N = P_sMVP_NL_N_opt*pinv(P_sMVP_NL_N*H_Src_N)*P_sMVP_NL_N*pinv(sqrtm(N));

      #+END_SRC

***** Optional clearing of selected filters

      #+NAME: ./tg_s05_Spatial_filtering___02_Definitions.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___02_Definitions.m :results silent

if(SETUP.fltREMOVE)
    rec_flt = rmfield(rec_flt,'RANDN');
    rec_flt = rmfield(rec_flt,'ZEROS');
end

      #+END_SRC

***** Cleanup

      #+NAME: ./tg_s05_Spatial_filtering___02_Definitions.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___02_Definitions.m :results silent

rawFixStrJoin

clearvars R N C C_NL S G P_* H_* jj;

      #+END_SRC

**** !opt Checkups

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos
rec_flt
rec_opt.ranks

     #+END_SRC

*** DONE Execution
**** !run Signal reconstruction

     Reconstruct signal in the sources using previously prepared
     spatial filters.

     #+NAME: ./tg_s05_Spatial_filtering___03_Execution.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___03_Execution.m :results silent

clearvars rec_sig ii jj kk nn tmp*
rec_sig.Original = sim_sig_SrcActiv.sigSRC_pst;
rec_sig.Dummy    = sim_sig_SrcActiv.sigSRC_pre;
tmp_fltFields    = fieldnames(rec_flt);

for nn = 1:length(tmp_fltFields),
    for kk = 1:SETUP.K00,
        rec_sig.(tmp_fltFields{nn})(:,:,kk) = (rec_flt.(tmp_fltFields{nn}) * y_Pst(:,:,kk)')';
    end,
end

clearvars ii jj kk nn tmp*

     #+END_SRC

***** !opt Checkups

      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos
rec_flt
rec_opt
rec_sig
whos

      #+END_SRC

**** !run Functional dependencies

     Calculate MVAR model coefficients and PDC.

     NB: Keep in mind that P00 is restricted to single value, i.e. the
     actual order of the MVAR model used to generate time-courses for
     signal of interest (it is not being estimated here).

     #+NAME: ./tg_s05_Spatial_filtering___03_Execution.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Spatial_filtering___03_Execution.m :results silent

clearvars ii jj kk nn rec_funDep_A00 rec_funDep_PDC
rec_funDep_A00.Original = sim_sig_SrcActiv.A00;

[~,rec_funDep_A00.Dummy,~,~,~,~] = arfit(sim_sig_SrcActiv.sigSRC_pre,sim_sig_SrcActiv.P00,sim_sig_SrcActiv.P00);

tmp_fltFields    = fieldnames(rec_flt);
for nn = 1:length(tmp_fltFields),[~,rec_funDep_A00.(tmp_fltFields{nn}),~,~,~,~] = arfit(rec_sig.(tmp_fltFields{nn}),SETUP.P00,SETUP.P00);end

tmp_allFields    = fieldnames(rec_funDep_A00);
for nn = 1:length(tmp_allFields),rec_funDep_PDC.(tmp_allFields{nn}) = abs(PDC(rec_funDep_A00.(tmp_allFields{nn}),SETUP.PDC_RES));end

     #+END_SRC

***** !opt Checkups

      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

whos
whos rec*
rec_flt
rec_opt
rec_opt.ranks
rec_sig
rec_funDep_A00
rec_funDep_PDC

      #+END_SRC

** DONE [s05] Localization
*** DONE Locallizers preparations
**** !run Preparations
     #+NAME: ./tg_s05_Localization___00_Preparations.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Localization___00_Preparations.m :results silent

       % pre-computed leadfields for source candidates
       load('./mat/sel_src.mat');    

     #+END_SRC
*** DONE Locallizers definitions 
**** !run MAI and MPZ Locallizers
     #+NAME: ./tg_s05_Localization___01_MAI_and_MPZ_Locallizers.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s05_Localization___01_MAI_and_MPZ_Locallizers.m :results silent

       structure = sel_src.lfg;
       % lf is the array of cells containing all leafields within ROI
       lf = structure.leadfield; 
       % position of sources
       pos_Src = sim_lfg_SrcActiv_orig.lfg.pos;
       [~, pos_Src_ind]=ismember(pos_Src,structure.pos,'rows');

       % Results will be stored in structs with fields: 
       % leadfield, source positions, localization error, rank selected
       % 1 = MAI(2011), 2 = MAI(2013), 3 = MPZ(2011), 4 = MPZ(2013)
       % 5 = MAI_RR_I, 6 = MAI_RR_C, 7 = MPZ_RR_I, 8 = MPZ_RR_C
       RES = repmat(struct('H',[],'sources',[]),8,1);

       % to be used for the proposed heuristic rank-selection criterion 
       % leading eigenvalues of R*pinv(N) are equal to those of C*G0+1
       x = sort(eig(R*pinv(N)),'descend');
       x = x(1:sum(SETUP.SRCS(:,1)));
       % percentage of variance explained
       y = x./sum(x);
       % percent of explained variance
       rs_opt = min(find(cumsum(y)>0.8));

       clear x y;      

       for i = 1:sum(SETUP.SRCS(:,1))
           tic
           i
           % for REAL WORLD APPLICABLE activity indices, rank is fixed
           RES_tmp = repmat({zeros(length(lf),1)}, 1, 8);
           rs = min([i rs_opt]);
           
           for k = 1:8
               
               for j = 1:length(lf)
                   
                   h_curr = lf{j};
                   
                   % we use H obtained for a given index at previous iteration with respect to i, 
                   % starting with null matrix
                   h = [RES(k).H h_curr];

                   G = h'*pinv(N)*h;
                   S = h'*pinv(R)*h;
                   T = h'*pinv(R)*N*pinv(R)*h;
                   
                   [U D ~] = svd(S);            
                   
                   MAI_eigenv = sort(real(eig(G*pinv(S))),'descend');
                   MPZ_eigenv = sort(real(eig(S*pinv(T))),'descend');

                   if ismember(k,[6 8]) & i>rs_opt
                       % Csq is rank-deficient if h = lf{j} matches previously found sources
                       Csq = real(sqrtm(pinv(S)-pinv(G)));            
                       h2 = h*Csq;
                       
                       G2 = h2'*pinv(N)*h2;
                       S2 = h2'*pinv(R)*h2;
                       T2 = h2'*pinv(R)*N*pinv(R)*h2;
                       
                       [U2 D2 ~] = svd(S2);
                   end
                   
                   % REAL-WORLD APPLICABLE family of indices
                   switch k
                     case 1
                       % MAI(2011) index
                       RES_tmp{k}(j) = sum(MAI_eigenv);
                     case 2             
                       % MAI(2013) index
                       RES_tmp{k}(j) = sum(MAI_eigenv(1:rs));
                     case 3
                       % MPZ(2011) index
                       RES_tmp{k}(j) = sum(MPZ_eigenv);    
                     case 4
                       % MPZ(2013) index
                       RES_tmp{k}(j) = sum(MPZ_eigenv(1:rs));
                     case 5
                       % MAI_RR_I index
                       P = U(:,1:rs)*U(:,1:rs)';
                       RES_tmp{k}(j) = trace(G*pinv(S)*P);
                       clear P;
                     case 6
                       % MAI_RR_C index
                       if i<=rs_opt
                           RES_tmp{k}(j) = RES_tmp{k-1}(j);
                       elseif i>rs_opt % in this case rs=rs_opt    
                           P2 = U2(:,1:rs_opt)*U2(:,1:rs_opt)';
                           RES_tmp{k}(j) = trace(G2*pinv(S2)*P2);
                           clear P2;
                       end
                     case 7
                       % MPZ_RR_I index
                       P = U(:,1:rs)*U(:,1:rs)';
                       RES_tmp{k}(j) = trace(S*pinv(T)*P);
                       clear P;                       
                     case 8
                       % MPZ_RR_C index
                       if i<=rs_opt
                           RES_tmp{k}(j) = RES_tmp{k-1}(j);
                       elseif i>rs_opt % in this case rs=rs_opt        
                           P2 = U2(:,1:rs_opt)*U2(:,1:rs_opt)';
                           RES_tmp{k}(j) = trace(S2*pinv(T2)*P2);
                           clear P2;
                       end
                   end
                   
                   % explicit clear to emphasize that all matrices 
                   % are created anew at each iteration for each j and k
                   clear h_curr h G S T Csq h2 G2 S2 T2 U D U2 D2 P P2;
                   
                   % below: end of j = 1:length(lf) loop    
               end

               [~, w]=max(RES_tmp{k});
               [dis dis_ind] = pdist2(pos_Src,structure.pos(w,:),'chebychev','Smallest',1);

               % store results
               RES(k).H = [RES(k).H lf{w}];
               RES(k).sources(i).candidate_number = w;
               RES(k).sources(i).rank_opt = rs;
               RES(k).sources(i).source_discovered = dis_ind;
               RES(k).sources(i).distance = dis;

               % below: end of k = 1:8 loop
           end
           
           % below: end of i = 1:SETUP.RANK_EIG loop            
           toc
       end

       clear structure lf pos_Src;
     #+END_SRC
** DONE [s06] Reconstruction error evaluation
*** !opt Code for testing of source reconstruction error evaluation

    Replace original signal with sine and the reconstructed signals
    with waves that produce easy to predict error values...

    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

if 0
    clearvars rec* ii jj kk nn tmp*

    tmp_arg = linspace(0,8*pi,SETUP.n00)';
    size(tmp_arg)

    rec_sig.Original = repmat(sin(tmp_arg),[1,sum(SETUP.SRCS(:,1)),SETUP.K00]);
    rec_sig.Dummy    = repmat(cos(tmp_arg),[1,sum(SETUP.SRCS(:,1)),SETUP.K00]);

    rec_sig.Test00   = rec_sig.Original
    rec_sig.Test01   = -rec_sig.Original
    rec_sig.Test02   = rec_sig.Original+0.2*randn(size(rec_sig.Original))
    rec_sig.Test03   = -rec_sig.Original+0.2*randn(size(rec_sig.Original))
    rec_sig.Test04   = zeros(size(rec_sig.Original))
end

clearvars ii jj kk nn tmp*
whos
whos rec*
rec_sig

    #+END_SRC
*** Preparation
    
    #+NAME: ./tg_s06_Error_evaluation.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s06_Error_evaluation.m :results silent

      clearvars ii jj kk nn tmp* 
      tmp_allFields  = fieldnames(rec_sig);

      for nn = 1:length(tmp_allFields)
          for ii = 1:SETUP.n00
              for jj = 1:SETUP.K00
 rec_sig.(tmp_allFields{nn})(ii,:,jj)=rec_sig.(tmp_allFields{nn})(ii,:,jj)/norm(rec_sig.(tmp_allFields{nn})(ii,:,jj));
              end
          end
      end

    #+END_SRC

*** =Euclid= for signal

    #+NAME: ./tg_s06_Error_evaluation.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s06_Error_evaluation.m :results silent

      clearvars ii jj kk nn tmp* rec_sigAmp_ErrEuclid
      tmp_allFields  = fieldnames(rec_sig);

      for nn = 1:length(tmp_allFields)
          tmp_error = zeros(SETUP.n00,SETUP.K00);
          for ii = 1:SETUP.n00
              for jj = 1:SETUP.K00
                  tmp_error(ii,jj)=norm(rec_sig.(tmp_allFields{1})(ii,:,jj)-rec_sig.(tmp_allFields{nn})(ii,:,jj))^2;                  
              end
          end
          rec_sigAmp_ErrEuclid.(tmp_allFields{nn})=mean(mean(tmp_error));
      end
    #+END_SRC

**** !opt Checkups

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :results silent

tmp_allFields  = fieldnames(rec_sig);
disp(tmp_allFields{1})

isequal(rec_sig.(tmp_allFields{1}),sim_sig_SrcActiv.sigSRC_pst)
norm(rec_sig.(tmp_allFields{1})(:)-sim_sig_SrcActiv.sigSRC_pst(:))

size(rec_sig.(tmp_allFields{1}))
size(sim_sig_SrcActiv.sigSRC_pre)
size(sim_sig_AdjSNRs.SrcActivPre)

     #+END_SRC

*** =CorrCf= for signal

    #+NAME: ./tg_s06_Error_evaluation.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s06_Error_evaluation.m :results silent

      clearvars ii jj kk nn tmp* rec_sigAmp_ErrCorrCf
      tmp_allFields  = fieldnames(rec_sig);

      for nn = 1:length(tmp_allFields)
          tmp_error = zeros(sum(SETUP.SRCS(:,1)),SETUP.K00);
          for kk = 1:sum(SETUP.SRCS(:,1))
              for jj = 1:SETUP.K00
                  tmp_x=corrcoef(rec_sig.(tmp_allFields{1})(:,kk,jj),rec_sig.(tmp_allFields{nn})(:,kk,jj));
                  tmp_error(kk,jj)=tmp_x(1,end);
              end
          end
          rec_sigAmp_ErrCorrCf.(tmp_allFields{nn}) = mean(mean(tmp_error));
      end

    #+END_SRC

*** =Euclid= for MVAR coefficients matrix

    #+NAME: ./tg_s06_Error_evaluation.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s06_Error_evaluation.m :results silent

  
      clearvars ii jj kk nn tmp* rec_funDep_A00_ErrMonday rec_funDep_A00_ErrEuclid
      tmp_allFields  = fieldnames(rec_sig);

      for nn = 1:length(tmp_allFields)
          rec_funDep_A00_ErrEuclid.(tmp_allFields{nn}) = (norm(rec_funDep_A00.(tmp_allFields{1}) - rec_funDep_A00.(tmp_allFields{nn}),'fro')^2);
      end

    #+END_SRC

*** =CorrCf= for MVAR coefficient matrix

    #+NAME: ./tg_s06_Error_evaluation.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s06_Error_evaluation.m :results silent

      clearvars ii jj kk nn tmp* rec_funDep_A00_ErrCorrCf
      tmp_allFields  = fieldnames(rec_sig);

      for nn = 1:length(tmp_allFields)
          tmp_error = zeros(sum(SETUP.SRCS(:,1)),1);
          for kk = 1:sum(SETUP.SRCS(:,1))
              tmp_x=corrcoef(rec_funDep_A00.(tmp_allFields{1})(kk,:),rec_funDep_A00.(tmp_allFields{nn})(kk,:));
              tmp_error(kk)=tmp_x(1,end);
          end
          rec_funDep_A00_ErrCorrCf.(tmp_allFields{nn}) = mean(tmp_error);
      end

    #+END_SRC

*** =Euclid= for PDC coefficient matrix

    #+NAME: ./tg_s06_Error_evaluation.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s06_Error_evaluation.m :results silent

      clearvars ii jj kk nn tmp* rec_funDep_PDC_err rec_funDep_PDC_ErrEuclid
      tmp_allFields  = fieldnames(rec_sig);

      for nn = 1:length(tmp_allFields)
          tmp_error = zeros(sum(SETUP.SRCS(:,1)),1);
          for ii = 1:sum(SETUP.SRCS(:,1))
              tmp_error(ii) = norm(squeeze(rec_funDep_PDC.(tmp_allFields{1})(ii,:,:)) - squeeze(rec_funDep_PDC.(tmp_allFields{nn})(ii,:,:)),'fro')^2;
          end
          rec_funDep_PDC_ErrEuclid.(tmp_allFields{nn}) = mean(tmp_error);
      end

    #+END_SRC
    
*** =CorrCf= for PDC coefficient matrix

    #+NAME: ./tg_s06_Error_evaluation.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s06_Error_evaluation.m :results silent

      clearvars ii jj kk nn tmp*  rec_funDep_PDC_ErrCorrCf
      tmp_allFields  = fieldnames(rec_sig);

      for nn = 1:length(tmp_allFields)
          tmp_error = zeros(sum(SETUP.SRCS(:,1)),1);
          for ii = 1:sum(SETUP.SRCS(:,1))
              tmp_x=corrcoef(reshape(squeeze(rec_funDep_PDC.(tmp_allFields{1})(ii,:,:)),[],1),reshape(squeeze(rec_funDep_PDC.(tmp_allFields{nn})(ii,:,:)),[],1));
              tmp_error(ii)=tmp_x(1,end);
          end
          rec_funDep_PDC_ErrCorrCf.(tmp_allFields{nn}) = mean(tmp_error);
      end

    #+END_SRC

*** Vectorize results

    #+NAME: ./tg_s06_Error_evaluation.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s06_Error_evaluation.m :results silent

clearvars ii jj kk nn tmp* rec*vec

rec_sigAmp_ErrEuclid_vec = cellfun(@(fn) rec_sigAmp_ErrEuclid.(fn), fieldnames(rec_sigAmp_ErrEuclid), 'UniformOutput', false);
rec_sigAmp_ErrEuclid_vec = squeeze(vertcat(rec_sigAmp_ErrEuclid_vec{:}));

rec_sigAmp_ErrCorrCf_vec = cellfun(@(fn) rec_sigAmp_ErrCorrCf.(fn), fieldnames(rec_sigAmp_ErrCorrCf), 'UniformOutput', false);
rec_sigAmp_ErrCorrCf_vec = squeeze(vertcat(rec_sigAmp_ErrCorrCf_vec{:}));

rec_funDep_A00_ErrEuclid_vec = cellfun(@(fn) rec_funDep_A00_ErrEuclid.(fn), fieldnames(rec_funDep_A00_ErrEuclid), 'UniformOutput', false);
rec_funDep_A00_ErrEuclid_vec = squeeze(vertcat(rec_funDep_A00_ErrEuclid_vec{:}));

rec_funDep_A00_ErrCorrCf_vec = cellfun(@(fn) rec_funDep_A00_ErrCorrCf.(fn), fieldnames(rec_funDep_A00_ErrCorrCf), 'UniformOutput', false);
rec_funDep_A00_ErrCorrCf_vec = squeeze(vertcat(rec_funDep_A00_ErrCorrCf_vec{:}));

rec_funDep_PDC_ErrEuclid_vec = cellfun(@(fn) rec_funDep_PDC_ErrEuclid.(fn), fieldnames(rec_funDep_PDC_ErrEuclid), 'UniformOutput', false);
rec_funDep_PDC_ErrEuclid_vec = squeeze(vertcat(rec_funDep_PDC_ErrEuclid_vec{:}));

rec_funDep_PDC_ErrCorrCf_vec = cellfun(@(fn) rec_funDep_PDC_ErrCorrCf.(fn), fieldnames(rec_funDep_PDC_ErrCorrCf), 'UniformOutput', false);
rec_funDep_PDC_ErrCorrCf_vec = squeeze(vertcat(rec_funDep_PDC_ErrCorrCf_vec{:}));

    #+END_SRC

*** Combine results in single array

    #+NAME: ./tg_s06_Error_evaluation.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s06_Error_evaluation.m :results silent

clearvars ii jj kk nn tmp* rec_res

if     ~isequal(fieldnames(rec_sigAmp_ErrEuclid),fieldnames(rec_sigAmp_ErrCorrCf))
    error('check fieldnames')
elseif ~isequal(fieldnames(rec_sigAmp_ErrEuclid),fieldnames(rec_funDep_A00_ErrEuclid))
    error('check fieldnames')
elseif ~isequal(fieldnames(rec_sigAmp_ErrEuclid),fieldnames(rec_funDep_A00_ErrCorrCf))
    error('check fieldnames')
elseif ~isequal(fieldnames(rec_sigAmp_ErrEuclid),fieldnames(rec_funDep_PDC_ErrEuclid))
    error('check fieldnames')
elseif ~isequal(fieldnames(rec_sigAmp_ErrEuclid),fieldnames(rec_funDep_PDC_ErrCorrCf))
    error('check fieldnames')
else
    rec_res.table_arrC = [ rec_sigAmp_ErrEuclid_vec,  rec_sigAmp_ErrCorrCf_vec,  rec_funDep_A00_ErrEuclid_vec,  rec_funDep_A00_ErrCorrCf_vec,  rec_funDep_PDC_ErrEuclid_vec,  rec_funDep_PDC_ErrCorrCf_vec];
    rec_res.table_varN = {'rec_sigAmp_ErrEuclid_vec','rec_sigAmp_ErrCorrCf_vec','rec_funDep_A00_ErrEuclid_vec','rec_funDep_A00_ErrCorrCf_vec','rec_funDep_PDC_ErrEuclid_vec','rec_funDep_PDC_ErrCorrCf_vec'};
    rec_res.table_rowN = fieldnames(rec_sigAmp_ErrEuclid);
    rec_res.table = array2table(rec_res.table_arrC,'VariableNames',rec_res.table_varN,'RowNames',rec_res.table_rowN);
end

    #+END_SRC

**** Checkups

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :results silent

rec_res.table_arrC = [ rec_sigAmp_ErrEuclid_vec,  rec_sigAmp_ErrCorrCf_vec ];
rec_res.table_varN = {'rec_sigAmp_ErrEuclid_vec','rec_sigAmp_ErrCorrCf_vec'};
rec_res.table_rowN = fieldnames(rec_sigAmp_ErrEuclid);
rec_res.table = array2table(rec_res.table_arrC,'VariableNames',rec_res.table_varN,'RowNames',rec_res.table_rowN);

     #+END_SRC

*** !opt Checkups (RESULTS CAN BE CHECKED HERE)

    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :results silent

whos rec*
whos rec*vec

rec_res

sortrows(rec_res.table,'rec_sigAmp_ErrEuclid_vec')
sortrows(rec_res.table,'rec_sigAmp_ErrCorrCf_vec','descend')

sortrows(rec_res.table,'rec_funDep_A00_ErrEuclid_vec')
sortrows(rec_res.table,'rec_funDep_A00_ErrCorrCf_vec','descend')
sortrows(rec_res.table,'rec_funDep_PDC_ErrEuclid_vec')
sortrows(rec_res.table,'rec_funDep_PDC_ErrCorrCf_vec','descend')

    #+END_SRC

*** !opt Checkups

    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :results silent

whos rec*
whos rec*vec

rec_sigAmp_ErrEuclid
rec_sigAmp_ErrCorrCf
rec_funDep_A00_ErrEuclid
rec_funDep_A00_ErrMonday
rec_funDep_PDC_ErrEuclid
rec_funDep_PDC_ErrMonday

    #+END_SRC

*** !opt Checkups

    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

rec_sig
rec_sigAmp_ErrEuclid
[r_opt2,r_opt3,size(vSRC_pst,2)]




TMP_COV = cov(sim_sig_SrcActiv.sigSRC_pst(:,:,1));
TMP_COR = corrcoef(sim_sig_SrcActiv.sigSRC_pst(:,:,1));

figure(1);clf;
rawImgSC(TMP_COV);
tmp_colorMapMat = rawHotColdColorMap(2048);
caxis([-max(abs(min(min(TMP_COV))),abs(max(max(TMP_COV)))) max(abs(min(min(TMP_COV))),abs(max(max(TMP_COV))))]);
colormap(tmp_colorMapMat);
colorbar;

figure(2);clf;
rawImgSC(TMP_COR);
tmp_colorMapMat = rawHotColdColorMap(2048);
caxis([-max(abs(min(min(TMP_COR))),abs(max(max(TMP_COR)))) max(abs(min(min(TMP_COR))),abs(max(max(TMP_COR))))]);
colormap(tmp_colorMapMat);
colorbar;





sort(eig(K))
figure(501);clf;set(gcf, 'Position', SETUP.DISP);rawImgSC(sort(eig(K)));colorbar;
rank(Hp)
svd(Hp)
svd(S1)
R



rec_sig
rec_sigAmp_ErrEuclid
[r_opt2,r_opt3,size(vSRC_pst,2)]
rec_funDep_A00
rec_funDep_A00_ErrEuclid
rec_funDep_PDC
rec_funDep_PDC_ErrEuclid
rec_funDep_PDC_ErrMonday



figure(100000)
set(gcf, 'Position', SETUP.DISP);
subplot(2,1,1)
plot(rec_funDep_PDC_ErrEuclid);colorbar;
subplot(2,1,2)
plot(rec_funDep_PDC_ErrMonday);colorbar;



figure(30000);clf;
set(gcf, 'Position', SETUP.DISP);
for nn = 1:length(tmp_allFields)
    subplot(length(tmp_allFields),1,nn)
    imagesc(rec_funDep_A00.(tmp_allFields{nn}));
    tmp_colorMapMat = jet(255);
    tmp_colorMapMat(128,:) = [1 1 1];
    caxis([-max(abs(min(min(rec_funDep_A00.(tmp_allFields{nn})))),abs(max(max(rec_funDep_A00.(tmp_allFields{nn}))))) max(abs(min(min(rec_funDep_A00.(tmp_allFields{nn})))),abs(max(max(rec_funDep_A00.(tmp_allFields{nn})))))]);
    colormap(tmp_colorMapMat);
    colorbar;
    hold on;
    tmp_stem = 0.5+sim_sig_SrcActiv.S00:sim_sig_SrcActiv.S00:sim_sig_SrcActiv.S00*(sim_sig_SrcActiv.P00-1)+0.5;
    tmp_vals = sim_sig_SrcActiv.S00+0.5*ones(size(tmp_stem));
    stem(tmp_stem,tmp_vals,'Color','k','LineWidth',0.5,'Marker', 'none');
    set(gca,'XTick',[1:sim_sig_SrcActiv.S00*sim_sig_SrcActiv.P00])
    title(tmp_allFields{nn})
end



rec_sigAmp_ErrEuclid
[r_opt2,r_opt3,size(vSRC_pst,2)]
rec_funDep_PDC_ErrEuclid
rec_funDep_PDC_ErrMonday



for nn = 1:length(tmp_allFields)
    figure(40000+nn);clf;
    set(gcf, 'Position', SETUP.DISP);
    rawPlotPDC(rec_funDep_PDC.(tmp_allFields{nn}),SETUP.PDC_RES)
    title(tmp_allFields{nn})
end

    #+END_SRC

** DONE [s07] BATCH
    
   Single run for tangled matlab-scripts.

*** [p00] RUN (this is the one for single simulations run)

    The below code block can be used to test simulations single run.

    *NB: before running the code below remember to tangle code blocks
    after any modifications !!!*
    
    #+NAME: ./tg_s07_BATCH___00_RUN.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s07_BATCH___00_RUN.m :results silent

      % if exist('fPath'), cd(fPath); else, try, cd('~/supFunSim/'); catch, warningMessage = 'Problem encoutered while trying to change working directory to ''~/supFunSim/''.'; end; end; disp(['CYBERCRAFT:: pwd is: ',pwd])

      run('./tg_s07_BATCH___01_Prelude.m'); whos
      SETUP
      chkSim___tg_s01_PRE___001(SETUP);
      chkSim___tg_s01_PRE___000_PARSE_SETUP(SETUP,sel_atl);
      run('./tg_s07_BATCH___02_Signals.m'); whos
      run('./tg_s07_BATCH___03_Leadfields.m'); whos

      % if BcgNoise (biological background  noise) is to be white Gaussian,
      % then please uncomment the following two lines. 
      % sim_sig_BcgNoise.sigSRC_pre = randn(size(sim_sig_BcgNoise.sigSRC_pre));
      % sim_sig_BcgNoise.sigSRC_pst = randn(size(sim_sig_BcgNoise.sigSRC_pst));
      % if one wishes the biological noise on sensors to be white Gaussian
      % sim_sig_BcgNoise.sigSNS_pre = randn(size(sim_sig_BcgNoise.sigSNS_pre));
      % sim_sig_BcgNoise.sigSNS_pst = randn(size(sim_sig_BcgNoise.sigSNS_pst));

      % rawImgSC(squareform(pdist(sim_lfg_SrcActiv_orig.lfg.pos,'euclidean')))

      run('./tg_s07_BATCH___04_Preparations.m');    whos

      switch SETUP.supSwitch
        case {'rec'}

          % Reconstruction
          disp('Reconstruction')
          run('./tg_s07_BATCH___05_Filters.m'); whos

          % Check the results
          sortrows(rec_res.table,'rec_sigAmp_ErrEuclid_vec')
          sortrows(rec_res.table,'rec_sigAmp_ErrCorrCf_vec','descend')
          sortrows(rec_res.table,'rec_funDep_A00_ErrEuclid_vec')
          sortrows(rec_res.table,'rec_funDep_A00_ErrCorrCf_vec','descend')
          sortrows(rec_res.table,'rec_funDep_PDC_ErrEuclid_vec')
          sortrows(rec_res.table,'rec_funDep_PDC_ErrCorrCf_vec','descend')
          rec_opt.ranks

        case {'loc'}

          % Localization
          disp('Localization')
          run('./tg_s07_BATCH___05_Localizers.m'); whos

        otherwise
          warning('Unexpected supSwitch value. No action taken, only data generated.')
      end

    #+END_SRC
    
**** !opt Checkups

     After execution of =tg_s07_BATCH___02_Signals.m= the generated
     timeseries can be swapped for simple waveforms (used for to
     facilitate *bug tracking* and for *explanatory/illustrative
     purpose*).

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

% (just trying simple signals over here)

tmp_arg = linspace(0,8*pi,SETUP.n00)';
size(tmp_arg)

sim_sig_SrcActiv.sigSRC_pre = repmat(sin(tmp_arg),    [1,sum(SETUP.SRCS(:,1)),SETUP.K00]);
sim_sig_SrcActiv.sigSRC_pst = repmat(sin(tmp_arg),    [1,sum(SETUP.SRCS(:,1)),SETUP.K00]);
sim_sig_IntNoise.sigSRC_pre = repmat(sin(tmp_arg*100),[1,sum(SETUP.SRCS(:,2)),SETUP.K00]);
sim_sig_IntNoise.sigSRC_pst = repmat(sin(tmp_arg*100),[1,sum(SETUP.SRCS(:,2)),SETUP.K00]);
sim_sig_BcgNoise.sigSRC_pre = repmat(sin(tmp_arg*10), [1,sum(SETUP.SRCS(:,3)),SETUP.K00]);
sim_sig_BcgNoise.sigSRC_pst = repmat(sin(tmp_arg*10), [1,sum(SETUP.SRCS(:,3)),SETUP.K00]);

size(sim_sig_SrcActiv.sigSRC_pre)
size(sim_sig_SrcActiv.sigSRC_pst)
size(sim_sig_BcgNoise.sigSRC_pre)
size(sim_sig_BcgNoise.sigSRC_pst)

close all
jj = 1
kk = 1
figure('Name','sim_sig_SrcActiv.sigSRC_pre'); plot(sim_sig_SrcActiv.sigSRC_pre(:,jj,kk))
figure('Name','sim_sig_BcgNoise.sigSRC_pre'); plot(sim_sig_BcgNoise.sigSRC_pre(:,jj,kk))

figure('Name','sim_sig_SrcActiv.sigSNS_pre'); plot(sim_sig_SrcActiv.sigSNS_pre(:,jj,kk))
figure('Name','sim_sig_SrcActiv.sigSNS_pst'); plot(sim_sig_SrcActiv.sigSNS_pst(:,jj,kk))

figure('Name','sim_sig_BcgNoise.sigSRC_pre'); plot(sim_sig_BcgNoise.sigSRC_pre(:,jj,kk))
figure('Name','sim_sig_BcgNoise.sigSRC_pst'); plot(sim_sig_BcgNoise.sigSRC_pst(:,jj,kk))

figure('Name','sim_sig_BcgNoise.sigSNS_pre'); plot(sim_sig_BcgNoise.sigSNS_pre(:,jj,kk))
figure('Name','sim_sig_BcgNoise.sigSNS_pst'); plot(sim_sig_BcgNoise.sigSNS_pst(:,jj,kk))

size(y_Pre)
size(y_Pst)

figure('Name','y_Pre'); plot(y_Pre(:,jj,kk))
figure('Name','y_Pst'); plot(y_Pst(:,jj,kk))

     #+END_SRC

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

% (just trying simple signals over here)

sim_sig_BcgNoise.sigSRC_pre = randn(size(sim_sig_BcgNoise.sigSRC_pre))
sim_sig_BcgNoise.sigSRC_pst = sim_sig_BcgNoise.sigSRC_pre

     #+END_SRC

     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle no :results silent

% (just trying simple signals over here)

sim_sig_BcgNoise.sigSRC_pst = sim_sig_BcgNoise.sigSRC_pre

     #+END_SRC

*** [p01] Prelude (with *optional settings overwrite*)
    :PROPERTIES:
    :ID:       6718d735-94d9-48ae-ae42-1e9f3b3e0799
    :END:

    The following settings will overwritte simulation default settings.
    The devault settings are optimized for *bug tracking* and for
    *explanatory/illustrative purpose*.  Please modify simulation
    settings only in this section.

    #+NAME: ./tg_s07_BATCH___01_Prelude.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s07_BATCH___01_Prelude.m :results silent :var fPath = (file-name-directory (or load-file-name buffer-file-name))

      % Tidy up and change working directory.
      clc; close all;clearvars('-except','fPath');
      if exist('fPath'),cd(fPath);else,try,cd('~/cc_overkill/git/supFunSim');catch,cd('~/supFunSim');end;end;
      clear all;
      close all;
      clc;

      run('./tg_s00_Prelude___01_Simulations_main_setup.m')
      run('./tg_s00_Prelude___02_SNR_Adjustments.m')

      % Rows of SETUP.SRCS reppresent ROIs.
      % Cols of SETUP.SRCS represent SrcActiv, IntNoise and BcgNoise, respectively.
      SETUP.rROI   = logical(1);       % random (1) or predefined (0) ROIs
      SETUP.rPNT   = logical(1);       % random (1) or predefined (0) candidate points for source locations: if 0, 
				       % number of sources as in SETUP.SRCS(1,1) will be fixed and in close locations
      SETUP.SRCS   = []; % Cortical sources (avoid placing more than 10 sources in single ROI)
      SETUP.SRCS   = [ SETUP.SRCS;  3  3  0 ];
      SETUP.SRCS   = [ SETUP.SRCS;  3  3  0 ];
      SETUP.SRCS   = [ SETUP.SRCS;  3  3  0 ];
      SETUP.SRCS   = [ SETUP.SRCS;  4  3  0 ];
      SETUP.SRCS   = [ SETUP.SRCS;  0  3  0 ];
      SETUP.SRCS   = [ SETUP.SRCS;  0  3  0 ];
      SETUP.SRCS   = [ SETUP.SRCS;  0  3  0 ];
      SETUP.SRCS   = [ SETUP.SRCS;  0  3  0 ];
      SETUP.SRCS   = [ SETUP.SRCS;  0  3  3 ];
      SETUP.SRCS   = [ SETUP.SRCS;  0  0  4 ];
      SETUP.DEEP   = [              0  0  20 ]; % deep sources

      SETUP.ERPs   = 0;       % Add ERPs (timelocked activity)

      % Basic setting s for sources and noise.
      SETUP.n00    = 500;     % number of time samples per trial
      SETUP.K00    = 1;       % number of independent realizations of signal and noise based on generated MVAR model
			      % note: covariance matrix R of observed signal and noise covariance matrix N are
			      % estimated from samples originating from all realizations of signal and noise
      SETUP.FIXED_SEED = 1; % Settings for seed selection
      if SETUP.FIXED_SEED, SETUP.SEED = rng(1964);else,SETUP.SEED = rng(round(1e3*randn()^2*sum(clock)));end
      SETUP.RANK_EIG = sum(SETUP.SRCS(:,1)); % rank of EIG-LCMV filter: set to number of active sources
      SETUP.fltREMOVE = 1;    % to keep (0) or remove (1) selected filters
      SETUP.SHOWori = 0; % to show (1) or do not show (0) Original and Dummy signals on Figures
      SETUP.IntLfgRANK = round(0.3*sum(SETUP.SRCS(:,2))); % rank of patch-constrained reduced-rank leadfield
      SETUP.supSwitch = 'rec'; % 'rec': run reconstruction of sources activity, 'loc': find active sources

      SETUP.CUBE           = 20;    % perturbation of the leadfields based on the shift of source position within a cube of given edge length (centered at the original leadfields positions), in [mm]
      SETUP.CONE           = pi/32; % perturbation of the leadfields based on the rotation of source orientation (azimuth TH, elevation PHI), in [rad]
      SETUP.H_Src_pert     = 1;     % use original (0) or perturbed (1) leadfield for signal reconstruction
      SETUP.H_Int_pert     = 0;     % use original (0) or perturbed (1) leadfield for nulling constrains
      SETUP.SINR           = 0;     % signal to interference noise power ratio expressed in dB (both measured on electrode level)
      SETUP.SBNR           = 0;    % signal to biological noise power ratio expressed in dB (both measured on electrode level)
      SETUP.SMNR           = 0;    % signal to measurment noise power ratio expressed in dB (both measured on electrode level)
      SETUP.WhtNoiseAddFlg = 1;     % white noise admixture in biological noise interference noise (FLAG)
      SETUP.WhtNoiseAddSNR = 3;     % SNR of BcgNoise and WhiNo (dB)
      SETUP.SigPre = 0;   SETUP.IntPre = 0;   SETUP.BcgPre = 1;   SETUP.MesPre = 1; % final signal components for pre-interval  (use zero or one for signal, interference noise, biological noise, measurement noise)
      SETUP.SigPst = 1;   SETUP.IntPst = 1;   SETUP.BcgPst = 1;   SETUP.MesPst = 1; % final signal components for post-interval (as above)

      % For localization, the default is to consider random locations of ROI with some fixed candidate points such that 
      % SETUP.SRCS(1,1) of them are in close locations; additionally, no interfering sources are considered
      % in the localization model. You may comment out the 'if' section below to experiment with other settings.
      if(SETUP.supSwitch == 'loc') 
	  SETUP.rROI   = logical(1);       % random (1) or predefined (0) ROIs
	  SETUP.rPNT   = logical(0);       % random (1) or predefined (0) candidate points for source locations: if 0, 
					   % number of sources as in SETUP.SRCS(1,1) will be fixed and in close locations
	  SETUP.SigPre = 0;   SETUP.IntPre = 0;   SETUP.BcgPre = 1;   SETUP.MesPre = 1; % final signal components for pre-interval  (use zero or one for signal, interference noise, biological noise, measurement noise)
	  SETUP.SigPst = 1;   SETUP.IntPst = 0;   SETUP.BcgPst = 1;   SETUP.MesPst = 1; % final signal components for post-interval (as above)
      end

      if SETUP.rPNT
	  disp('CC: using random source locations')
      else
	  disp('CC: using predefined source locations')
      end

    #+END_SRC

*** [p02] Signals

    #+NAME: ./tg_s07_BATCH___02_Signals.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s07_BATCH___02_Signals.m :results silent

disp('CYBERCRAFT:: Generation of timeseries for bioelectrical activity of interest')
run('./tg_s01_Timeseries___01_SrcActiv.m')

disp('CYBERCRAFT:: Generation of timeseries for bioelectrical interference noise')
run('./tg_s01_Timeseries___02_IntNoise.m')

disp('CYBERCRAFT:: Generation of timeseries for bioelectrical background noise')
run('./tg_s01_Timeseries___03_BcgNoise.m')

disp('CYBERCRAFT:: Measurement noise generation')
run('./tg_s01_Timeseries___04_MesNoise.m')

    #+END_SRC

*** [p03] Leadfields

    #+NAME: ./tg_s07_BATCH___03_Leadfields.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s07_BATCH___03_Leadfields.m :results silent

run('./tg_s02_Geometry___01_RandSamp.m')
run('./tg_s02_Geometry___02_Indices.m')
run('./tg_s02_Geometry___03_Coordinates.m')
run('./tg_s02_Geometry___04_DeepSrc.m')
run('./tg_s02_Geometry___05_Perturbation.m')
run('./tg_s02_Geometry___06_lfg.m')
run('./tg_s03_Forward_modeling.m')

    #+END_SRC

*** [p04] Preparations

    #+NAME: ./tg_s07_BATCH___04_Preparations.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s07_BATCH___04_Preparations.m :results silent

run('./tg_s04_Preparation___00_SNRs_Adjustment.m')
run('./tg_s04_Preparation___01_Signal_Covariances.m')


    #+END_SRC
    
*** [p05] Filters

    #+NAME: ./tg_s07_BATCH___05_Filters.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s07_BATCH___05_Filters.m :results silent

run('./tg_s05_Spatial_filtering___01_Preparation.m')
run('./tg_s05_Spatial_filtering___02_Definitions.m')
run('./tg_s05_Spatial_filtering___03_Execution.m')
run('./tg_s06_Error_evaluation.m')

    #+END_SRC

*** [p05] Localizers
    
    #+NAME: ./tg_s07_BATCH___05_Localizers.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s07_BATCH___05_Localizers.m :results silent

      run('./tg_s05_Localization___00_Preparations.m')
      run('./tg_s05_Localization___01_MAI_and_MPZ_Locallizers.m')

    #+END_SRC
    
** DONE [s08] LOOP

   Multiple runs for tangled matlab-scripts.

*** [p00] RUN All
**** Execute loop
     
     - Simulation settings can be modified in:
     - [[*%5Bs00%5D%20Prelude][{s00} Prelude]]
     - [[id:6718d735-94d9-48ae-ae42-1e9f3b3e0799][Prelude (with *optional settings overwrite*)]], above
     - the code block below

     #+NAME: ./tg_s08_LOOP___00_RUN.m
     #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s08_LOOP___00_RUN.m :results silent

    if exist('fPath'), cd(fPath); else, try, cd('~/supFunSim/'); catch, warningMessage = 'Problem encoutered while trying to change working directory to ''~/supFunSim/''.'; end; end; disp(['CYBERCRAFT:: pwd is: ',pwd])

    clearvars ii jj kk nn tmp* LOOP* Loop*;
    run('./tg_s07_BATCH___01_Prelude.m');
    SETUP

    % Get some details for simulation output filename (*.mat file)
    LOOP.DATE = datestr(now,'yyyymmdd_HHMMSS');
    LOOP.NAME = tempname; [~, LOOP.NAME] = fileparts(LOOP.NAME); % simulation unique name

    % Number of simulation runs for each SNRs combination
    LOOP.totSimCount = 1000;
    LOOP.rngSimCount = 1:LOOP.totSimCount;

    % Range of SNRs 
    LOOP.rngMesNoise = 10;
    LOOP.rngBcgNoise = 0;
    LOOP.rngIntNoise = [0:5:10];

    LOOP
    SETUP

    run('./tg_s08_LOOP___01_Main_LOOP.m'); whos

    switch SETUP.supSwitch
      case {'rec'}

	run('./tg_s08_Loop_Reconstruction_Plot___02_.m')

      case {'loc'}

	run('./tg_s08_Loop_Localization_Plot___02_.m')

      otherwise
	warning('Unexpected supSwitch value. No action taken, only data generated.')
    end
     #+END_SRC
**** Extract and visualize errors
***** Reconstruction
      #+NAME: ./tg_s08_Loop_Reconstruction_Plot___02_.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s08_Loop_Reconstruction_Plot___02_.m :results silent

       LOOP
       LOOP.table_varN'
       LOOP.table_rowN
       size(LOOP.table_arrC)

       % avg and std arrays #filters x #error_measures x #SXNR_levels 
       LOOP.table_arrC_avgSimCount = squeeze(mean(LOOP.table_arrC,3))
       LOOP.table_arrC_stdSimCount = squeeze(std(LOOP.table_arrC,[],3))
       size(LOOP.table_arrC_avgSimCount)

       % SAVE THE RESULTS
       mySave = ['./res___',LOOP.DATE,'___',LOOP.NAME,'___',datestr(now,'yyyymmdd_HHMMSS')];
       save(mySave)
       LOOP

       close all
       for ii = 1:length(LOOP.table_varN)
           LOOP.table_varN'
           tmp_pickVar = ii;
           disp(LOOP.table_varN(tmp_pickVar))

           TMP_res.table_arrC = squeeze(LOOP.table_arrC_avgSimCount(:,tmp_pickVar,:))
           TMP_res.table_varN = strrep(strrep(strsplit(num2str(LOOP.(LOOP.NamNoise),'%+d ')),'-','n'),'+','p')
           TMP_res.table_rowN = LOOP.table_rowN
           TMP_res.table = array2table(TMP_res.table_arrC,'VariableNames',TMP_res.table_varN,'RowNames',TMP_res.table_rowN);
           sortrows(TMP_res.table,1)

           figure('Name',[LOOP.table_varN{tmp_pickVar},' for ',LOOP.NamNoise,' in ',num2str(LOOP.(LOOP.NamNoise),'%+d ')])
           bar(TMP_res.table_arrC)
           set(gca, 'XTick', 1:length(TMP_res.table_rowN), 'XTickLabel', strrep(TMP_res.table_rowN,'_','\_'));
           h1 = legend(strsplit(num2str(LOOP.(LOOP.NamNoise),'%+d ')))
           myTitle = strrep([LOOP.table_varN{tmp_pickVar},' for ',LOOP.NamNoise,' in ',num2str(LOOP.(LOOP.NamNoise),'%+d ')],'_','\_');
           % title(myTitle)
           myName = ['./fig/',LOOP.DATE,'___',LOOP.NAME,'___fig_',num2str(ii),'___',strrep(strrep(myTitle,'\',''),' ','_')];

           v1 = get(h1,'title')

           % fix me: it does not have to be SINR
           set(v1,'string','SINR')

           set(gcf, 'PaperUnits', 'centimeters')
           set(gcf, 'PaperOrientation', 'portrait')
           set(gcf, 'PaperSize', [20, 20])
           set(gcf, 'Position', [100, 100, 800, 750])
           legend('Location','northeastoutside')
           xtickangle(45)

           savefig(gcf,myName,'compact')
           savefig(gcf,myName)
           saveas(gcf,myName,'pdf')
       end

      #+END_SRC
***** Localization
      #+NAME: ./tg_s08_Loop_Localization_Plot___02_.m
      #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s08_Loop_Localization_Plot___02_.m :results silent

	LOOP

	% USED FOR SAVING FIGURES
	myName = ['./fig___',LOOP.DATE,'___',LOOP.NAME,'___',datestr(now,'yyyymmdd_HHMMSS')];

	% use LOOP.RES.ITERATION{x}{y}(z) to access results struct for 
	% z-th activity index (z=1,2,...,8 with z=1 for MAI(2011) etc., as defined in !run MAI and MPZ Locallizers) 
	% y-th level of SNR of type considered (interference, biological, measurement)
	% x-th iteration of the experiment

	% XRANK is an array of size |x| by |y|, representing ranks selected for reduced-rank indices
	XRANK=zeros(LOOP.totSimCount,length(LOOP.(LOOP.NamNoise)));

	% number of simulation runs
	for i=1:LOOP.totSimCount
	    % number of levels of SNR
	    for j=1:length(LOOP.(LOOP.NamNoise))
		XRANK(i,j)=max([LOOP.RES.ITERATION{i}{j}(1).sources.rank_opt]);
	    end
	end

	% XRES is an array of size |x| by |y| by |z| by |#sources|, representing errors for each source
	% found using iterative procedure
	XRES=zeros(LOOP.totSimCount,length(LOOP.(LOOP.NamNoise)),8,sum(SETUP.SRCS(:,1)));

	% number of simulation runs
	for i=1:LOOP.totSimCount
	    % number of levels of SNR
	    for j=1:length(LOOP.(LOOP.NamNoise))
		% number of activity indices
		for k=1:8
		    XRES(i,j,k,:)=[LOOP.RES.ITERATION{i}{j}(k).sources.distance]';
		end 
	    end 
	end

	% RRES is same like XRES, with the trailing few sources considered
	RRES=XRES(:,:,:,end-1:end);

	% |x| by |y| by |z| array, with last dimension representing mean of errors across sources discovered
	XRES_s = mean(XRES,4);
	RRES_s = mean(RRES,4);

	% |y| by |z| matrix, with entries representing averaged mean error and its std across experiment repetitions
	XRES_sa = squeeze(mean(XRES_s,1));
	XRES_sstd = squeeze(std(XRES_s,1));
	RRES_sa = squeeze(mean(RRES_s,1));
	RRES_sstd = squeeze(std(RRES_s,1));

	% change order of and remove some activity indices as:
	% 1 = MAI(2011), 2 = MAI(2013), 3 = MAI_RR_I
	% 4 = MPZ(2011), 5 = MPZ(2013), 6 = MPZ_RR_I
	XRES_s = XRES_s(:,:,[1 2 5 3 4 7]);
	RRES_s = RRES_s(:,:,[1 2 5 3 4 7]);

	XRES_sa = XRES_sa(:,[1 2 5 3 4 7]);
	XRES_sstd = XRES_sstd(:,[1 2 5 3 4 7]);

	RRES_sa = RRES_sa(:,[1 2 5 3 4 7]);
	RRES_sstd = RRES_sstd(:,[1 2 5 3 4 7]);

	% Mann-Whitney U-test to check whether proposed activity indices achieve 
	% significantly lower average localization error

	% XRES_U and RRES_U are |y| by 4 (# of comparisons) storing resulting p-values of performed Mann-Whitney U test
	XRES_U=zeros(length(LOOP.(LOOP.NamNoise)),4);
	RRES_U=zeros(length(LOOP.(LOOP.NamNoise)),4);

	for j=1:length(LOOP.(LOOP.NamNoise))
	    XRES_U(j,1)=ranksum(XRES_s(:,j,1),XRES_s(:,j,3),'tail','right');
	    XRES_U(j,2)=ranksum(XRES_s(:,j,2),XRES_s(:,j,3),'tail','right');
	    XRES_U(j,3)=ranksum(XRES_s(:,j,4),XRES_s(:,j,6),'tail','right');
	    XRES_U(j,4)=ranksum(XRES_s(:,j,5),XRES_s(:,j,6),'tail','right');
	end

	for j=1:length(LOOP.(LOOP.NamNoise))
	    RRES_U(j,1)=ranksum(RRES_s(:,j,1),RRES_s(:,j,3),'tail','right');
	    RRES_U(j,2)=ranksum(RRES_s(:,j,2),RRES_s(:,j,3),'tail','right');
	    RRES_U(j,3)=ranksum(RRES_s(:,j,4),RRES_s(:,j,6),'tail','right');
	    RRES_U(j,4)=ranksum(RRES_s(:,j,5),RRES_s(:,j,6),'tail','right');
	end

	% pad with zero vectors to increase spacing between bar groups
	XRES_sa = [XRES_sa; zeros(2,6)];
	XRES_sstd = [XRES_sstd; zeros(2,6)];

	RRES_sa = [RRES_sa; zeros(2,6)];
	RRES_sstd = [RRES_sstd; zeros(2,6)];

	clear gca;      
	figure
	name = {'MAI';'MAI_{ext}';'MAI_{RR-I}';'MPZ';'MPZ_{ext}';'MPZ_{RR-I}'};
	bar(XRES_sa',1.2);
	hold on
	bar(XRES_sstd',1.2,'FaceColor',[0 0 0])
	hold off
	f=get(gca,'Children');
	f=flipud(f(3:end));
	lgd=strsplit(num2str(LOOP.(LOOP.NamNoise),'%d '));
	% need to skip items in legend
	for i=1:length(lgd)
	    legendInfo{i} = ['SNR = ', lgd{i}, ' dB']; 
	end
	l=legend(f,legendInfo);	
	l.Location='northwest';
	set(gca,'XTickLabel',name);
	ylabel('Average localization error [mm]')
	% ax = findobj(gcf,'type','axes'); 
	print(gcf,[fullfile(pwd,filesep), strcat(myName,'_XRES_')],'-depsc')

	figure
	bar(RRES_sa',1.2);
	hold on
	bar(RRES_sstd',1.2,'FaceColor',[0 0 0])
	hold off
	f=get(gca,'Children');
	f=flipud(f(3:end));
	l=legend(f,legendInfo);
	l.Location='northwest';
	set(gca,'XTickLabel',name);
	ylabel('Average localization error [mm]')
	print(gcf,[fullfile(pwd,filesep), strcat(myName,'_RRES_')],'-depsc')

	clear gca;
	figure
	plot(XRANK,'-o','MarkerIndices',1:LOOP.totSimCount,'LineWidth',2);
	f=get(gca,'Children');
	f=flipud(f);
	l=legend(f,legendInfo);
	l.Location='northwest';
	axis([1 LOOP.totSimCount 1 sum(SETUP.SRCS(:,1))])
	xlabel('Simulation run')
	xticks([1:LOOP.totSimCount])
	ylabel('Rank selected')
	yticks([1:sum(SETUP.SRCS(:,1))])
	print(gcf,[fullfile(pwd,filesep), strcat(myName,'_XRANK_')],'-depsc')

      #+END_SRC
*** [p01] Main LOOP

    #+NAME: ./tg_s08_LOOP___01_Main_LOOP.m
    #+BEGIN_SRC matlab :session *MATLAB* :eval yes :tangle ./tg_s08_LOOP___01_Main_LOOP.m :results silent

      [LOOP.c01RR,LOOP.c02MM,LOOP.c03BB,LOOP.c04II] = ndgrid(LOOP.rngSimCount,LOOP.rngMesNoise,LOOP.rngBcgNoise,LOOP.rngIntNoise);
      LOOP.coords = [LOOP.c01RR(:),LOOP.c02MM(:),LOOP.c03BB(:),LOOP.c04II(:)];

      2 == length(size(squeeze(double.empty([[length(LOOP.rngMesNoise),length(LOOP.rngBcgNoise),length(LOOP.rngIntNoise)],0]))))
      if 2 == length(size(squeeze(double.empty([[length(LOOP.rngMesNoise),length(LOOP.rngBcgNoise),length(LOOP.rngIntNoise)],0]))))
	  if max([length(LOOP.rngMesNoise),length(LOOP.rngBcgNoise),length(LOOP.rngIntNoise)]) == min([length(LOOP.rngMesNoise),length(LOOP.rngBcgNoise),length(LOOP.rngIntNoise)]),
	      LOOP.IdxNoise = 0;
	  else,
	      [~,LOOP.IdxNoise] = max([length(LOOP.rngMesNoise),length(LOOP.rngBcgNoise),length(LOOP.rngIntNoise)]);
	  end;
	  if LOOP.IdxNoise == 1,
	      LOOP.NamNoise = 'rngMesNoise';
	  elseif LOOP.IdxNoise == 2,
	      LOOP.NamNoise = 'rngBcgNoise';
	  elseif LOOP.IdxNoise == 3,
	      LOOP.NamNoise = 'rngIntNoise';
	  else,
	      LOOP.NamNoise = 'AnyNoise';
	  end;
      else,
	  error('CYBERCRAFT:: ensure that only one noise type is vector and the remaining are scalars');
      end;
      LOOP.progress = []; LOOP.table_arrC = []; LOOP.table_varN = {}; LOOP.table_rowN = {};

      for LoopSimCount = 1:LOOP.totSimCount
	  LoopIterCount = [];      
	  run('./tg_s07_BATCH___02_Signals.m');
	  run('./tg_s07_BATCH___03_Leadfields.m');
	  for LoopIntNoise = 1:length(LOOP.rngIntNoise)
	      for LoopBcgNoise = 1:length(LOOP.rngBcgNoise)
		  for LoopMesNoise = 1:length(LOOP.rngMesNoise)

		      % used by localization
		      LoopIterCount = [LoopIterCount;LoopMesNoise,LoopBcgNoise,LoopIntNoise];
		      % used by reconstruction 
		      LOOP.progress   = [LOOP.progress;LoopSimCount,LOOP.rngMesNoise(LoopMesNoise),LOOP.rngBcgNoise(LoopBcgNoise),LOOP.rngIntNoise(LoopIntNoise)];
		      fprintf('\n');
		      disp('============================')
		      disp(['CurrIter: ',num2str(size(LOOP.progress,1)),' of ',num2str(length(LOOP.coords))])
		      disp('----------------------------')
		      disp(['SimCount: ',num2str(LoopSimCount),' of ',num2str(LOOP.totSimCount)])
		      disp(['MesNoise: ',num2str(LoopMesNoise),' of ',num2str(length(LOOP.rngMesNoise)),' (',num2str(LOOP.rngMesNoise(LoopMesNoise)),' dB in ' ,'[ ',num2str(LOOP.rngMesNoise),' ])'])
		      disp(['BcgNoise: ',num2str(LoopBcgNoise),' of ',num2str(length(LOOP.rngBcgNoise)),' (',num2str(LOOP.rngBcgNoise(LoopBcgNoise)),' dB in ' ,'[ ',num2str(LOOP.rngBcgNoise),' ])'])
		      disp(['IntNoise: ',num2str(LoopIntNoise),' of ',num2str(length(LOOP.rngIntNoise)),' (',num2str(LOOP.rngIntNoise(LoopIntNoise)),' dB in ' ,'[ ',num2str(LOOP.rngIntNoise),' ])'])
		      disp('----------------------------')
		      disp(['SETUP.n00: ',num2str(SETUP.n00)])
		      disp(['SETUP.K00: ',num2str(SETUP.K00)])
		      % chkSim___tg_s01_PRE___001(SETUP)
		      disp('============================')
		      SETUP.SINR    = LOOP.rngIntNoise(LoopIntNoise);
		      SETUP.SBNR    = LOOP.rngBcgNoise(LoopBcgNoise);
		      SETUP.SMNR    = LOOP.rngMesNoise(LoopMesNoise);

		      run('./tg_s07_BATCH___04_Preparations.m');

		      switch SETUP.supSwitch
			case {'rec'}

			  % Reconstruction
			  disp('Reconstruction')

			  run('./tg_s07_BATCH___05_Filters.m');

			  % removes Original and Dummy signals from Figures
			  if(SETUP.SHOWori==0)
			      rec_res.table_rowN = rec_res.table_rowN(3:end);
			      rec_res.table_arrC = rec_res.table_arrC(3:end,:);
			  end

			  LOOP.table_rowN{size(LOOP.progress,1)} = rec_res.table_rowN;
			  LOOP.table_varN{size(LOOP.progress,1)} = rec_res.table_varN;
			  LOOP.table_arrC(:,:,LoopSimCount,LoopMesNoise,LoopBcgNoise,LoopIntNoise) = rec_res.table_arrC;

			case {'loc'}

			  % Localization
			  disp('Localization')

			  run('tg_s07_BATCH___05_Localizers.m');

			  % use LOOP.RES.ITERATION{x}{y}(z) to access results struct for 
			  % z-th activity index (z=1,2,...,8 with z=1 for MAI(2011) etc., as defined in !run MAI and MPZ Locallizers) 
			  % y-th level of SNR of type considered (interference, biological, measurement)
			  % x-th iteration of the experiment
			  LOOP.RES.ITERATION{LoopSimCount}{size(LoopIterCount,1)} = RES;

			  % temporary helper code
			  % save(['RES__','ITERATION_',num2str(LoopSimCount),'__SINR_',num2str(SETUP.SINR),'__SBNR_',num2str(SETUP.SBNR),'__SMNR_',num2str(SETUP.SMNR)],'RES')

			otherwise
			  warning('Unexpected supSwitch value. No action taken, only data generated.')
		      end
		  end
	      end
	  end
      end

      switch SETUP.supSwitch
	case {'rec'}

	  if isequal(LOOP.table_rowN{:}), LOOP.table_rowN = LOOP.table_rowN{1}; else, error('CYBERCRAFT:: row names are inconsistent'); end;
	  if isequal(LOOP.table_varN{:}), LOOP.table_varN = LOOP.table_varN{1}; else, error('CYBERCRAFT:: var names are inconsistent'); end;

	case {'loc'}

	otherwise
	  warning('Unexpected supSwitch value. No action taken, only data generated.')
      end

    #+END_SRC

* DONE Supplementaries
** DOIN [z00] Auxiliary functions
*** =rawTimeSeries3d2d=

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawTimeSeries3d2d.m :results silent

function y = rawTimeSeries3d2d(x)

% y = rawTimeSeries3d2d(x)
% NB.: this function asumes that trials are along third dimention.

    y = reshape(permute(x,[1,3,2]),[],size(x,2),1);

end

    #+END_SRC

*** =rawRad2Deg=

    Convert angles from radians to degrees.

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawRad2Deg.m :results silent

function angDeg = rawRad2Deg(angRad)
% Convert angs from radians to degrees
angDeg = (180/pi) * angRad;

    #+END_SRC

*** =rawIsStableMVAR=

    Check whether the MVAR model is stable.

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawIsStableMVAR.m :results silent

function [H,varargout] = rawIsStableMVAR(A00,varargin)
    if isempty(varargin)
        STAB = 1;
        disp(['Checking if eigenvalues are inside unit circle.']);
    elseif length(varargin) == 1
        STAB = varargin{1};
        disp(['Checking if eigenvalues are inside circle with radius of ',num2str(STAB) ,'.']);
    else
        error('Too many input arguments');
    end
    S00          = size(A00,1);                  % dimension of state vectors
    P00          = size(A00,2)/S00;                % order of process
    tmp_lambda   = eig([A00; eye((P00-1)*S00) zeros((P00-1)*S00,S00)]);
    H            = ~any(abs(tmp_lambda)>STAB);
    varargout{1} = max(abs(tmp_lambda));
end

    #+END_SRC

*** =rawSize=

    Get size of an array (plain).

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawSize.m :results silent

function [s] = rawSize(x,varargin)
    if isempty(varargin)
        s = size(x);
    else
       s = [];
       for ii = 1:length(varargin{1})
           s = [s,size(x,ii)];
       end
    end
end

    #+END_SRC

*** =rawImgSC=

    Pretty table made from 2d array.

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawImgSC.m :results silent

function [h1,h2] = rawImgSC(mat,varargin)
    h1 = imagesc(mat);
    [x,y] = ndgrid(1:size(mat,2),1:size(mat,1));
    mat = mat';
    if ~isempty(varargin)
        fSize = varargin{1};
    else
       fSize = 18;
    end
    h2 = text(x(:),y(:),strtrim(cellstr(num2str(mat(:),'%.2f'))),'HorizontalAlignment','center','FontSize',fSize);
    if 0
        mat = [reshape(-12:12,5,5);reshape(-12:12,5,5);reshape(-12:12,5,5)]
        figure(17);clf;rawImgSC(mat);colorbar;
    end
end

    #+END_SRC

*** =rawHotColdColorMap=

    Pretty color-map for MATLAB figures.


    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawHotColdColorMap.m :results silent

function [y] = rawHotColdColorMap(x)
% [y] = rawHotColdColorMap(x)
% Suggested rationalization: base/2
    if 0
        x = 255
        y = rawHotColdColorMap(x)
    end
    y = [linspace(0,1,x)',linspace(0,1,x)',linspace(1,1,x)'];
    y = [y;fliplr(flipud(y))];
    % imagesc(permute(y,[1,3,2]))
end

    #+END_SRC

*** =rawPlotPDC=

    Plot PDC values.

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawPlotPDC.m :results silent

function rawPlotPDC(pdcData,accf)
    tmp_row  = size(pdcData,1);
    tmp_col  = size(pdcData,2);
    tmp_cnt  = 0;
    for i = 1:tmp_row
        for j = 1:tmp_col
            tmp_cnt = tmp_cnt+1;
            subplot(tmp_row,tmp_col,tmp_cnt)
            area(accf,squeeze(pdcData(i,j,:)))
            axis tight
            ylim([0,1])
        end
    end
end

    #+END_SRC

*** =rawAdjTotSNRdB=

    Adjust the SNR by changing the amplitude of the noise time-series.

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawAdjTotSNRdB.m :results silent

function [y] = rawAdjTotSNRdB(x01,x02,newSNR)
    y = ((x02 / rawNrm(x02)) * rawNrm(x01)) / (db2pow(0.5*newSNR));
end

    #+END_SRC

*** =rawNrm=

    Compute norm of the random vector.

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawNrm.m :results silent

function [y] = rawNrm(x)
 % Compute norm of the random vector
     y = sqrt(sum(rawMom(x,2),2));
 end

    #+END_SRC

*** =rawMom=

    Compute raw moment of the k-th order for the random vector.

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawMom.m :results silent

function [M] = rawMom(x01,k)
% Compute raw moment of the k-th order for the random vector
    M = mean(x01.^k, 1);
end

    #+END_SRC

*** =rawSNRdB=

    Compute raw moment of the k-th order for the random vector.

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawSNRdB.m :results silent

function [y] = rawSNRdB(x1,x2)
% Compute SNR in dB for two random vectors (signal and noise)
    y = 20.*log10(norm(x1,'fro')./norm(x2,'fro'));
end

    #+END_SRC

*** =ccrender=

    Prettify 3d MATLAB figures.

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/ccrender.m :results silent

function [] = ccrender(varargin)

%%
% CCRENDER --- figure rendering settings, usage:
%
% ccrender( axisLimits, finish, origin, originLimits )
%
%  - axisLimits [optional]:
%      limits display to specified range defined using either:
%        - vector of length 2 (all axes have same limits)
%          [ xyzLim_1, xyzLim_2 ],
%        - vector of length 6 (each axis has its own limits)
%          [ xLim_1, xLim_2, yLim_1, yLim_2, zLim_1, zLim_2 ];
%  - finish [optional]:
%      changes rendering scheme;
%  - origin [optional]:
%      true/false;
%  - originLim [optional]:
%      a number.
%
% Questions and comments: nikadon@gmail.com
%
% TODO:
% - add "try" envelopes for major blocks
%
%

%%

p = inputParser;

defAxisLimits = [];
%   chkAxisLimits = @(x) validateattributes(x,{'numeric'},{'vector','numel',2});
chkAxisLimits = @(x) isnumeric(x) && isvector(x) && ( length(x) == 2 || length(x) == 6 );
addOptional(p,'axisLimits',defAxisLimits,chkAxisLimits)

defOriginLimits = [];
chkOriginLimits = @(x) isnumeric(x) && isvector(x) && ( length(x) == 2 || length(x) == 6 );
addOptional(p,'originLimits',defOriginLimits,chkOriginLimits)

defEqualize = false;
chkEqualize = @(x) validateattributes(x,{'logical'},{});
addParameter(p,'equalize',defEqualize,chkEqualize)

defOrigin = true;
chkOrigin = @(x) validateattributes(x,{'logical'},{});
addParameter(p,'origin',defOrigin,chkOrigin)

defLabels = true;
chkLabels = @(x) validateattributes(x,{'logical'},{});
addParameter(p,'labels',defLabels,chkLabels)

defFinish = 'matte';
vldFinish = {'glossy','matte'};
chkFinish = @(x) any(validatestring(x,vldFinish));
addParameter(p,'finish',defFinish,chkFinish)




s = dbstack;


p.KeepUnmatched = true;


parse(p,varargin{:})

tmp_gcf = gcf

ccdisp('INIT: ccrender')
ccdisp(['Updating figure: ',num2str(tmp_gcf.Number),'...'])




if ~isempty(p.UsingDefaults)
    ccdisp([ s(1).name,': Using defaults:'])
    disp(p.UsingDefaults')
end
if ~isempty(p.Results)
    ccdisp([ s(1).name,': Results:'])
    disp(p.Results')
end
if ~isempty(fieldnames(p.Unmatched))
    ccdisp([ s(1).name,': Extra (unmatched) inputs:'])
    disp(p.Unmatched')
end





%%

if length(p.Results.axisLimits) == 2
    if p.Results.axisLimits(1) >= p.Results.axisLimits(2)
        help ccrender;
        error('CYBERCRAFT: axis limits [xyzLim_1,xyzLim_2] must be increasing.');
    else
        axis([ p.Results.axisLimits(1) p.Results.axisLimits(2) p.Results.axisLimits(1) p.Results.axisLimits(2) p.Results.axisLimits(1) p.Results.axisLimits(2) ]);
    end
elseif length(p.Results.axisLimits) == 6
    if p.Results.axisLimits(1) >= p.Results.axisLimits(2)  ||  p.Results.axisLimits(3) >= p.Results.axisLimits(4)  ||  p.Results.axisLimits(5) >= p.Results.axisLimits(6)
        help ccrender;
        error('CYBERCRAFT: axis limits [xLim_1,xLim_2,yLim_1,yLim_2,zLim_1,zLim_2] must be increasing for each axis.');
    else
        axis([ p.Results.axisLimits(1) p.Results.axisLimits(2) p.Results.axisLimits(3) p.Results.axisLimits(4) p.Results.axisLimits(5) p.Results.axisLimits(6) ]);
    end
end
%{
xlim([-15,15])
ylim([-15,15])
zlim([-15,15])
axis square
axis equal
axis auto
%}

if p.Results.equalize == true
    xMin = min(xlim);
    xMax = max(xlim);
    yMin = min(ylim);
    yMax = max(ylim);
    zMin = min(zlim);
    zMax = max(zlim);
    xyzMin = min( [ xMin yMin zMin ] );
    xyzMax = max( [ xMax yMax zMax ] );
    axis([ xyzMin xyzMax xyzMin xyzMax xyzMin xyzMax ]);
end


if p.Results.labels == true
    xlabel('X-axis')
    ylabel('Y-axis')
    zlabel('Z-axis')
end

set(gca,'DataAspectRatio',   [1 1 1])
set(gca,'PlotBoxAspectRatio',[1 1 1])

camproj('perspective'); % perspective | ortographic

%{
set(gca,'Projection','ortographic') % problem here?
                                    %}

set(gca,'CameraViewAngle',10)

az = 60;
el = 30;
view(az, el);

if strcmp(p.Results.finish,'glossy')
    lighting gouraud;
    material shiny;
    camlight;
    %{
    shading interp;
    light;
    lighting phong;
    %}
end

axis on;
box on;

%%
if p.Results.origin
    if length(p.Results.originLimits) == 2
        if p.Results.originLimits(1) > 0
            help ccrender;
            error('CYBERCRAFT: origin axis low limit must be negative number.')
        elseif p.Results.originLimits(2) < 0
            help ccrender;
            error('CYBERCRAFT: origin axis higher limit must be positive number.')
        else
            xMin = p.Results.originLimits(1);
            xMax = p.Results.originLimits(2);
            yMin = p.Results.originLimits(1);
            yMax = p.Results.originLimits(2);
            zMin = p.Results.originLimits(1);
            zMax = p.Results.originLimits(2);
        end
    elseif length(p.Results.originLimits) == 6
        if p.Results.originLimits(1) > 0  ||  p.Results.originLimits(3) > 0  || p.Results.originLimits(5) > 0
            help ccrender;
            error('CYBERCRAFT: origin axis low limit must be negative number.')
        elseif p.Results.originLimits(2) < 0  ||  p.Results.originLimits(4) < 0  ||  p.Results.originLimits(6) < 0
            help ccrender;
            error('CYBERCRAFT: origin axis higher limit must be positive number.')
        else
            xMin = p.Results.originLimits(1);
            xMax = p.Results.originLimits(2);
            yMin = p.Results.originLimits(3);
            yMax = p.Results.originLimits(4);
            zMin = p.Results.originLimits(5);
            zMax = p.Results.originLimits(6);
        end
    else
        xMin = min(xlim);
        xMax = max(xlim);
        yMin = min(ylim);
        yMax = max(ylim);
        zMin = min(zlim);
        zMax = max(zlim);
    end


    hold on




    try
        ccfgFigH = evalin( 'base', 'ccfgFigH' );
    catch
        ccfgFigH = [];
    end



    try
        tmp_fields = fieldnames(ccfgFigH(tmp_gcf.Number).origin);
        for ii = 1:length(tmp_fields)
            try
                delete(ccfgFigH(tmp_gcf.Number).origin.(tmp_fields{ii}));
            catch tmp_catched_2
                fprintf('\n\nPSEUDO-WARNING[2]: %s\n\n',tmp_catched_2.message);
            end
        end
    catch tmp_catched_1
        fprintf('\n\nPSEUDO-WARNING [1]: %s\n\n',tmp_catched_1.message);
    end


    ccfgFigH(tmp_gcf.Number).origin.xp = quiver3( 0, 0, 0, 1.1*xMax, 0,    0,    'color', [1 0 0], 'LineWidth', 1.0, 'LineStyle', '-' ); hold on;
    ccfgFigH(tmp_gcf.Number).origin.yp = quiver3( 0, 0, 0, 0,    1.1*yMax, 0,    'color', [0 1 0], 'LineWidth', 1.0, 'LineStyle', '-' ); hold on;
    ccfgFigH(tmp_gcf.Number).origin.zp = quiver3( 0, 0, 0, 0,    0,    1.1*zMax, 'color', [0 0 1], 'LineWidth', 1.0, 'LineStyle', '-' ); hold on;

    ccfgFigH(tmp_gcf.Number).origin.xn = quiver3( 0, 0, 0, xMin, 0,    0,    'color', [1 0 0], 'LineWidth', 1.0, 'LineStyle', '-.', 'ShowArrowHead', 'off' ); hold on;
    ccfgFigH(tmp_gcf.Number).origin.yn = quiver3( 0, 0, 0, 0,    yMin, 0,    'color', [0 1 0], 'LineWidth', 1.0, 'LineStyle', '-.', 'ShowArrowHead', 'off' ); hold on;
    ccfgFigH(tmp_gcf.Number).origin.zn = quiver3( 0, 0, 0, 0,    0,    zMin, 'color', [0 0 1], 'LineWidth', 1.0, 'LineStyle', '-.', 'ShowArrowHead', 'off' ); hold on;

    ccfgFigH(tmp_gcf.Number).origin.xt = text(   1.2*xMax, 0,        0,        'OX', 'color', [1 0 0] );
    ccfgFigH(tmp_gcf.Number).origin.yt = text(   0,        1.2*yMax, 0,        'OY', 'color', [0 1 0] );
    ccfgFigH(tmp_gcf.Number).origin.zt = text(   0,        0,        1.2*zMax, 'OZ', 'color', [0 0 1] );

    assignin( 'base', 'ccfgFigH', ccfgFigH );

end


%%


xLimits = xlim;
yLimits = ylim;
zLimits = zlim;




ccdisp( [ 'xLimits: [', num2str(xLimits(1)), ', ', num2str(xLimits(2)), ']' ] )
ccdisp( [ 'yLimits: [', num2str(yLimits(1)), ', ', num2str(yLimits(2)), ']' ] )
ccdisp( [ 'zLimits: [', num2str(zLimits(1)), ', ', num2str(zLimits(2)), ']' ] )

grid on;
rotate3d on;
hold on;



ccdisp(['Figure: ',num2str(tmp_gcf.Number),' has been updated!'])


ccdisp('DONE: ccrender')

return

    #+END_SRC

*** =ccdisp=

    Just throw some text to MATLAB console.

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/ccdisp.m :results silent

function [ ] = ccdisp( str )

% CCDISP display as CYBERCRAFT
%
% Use
%
%   ccdisp( str )
%

% add here "try + catch" OR "ifexists"

disp( ['CYBERCRAFT: ', str ] );

return

    #+END_SRC

*** =rawCartToSph=

    Convert cartesian to spherical coordinates (row-wise).

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawCartToSph.m :results silent

function [y] = rawCartToSph(x)
% Convert cartesian to spherical coordinates (row-wise)
    if ~isempty(x)
        for ii = 1:size(x,1)
            [y(ii,1),y(ii,2),y(ii,3)] = cart2sph(x(ii,1),x(ii,2),x(ii,3));
        end
    else
        y = x;
    end
end

    #+END_SRC

*** =rawSphToCart=

    Convert spherical to cartesian coordinates (row-wise).

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawSphToCart.m :results silent

function [y] = rawSphToCart(x)
% Convert spherical to cartesian coordinates (row-wise)
    if ~isempty(x)
        for ii = 1:size(x,1)
            [y(ii,1),y(ii,2),y(ii,3)] = sph2cart(x(ii,1),x(ii,2),x(ii,3));
        end
    else
        y=x;
    end
end

    #+END_SRC

*** =rawFixStrJoin=

    Some EEGLab plugins (e.g., cleanline) contain function =strjoin=
    which overwrites the original MATLAB function. The following
    function overshadows the any external =strjoin= functions with the
    original one.

    #+BEGIN_SRC matlab :session *MATLAB* :eval no :tangle ./fun/rawFixStrJoin.m :results silent

function [] = rawFixStrJoin()
    disp('=== Before ===')
    which('strjoin','-all')
    addpath([matlabroot,'/toolbox/matlab/strfun/'])
    disp('=== After ===')
    which('strjoin','-all')
end

    #+END_SRC

** DOIN [z01] Handy scripts

   #+NAME: ./tg_z01_hs___rem_all_tangled_files.sh
   #+BEGIN_SRC sh :shebang "#!/bin/bash" :eval yes :results silent :tangle ./tg_z01_hs___rem_all_tangled_files.sh

rm ./tg_*
rm ./fun/*

   #+END_SRC

** OPEN [z02] Remarks and observations
*** OPEN General

    - Principal sub-matrices in MVAR model are reconstructed quite
      robustly from IntNoise signal (that may have lowered dimensionality)
      when compared with SigIn
    - The above also aplies to PDC
** DONE [z91] GNU Emacs specific settings

   #+BEGIN_SRC emacs-lisp :eval yes :results output silent

(setq org-image-actual-width '(800))

(setq org-todo-keywords '(
                           (sequence "TODO(t)" "|" "DONE(d)" "SRTD" "NRML")
                           (sequence "NEW!" "FIX!" "WTF!" "DUE!" "CHCK" "REDO" "TOSS" "TEST" "OPEN" "DOIN"  "|" "CNCL" "NVRM")
                           (sequence "SOON" "LATE" "NEXT" "|" "OPTN" "WAIT" "HOLD")
                           (sequence "!run" "!std" "!opt" "!fun" "!fig" "|")
                           (sequence "CYCL" "|")
                           ))

   #+END_SRC
* TODO Things to be done shortly
** Major
*** DONE check if clear all is not disturbing batch-run /// loop-run !!!
*** DONE add evoked time-locked response

    As per:
    - cite:g2010---Hui-et-al---Identifying-true-cortical-interactions-in-MEG-using-the-nulling-beamformer

*** DONE SETTINGS parser
*** DONE Include  =SETUP.ERPs= in BATCH and LOOP
*** DONE Fix SETUP.ERPs generation (match array sizes)

** Minor
*** TODO ERPs figure
*** TODO Simulations work flow diagram

    #+NAME: ./tg_i01_Introduction___01_Simul_diagram.m
    #+BEGIN_SRC ditaa :file ./tg_i01_Gen_intro___01_Simul_diagram.png :cmdline -E :tangle ./tg_i01_Gen_intro___01_Simul_diagram.ditaa

/--------------------------\       /--------------------------\
| sel_atl.mat              |       | ROIs                     |
|--------------------------|       |--------------------------|
| sel_ele.mat              *       |                          *
|--------------------------|       |--------------------------|
| sel_ele.mat              |       |                          |
|--------------------------|       |--------------------------|
| sel_geo_deep_thalami.mat |       |                          |
|--------------------------|       |--------------------------|
| sel_mri00.mat            |       |                          |
|--------------------------|       |--------------------------|
| sel_msh.mat              |       |                          |
|--------------------------|       |--------------------------|
| sel_vol.mat              |       |                          |
\--------------------------/       \--------------------------/






/----------------------------\
|                            |
| Brain triangulation        |
|                            |
+----------------------------+                   /----------------------------\
|                            |                   |                            |
| Skull triangulation        *------------------>* Volume conduction model    |
|                            |                   |                            |
+----------------------------+                   +----------*-----------------+
|                            |                              |
| Scalp Triangulation        |                              |
|                            |                              |
\----------------------------/                              |
                                                            |
/----------------------------------------\                  |
|                                        |                  |
| Cortex parcellation                    |                  |
|                                        |                  |
\----------------*-----------------------/                  |
                 |                                          |
                 |                                          |
                 V                                          V
/----------------------------------------\       /----------*--------------------\
|                                        |       |                               |
| ROIs random sampling                   *------>*  Leadfield computation        |
|                                        |       |                               |
\-----------*----------------------------/       \-------------------------------/
            |                                               |
            |   /------------------------\                  |
            |   |                        |                  |
            +-->*   Volume cond. model   |                  |
                |                        |                  V
                +------------------------+       /----------*--------------------\
                |                        |       |                               |
                | Grid (ROI based)       *<------* Volume lookup (ROI)           |
                |                        |       |                               |
                +------------------------+       \-------------------------------/
                |                        |
                | Leadfields             |
                |                        |
                \------------------------/

    #+END_SRC
    #+CAPTION: Diagram describing outline of the simulations work flow
    #+LABEL: fig:org_fig
    #+ATTR_LATEX: :width 0.5\textwidth

    #+RESULTS:
    [[file:./tg_i01_Gen_intro___01_Simul_diagram.png]]
*** TODO Neuro-scientific background introduction