clear

ft_defaults
% this returns the location where FieldTrip is installed
[ftver, ftdir] = ft_version;


%% Loading some preprocessed data

isubj = 51
%% Select your parameters to save the data in one file

format.baseline = [-0.2 0];
format.lpfreq = 30;
%format.hpfreq = 2;
format.resamplefs = 200;
format.detrend = 'no';
format.trigger_ref = 'target'; % what is the trigger used as a reference


data_file_name = ['dformat_s',num2str(isubj)];
% Define globals
global G;
G = [];

%% General paths
exp_main_path = '/project/3018029.06';

SUBJS = {['0',num2str(isubj)]}
curSubj = SUBJS{1}

addpath(genpath( [exp_main_path filesep 'Matlab' filesep 'MEG3_scripts'])); % general processing scripts in this exp.

G.subjPath = [exp_main_path filesep 'MEG' filesep 'AVE_MEG3' filesep curSubj filesep]; % main subject folder
G.outputPath = [G.subjPath ,'func', filesep]; % subject saving artifact and preprocessing matrixes
G.cfgbehavPath = [exp_main_path filesep 'Behavioral' filesep 'AVE_MEG3' filesep curSubj filesep] ; % behavioral data of participant

subjscriptsPath = [G.subjPath ,'analysis', filesep]; %  participant specific scripts


%% loading preproc data
load([G.outputPath 'data_preproc_clean.mat'])




%% Reading behavioral data to append in trialinfo
dirs = dir([G.cfgbehavPath, '*_CM_exp.mat']);
if (length(dirs) ~= 1)
    warning(['Not exactly one MEG-file found for Subj', num2str(curSubj), '!']);
end

load([G.cfgbehavPath, dirs.name])

%% The 
part below is just to recover some time information (not used by the
% moment in the source reconstruction protocol)
% concatenating blocks
behav_matrix = [];
for i = 1 : size(Experimentalsession,2)
    behav_matrix = [behav_matrix ;   Experimentalsession(i).results.responses_vector];
end


time_data = [];

for iblock =   1 : size(Experimentalsession,2)
       % select data from iblock in preprocessed data
       trials_block = find(data_preproc_clean.trialinfo(:,2) == iblock);
        % select trial indexes
       trials_sel = data_preproc_clean.trialinfo(trials_block,3);
       % select data from iblock sin total behav data
       trials_block = find(cell2mat(behav_matrix(:,1)) == iblock);
       total_trials = cell2mat(behav_matrix(trials_block,2));
       % find those trials that were selected during preproc
       filtered_trials = ismember(total_trials,trials_sel);
       d = cell2mat(behav_matrix(filtered_trials,3:7));
       time_data = [time_data; d];
       size(time_data)    
end
%%%%%%%%%%%%%%%


%% Process dataclean and downdsampled here 



% Downsampling the data
cfg                 = [];
cfg.resamplefs = 300; % 200 hz is 5 ms time-res, 20Hz is 50 ms
cfg.detrend    = 'no';
data = ft_resampledata(cfg, data_preproc_clean);


% Baseline correction 
% cfg                 = [];
% cfg.baseline = [-0.2, 0]; % 200 hz is 5 ms time-res, 20Hz is 50 ms
% data = ft_timelockbaseline(cfg, data);

cfg                 = [];  % to keep the data in raw format, use ft_preproc instead of baseline correction
cfg.demean = 'yes';
cfg.baselinewindow = [-0.2, 0]; % 200 hz is 5 ms time-res, 20Hz is 50 ms
data = ft_preprocessing(cfg, data);


cfg                 = [];
cfg.channel       = {'MEG'};
data = ft_selectdata(cfg, data)

% dataframe 4 column % 0 is noise 1 is signal 2  is signal2x
% dataframe 5 column %  % 0 is V, 1 is AV
% dataframe 7 column % % contrast
% dataframe 8 column %  % 1 is NO, 2 is YEs and, 0 is NoResponse
% dataframe 9 column % % incorrect correct
% dataframe 10 column %  % RT
% dataframe 11 column %  % time template
% dataframe 12 column %  % time 1stcue
% dataframe 13 column %  % time stim  (you can calculate the time of the cues relative to the stim time)
% dataframe 14 column %  % time 2ndcue
% Validate trialinfo


%% TL data selection to verification

cfg = [];
%cfg.keepindividual = 'yes';
cfg.covariance       = 'yes';
%cfg.keeptrials = 'yes';
cfg.latency = [-0.2, 1.25]
cfg.channel = {'MEG'}
tl_rf_data = ft_timelockanalysis(cfg, data);


figure;
cfg = [];
cfg.channel =  {'MLO*','MZO*','MRO*'} %{'MEG'}; %'MEG'; {'MLO*','MZO*','MRO*'}
cfg.baseline = [-0.2, 0];
cfg.linewidth  = 1;
cfg.parameter = 'avg';
cfg.linewidth = [2.0];
%  cfg.ylim = [-6e-14,14e-14]
cfg.xlim = [-0.2,1.5];
cfg.layout = 'CTF275_helmet.lay';
ft_singleplotER(cfg, tl_rf_data);



%% Contrasting high contrast vs no visual stim condition

time_win = [-0.2, 1]
cfg                = [];
cfg.channel = {'MEG'}
cfg.latency = time_win; 
cfg.trials  = data.trialinfo(:, 4) == 0; %noise
v_data      = ft_selectdata(cfg, data);

cfg.trials  = data.trialinfo(:, 4) == 2; % high_contrast
cfg.latency = time_win;
av_data     = ft_selectdata(cfg, data);

%v_data.time = av_data.time; % alling time dimension

cfg      = [];
data_cmb = ft_appenddata(cfg, av_data, v_data);



%% SOund presentation
% 
 cfg                = [];
 %cfg.latency = [-0.2, 0.0]
 cfg.trials         = data.trialinfo(:, 5) == 0; %data.trialinfo == 0; trials = find(data.trialinfo(:, 4) == conds(icond)
 v_data           = ft_selectdata(cfg, data);
% 
% 
 cfg.trials         = data.trialinfo(:, 5) == 1;
% cfg.latency = [1.3, 1.5]
 av_data            = ft_selectdata(cfg, data);
 
% v_data.time = av_data.time % alling time dimension

cfg      = [];
data_cmb = ft_appenddata(cfg, av_data, v_data);

% give a number to each trial: 0 = baseline, 1 = experimental condition
%data_cmb.trialinfo = [zeros(size(av_data.trial,1), 1); ones(size(v_data.trial,1), 1)];


cfg = [];
%cfg.keepindividual = 'yes';
%cfg.covariance       = 'yes';
%cfg.keeptrials = 'yes';
%cfg.latency = [1.25, 0.25]
cfg.channel = {'MEG'}
tl_cmb_data = ft_timelockanalysis(cfg, data_cmb);



cfg = [];
%cfg.latency = [0.05,0.25]
cfg.covariance       = 'yes';
cfg.channel = {'MEG'}
tl_av_data = ft_timelockanalysis(cfg, av_data);
tl_v_data = ft_timelockanalysis(cfg, v_data);


tl_av_cmb_data = AVE_combine_planar(tl_av_data );

tl_v_cmb_data = AVE_combine_planar(tl_v_data );


cfg = [];
cfg.operation = 'subtract';
cfg.parameter = 'avg';
% cfg.parameter = 'individual';
diffav_tl_cm_data = ft_math(cfg, tl_av_cmb_data, tl_v_cmb_data);

%cfg = []
%cfg.avgovertime = 'yes'
%diffav_tl_data = ft_selectdata(cfg, diffav_tl_data)

figure;
cfg = [];
cfg.layout = 'CTF275.lay';
cfg.baseline = [-0.2, 0];
cfg.xlim = [-0.1 : 0.1 : 0.6];  
%cfg.xlim = [-0. : 0.1 : 0.9];  
cfg.interactive   = 'yes';
cfg.zlim = [-0.5e-13 0.5e-13];      % Set the 'color' limits.
ft_topoplotER(cfg, diffav_tl_cm_data);





% %% Prepare_data_in_different_formats for analyses in each subject
% cfg = [];
% cfg.filename = [analysis_Path,'mri_ctf']
% cfg.filetype = 'nifti';
% cfg.parameter = 'anatomy';
% ft_volumewrite(cfg, mri_ctf);
% %segmentedmri = ft_read_mri([outpath, filesep,'segmented.nii']);
% 
% %% saving to source folder
% cfg = [];
% cfg.filename = [analysis_Path,'segmented']
% cfg.filetype = 'nifti';
% cfg.parameter = 'anatomy';
% ft_volumewrite(cfg, segmentedmri);


%% Starting source reconstruction

exp_main_path = '/project/3018029.06';
analysis_Path = [exp_main_path filesep 'Matlab' filesep  'MEG3_scripts' filesep 'source_reconstruction' filesep]; % subjecf MEG data
path = '/Users/alex/Dropbox/TOOLBOXES/fieldtrip/template/anatomy/'

mri = ft_read_mri([ftdir, filesep, 'template', filesep, 'anatomy',filesep 'single_subj_T1_1mm.nii']);

%% allingning to ctf coordinates
cfg = [];
cfg.method = 'interactive';
cfg.coordsys = 'ctf';
mri_ctf = ft_volumerealign(cfg, mri);

save([analysis_Path 'mri_ctf.mat'], 'mri_ctf');

%% or load instead

load([analysis_Path 'mri_ctf.mat']);


%% segmenting t1
cfg          = [];
cfg.coordsys = 'ctf';
segmentedmri = ft_volumesegment(cfg,mri_ctf);


%segmentedmri = ft_read_mri([analysis_Path 'segmented.nii']);

% add anatomical information to the segmentation
segmentedmri.transform = mri_ctf.transform;
segmentedmri.anatomy   = mri_ctf.anatomy;

% Visualizing correct segmentation
cfg              = [];
cfg.funparameter = 'gray';
ft_sourceplot(cfg, segmentedmri);

save([analysis_Path 'segmentedmri_ctf.mat'], 'segmentedmri');

%% or load instead

load([analysis_Path 'segmentedmri_ctf.mat']);

%% preparing headmodel
cfg        = [];
%cfg.spmversion = 'spm12',
cfg.method = 'singleshell';
hdm        = ft_prepare_headmodel(cfg, segmentedmri);

save([analysis_Path  'hdm'], 'hdm')

%% or load instead

load([analysis_Path 'hdm.mat']);


%% and this is where the template source models are
templatedir = fullfile(ftdir, 'template', 'sourcemodel');

template = load(fullfile(templatedir, 'standard_sourcemodel3d10mm')); % 10mm spacing grid

%% inverse-warp the template grid to subject specific coordinates
cfg                = [];
%cfg.spmversion = 'spm12',
cfg.warpmni   = 'yes';
cfg.template  = template.sourcemodel;
cfg.nonlinear = 'yes'; % use non-linear normalization
cfg.mri            = mri_ctf;
%cfg.symmetry      = 'x'
%cfg.unit = 'mm'
sourcemodel        = ft_prepare_sourcemodel(cfg);


%% Lets visualize if the corregistration went well.

hdm_cm = ft_convert_units(hdm, 'cm');
figure; hold on     % plot all objects in one figure
%ft_plot_headshape(shape);

ft_plot_headmodel(hdm_cm, 'edgecolor', 'none')
alpha 0.4           % make the surface transparent
ft_plot_mesh(sourcemodel.pos(sourcemodel.inside,:));
ft_plot_sens(data.grad, 'style', 'ob');

% All looks good so far! (ES)

%% Preparing source model (leadfield)


cfg             = [];
cfg.grid        = sourcemodel;
cfg.headmodel   = hdm;
cfg.channel     = {'MEG'};
cfg.grad        = data_cmb.grad;
%cfg.normalize   = 'yes' % if not comparing two conditions then yes

sourcemodel_lf  = ft_prepare_leadfield(cfg,  data_cmb );


%% What I want to beam in time

% note: overwriting previously computed TL structures

cfg = [];
cfg.covariance       = 'yes';
cfg.channel = {'MEG'};
cfg.covariancewindow = [0.1 0.25];
tl_av_data = ft_timelockanalysis(cfg, av_data);
tl_v_data = ft_timelockanalysis(cfg, v_data);
tl_cmb_data = ft_timelockanalysis(cfg, data_cmb);


%% create spatial filter using the lcmv beamformer
cfg                  = [];
cfg.method           = 'lcmv';
cfg.sourcemodel      = sourcemodel_lf; % leadfield
cfg.headmodel        = hdm; % volume conduction model (headmodel)
cfg.lcmv.keepfilter  = 'yes';
%cfg.lcmv.fixedori    = 'yes'; % project on axis of most variance using SVD O
%cfg.lcmv.projectnoise  = 'yes'; %only when you dont waht to contrast two
%conditions
%cfg.lcmv.weightnorm    = 'nai'; %only when you dont waht to contrast two
%conditions
cfg.lcmv.lambda        = '15%'; % smoothness (high because not real structural data)
source               = ft_sourceanalysis(cfg,  tl_cmb_data);

%% beam pre- and poststim by using the common filter
cfg.sourcemodel.filter  = source.avg.filter;
source_av       = ft_sourceanalysis(cfg, tl_av_data);
source_v       = ft_sourceanalysis(cfg, tl_v_data);
%cfg.sourcemodel.filter  = source.avg.filter;
%source_tl_data       = ft_sourceanalysis(cfg,  tl_rf_data);

%%

cfg = [];
cfg.operation = 'log10(x1/x2)'; %'(x1-x2)/(x1+x2)' %'log10(x1/x2)';
cfg.parameter = 'avg.pow';
% cfg.parameter = 'individual';
diff_av_source = ft_math(cfg, source_av, source_v);

%diff_av_source.pos = template.sourcemodel.pos;
%diff_av_source.dim = template.sourcemodel.dim;


%%
cfg              = [];
cfg.parameter    = 'pow';
cfg.interpmethod = 'nearest';
source_int  = ft_sourceinterpolate(cfg,  diff_av_source  , mri_ctf);

%%   
cfg               = [];
cfg.method        = 'ortho';
cfg.funparameter  = 'pow';
cfg.location = 'min';%diff_av_source.pos(maxval_erf,:);
%cfg.latency  = 'prestim'
cfg.maskparameter = cfg.funparameter;
% cfg.funcolorlim   = 'auto';
% cfg.funcolorlim   = [0 2.3063e-07];
%cfg.opacitylim    = [0.03 0.06];
%cfg.opacitylim    = [1 3];
cfg.interactive   = 'yes';
% cfg.opacitymap    = 'rampup';
%cfg.avgovertime = 'yes'
ft_sourceplot(cfg,diff_av_source);

%% selecting grid points based on max/min values in power contrast

[maxval, maxvalindx] = max(diff_av_source.pow )

%[nullval, nullvalindx] = find(diff_av_source.pow > -0.0001 & diff_av_source.pow < 0.0001)

[minval, minvalindx] = min(diff_av_source.pow )

%% selecting grid points based on max/min values in absolute ERF contrast

% idea: compute absolute ERF contrast, average value in time window of
% interest, project this average difference value to the whole source
% space, select grid points there

% note that projecting to source space is a linear operation, so projecting
% a difference is the same as projecting both separately and then
% subtracting

% ERF difference at channel level
erf_diff = tl_av_data.avg - tl_v_data.avg;

% average in time window of interest
timinds = find(tl_av_data.time >= 0.25 & tl_av_data.time <= 0.5);
erf_diff = mean(erf_diff(:,timinds), 2);
erf_v = mean(tl_v_data.avg(:,timinds), 2);
erf_av = mean(tl_av_data.avg(:,timinds), 2);

% loop over grid points
% note we do this in a loop to avoid having to deal with source.inside indexing etc.
ngridpoint = numel(source.avg.filter);
sourcediffs = zeros(ngridpoint,1);
for k = 1:ngridpoint
  if isempty(source.avg.filter{k})
    % outside the brain
    sourcediffs(k) = nan;
  else
    thisfilt = source.avg.filter{k};
    % norm of source-projected difference score
    sourcediffs(k) = sqrt(sum((thisfilt * erf_diff).^2));
    
    % normalize by total power in two conditions
    %sourcediffs(k) = sourcediffs(k) ./ sqrt((sum((thisfilt * erf_v).^2) + sum((thisfilt * erf_av).^2))./2);
  end
end

[~,maxval_erf] = max(sourcediffs);
[~,minval_erf] = min(sourcediffs);



%% ERFs in a bunch of grid points of interest

% use indices of max and min val in power contrast, but also 3 randomly
% chosen brain grid points
randpoints = randsample(find(diff_av_source.inside), 2);
frontalnull = 4132;
allpoints = [maxvalindx; minvalindx; maxval_erf; minval_erf; frontalnull;  randpoints];
npoint = numel(allpoints);
nchan = numel(source.avg.label);

%% Recomputing source model taking into account the whole trial length
% note: overwriting previously computed TL structures

cfg = [];
cfg.covariance       = 'yes';
cfg.channel = {'MEG'};
%cfg.covariancewindow = [0 0.3];
tl_av_data = ft_timelockanalysis(cfg, av_data);
tl_v_data = ft_timelockanalysis(cfg, v_data);
tl_cmb_data = ft_timelockanalysis(cfg, data_cmb);


%% create spatial filter using the lcmv beamformer
cfg                  = [];
cfg.method           = 'lcmv';
cfg.sourcemodel      = sourcemodel_lf; % leadfield
cfg.headmodel        = hdm; % volume conduction model (headmodel)
cfg.lcmv.keepfilter  = 'yes';
%cfg.lcmv.fixedori    = 'yes'; % project on axis of most variance using SVD O
%cfg.lcmv.projectnoise  = 'yes'; %only when you dont waht to contrast two
%conditions
%cfg.lcmv.weightnorm    = 'nai'; %only when you dont waht to contrast two
%conditions
cfg.lcmv.lambda        = '15%'; % smoothness (high because not real structural data)
source               = ft_sourceanalysis(cfg,  tl_cmb_data);


% find the orientation of maximum power for all grid points of interest
% note that we multiply the max-power SVD with the spatial filter
% immediately here, so no need to project to max direction later on
single_ori_filters = zeros(npoint,nchan);
alldat = cat(2, data.trial{:});
for k = 1:npoint
  thisfilt = source.avg.filter{allpoints(k)};
  alldat_filt = thisfilt * alldat;
  [u,s,v] = svd(alldat_filt, 'econ');
  single_ori_filters(k,:) = u(:,1)' * thisfilt;
end
clear alldat; % save a lot of memory


% create timelock structure
tl_poi_v = [];
tl_poi_v.label = {'maxval pow', 'minval pow', 'maxval erf', 'minval erf', 'frontalnull', 'rand1', 'rand2'};
tl_poi_v.time = tl_v_data.time;
tl_poi_v.dimord = 'chan_time';

% project event-related data through filter
tl_poi_v.avg = single_ori_filters * tl_v_data.avg;

tl_poi_av = tl_poi_v;
tl_poi_av.avg = single_ori_filters * tl_av_data.avg;

%% plot


figure();
allax = [];
for k = 1:npoint
  allax(end+1) = subplot(3,3,k);
  plot(tl_poi_v.time, tl_poi_v.avg(k,:), tl_poi_av.time, tl_poi_av.avg(k,:));
  legend('absent', 'high contrast');
  title(tl_poi_v.label{k});
end

linkaxes(allax);


















%% Contrasting high contrast vs no visual stim condition

time_win_all = [-0.2, 0.75]
cfg         = [];
cfg.channel = {'MEG'}
cfg.latency = time_win_all; 
cfg.trials  = data.trialinfo(:, 4) == 0; %noise
x_data      = ft_selectdata(cfg, data);

cfg.trials  = data.trialinfo(:, 4) == 2; % high_contrast
cfg.latency = time_win_all;
xv_data     = ft_selectdata(cfg, data);

x_data.time = xv_data.time % alling time dimension

cfg      = [];
datax_cmb = ft_appenddata(cfg, xv_data, x_data);

cfg = [];
%cfg.keepindividual = 'yes';
cfg.covariance       = 'yes';
%cfg.keeptrials = 'yes';
%cfg.latency = [1.25, 0.25]
cfg.channel = {'MEG'}
tl_cmb_datax = ft_timelockanalysis(cfg, datax_cmb);


cfg = [];
%cfg.latency = [0.05,0.25]
cfg.covariance       = 'yes';
cfg.channel = {'MEG'}
tl_xv_data = ft_timelockanalysis(cfg, xv_data);
tl_x_data = ft_timelockanalysis(cfg, x_data);


% Now lets use the virtual channel to select activity from the whole data
% set. We have to generate a source model to do that
cfg              = [];
cfg.method       = 'lcmv';
cfg.headmodel    = hdm;
cfg.sourcemodel.pos     = sourcemodel.pos(minvalindx, :); % sourcemodel.pos([nullval(1)],:) %sourcemodel.pos([maxvalindx], :); %[7.1 -0.6 10.9]; %sourcemodel.pos([maxvalindx], :);
%cfg.sourcemodel.inside  = true(1,1); %true(2,1); % dependeing on the number of sources to localice
cfg.unit    = sourcemodel.unit;
cfg.lcmv.keepfilter = 'yes';
cfg.lcmv.lambda        = '15%'; % smoothness (high because not real structural data)
source_idx       = ft_sourceanalysis(cfg, tl_cmb_datax );


beamformer_vchan = source_idx.avg.filter{1};
%beamformer_gam_pow = source_idx.avg.filter{2};

cfg      = [];
data_cell = ft_appenddata(cfg, datax_cmb);



chansel = ft_channelselection('MEG', data_cell.label); % find MEG sensor names
chansel = match_str(data_cell.label, chansel);         % find MEG sensor indices
source_chan = [];
source_chan.label = {'av_diff__x', 'av_diff__y', 'av_diff__z'};
source_chan.time = data_cell.time;



for i=1:length(data_cell.trial)
   source_chan.trial{i} =  beamformer_vchan * data_cell.trial{i}(chansel,:);
end

%for i=1:length(data_cell.trial)
%   source_chan.trial{i} =  beamformer_vchan * squeeze(data_cell.trial{i});
%end

%cfg = [];
%cfg.viewmode = 'vertical';  % you can also specify 'butterfly'
%ft_databrowser(cfg, source_chan);

% Decomposinh the projections
timeseries = cat(2, source_chan.trial{:});
[u, s, v] = svd(timeseries, 'econ');

% Recombining the filtered data in one virtual channel

% this is equal to the first column of matrix V, apart from the scaling with s(1,1)
timeseriesmaxproj = u(:,1)' * timeseries;

vchanneldata = [];
vchanneldata.label = {'cortex'};
vchanneldata.time = data_cell.time;
for i=1:length(data_cell.trial)
  vchanneldata.trial{i} = u(:,1)' * beamformer_vchan * data_cell.trial{i}(chansel,:);
end



%cfg = [];
%cfg.viewmode = 'vertical';  % you can also specify 'butterfly'
%ft_databrowser(cfg, vchanneldata);


%% Lets plot the results

cfg         = [];
%cfg.channel = {'MEG'}
cfg.trials  = datax_cmb.trialinfo(:, 4) == 0; %noise
x_vchannel      = ft_timelockanalysis(cfg, vchanneldata)
cfg.trials  = datax_cmb.trialinfo(:, 4) == 2; % high_contrast
xv_vchannel     =  ft_timelockanalysis(cfg, vchanneldata);


figure;
cfg = [];
cfg.graphcolor    =  [0 0 0;0.0 0.6 0.9];
cfg.parameter = 'avg';
%cfg.baseline = [-0.2, 0];
%cfg.xlim = [-0. : 0.1 : 0.9]; 
cfg.ylim = [-0.03, 0.04]; 
%cfg.zlim = [-0.5e-13 0.5e-13];      % Set the 'color' limits.
ft_singleplotER(cfg, xv_vchannel, x_vchannel );



%av_vchannel.label = {'vchan'}


% give a number to each trial: 0 = baseline, 1 = experimental condition
% data_cmb.trialinfo = [zeros(size(av_data.trial,1), 1); ones(size(v_data.trial,1), 1)];


cfg = [];
cfg.operation = 'subtract';
cfg.parameter = 'avg';
% cfg.parameter = 'individual';
diffav_tl_data = ft_math(cfg, tl_av_data, tl_v_data);

cfg = []
cfg.avgovertime = 'yes'
diffav_tl_data = ft_selectdata(cfg, diffav_tl_data)

tl_rf_cmb_data = AVE_combine_planar(diffav_tl_data);




%%% TRYING TO REPEAT THE RESULTS USING THE SAME AMMOUNT OF TIME USED TO
%%% GENERATE THE SOURCE MODEL

[minval, minvalindx] = min(diff_av_source.pow )
[maxval, maxvalindx] = max(diff_av_source.pow )


%% Contrasting high contrast vs no visual stim condition



% Now lets use the virtual channel to select activity from the whole data
% set. We have to generate a source model to do that
cfg              = [];
cfg.method       = 'lcmv';
cfg.headmodel    = hdm;
cfg.sourcemodel.pos     = sourcemodel.pos(maxvalindx, :); % sourcemodel.pos([nullval(1)],:) %sourcemodel.pos([maxvalindx], :); %[7.1 -0.6 10.9]; %sourcemodel.pos([maxvalindx], :);
%cfg.sourcemodel.inside  = true(1,1); %true(2,1); % dependeing on the number of sources to localice
cfg.unit    = sourcemodel.unit;
cfg.lcmv.keepfilter = 'yes';
cfg.lcmv.lambda        = '15%'; % smoothness (high because not real structural data)
source_idx       = ft_sourceanalysis(cfg, tl_cmb_data);


beamformer_vchan = source_idx.avg.filter{1};
%beamformer_gam_pow = source_idx.avg.filter{2};

cfg      = [];
data_cell = ft_appenddata(cfg, data_cmb);



chansel = ft_channelselection('MEG', data_cell.label); % find MEG sensor names
chansel = match_str(data_cell.label, chansel);         % find MEG sensor indices
source_chan = [];
source_chan.label = {'av_diff__x', 'av_diff__y', 'av_diff__z'};
source_chan.time = data_cell.time;



for i=1:length(data_cell.trial)
   source_chan.trial{i} =  beamformer_vchan * data_cell.trial{i}(chansel,:);
end


%cfg = [];
%cfg.viewmode = 'vertical';  % you can also specify 'butterfly'
%ft_databrowser(cfg, source_chan);

% Decomposinh the projections
timeseries = cat(2, source_chan.trial{:});
[u, s, v] = svd(timeseries, 'econ');

% Recombining the filtered data in one virtual channel

% this is equal to the first column of matrix V, apart from the scaling with s(1,1)
timeseriesmaxproj = u(:,1)' * timeseries;

vchanneldata = [];
vchanneldata.label = {'cortex'};
vchanneldata.time = data_cell.time;
for i=1:length(data_cell.trial)
  vchanneldata.trial{i} = u(:,1)' * beamformer_vchan * data_cell.trial{i}(chansel,:);
end



%cfg = [];
%cfg.viewmode = 'vertical';  % you can also specify 'butterfly'
%ft_databrowser(cfg, vchanneldata);


%% Lets plot the results

cfg         = [];
%cfg.channel = {'MEG'}
cfg.trials  = data_cmb.trialinfo(:, 4) == 0; %noise
x_vchannel      = ft_timelockanalysis(cfg, vchanneldata)
cfg.trials  = data_cmb.trialinfo(:, 4) == 2; % high_contrast
xv_vchannel     =  ft_timelockanalysis(cfg, vchanneldata);


figure;
cfg = [];
cfg.graphcolor    =  [0 0 0;0.0 0.6 0.9];
cfg.parameter = 'avg';
%cfg.baseline = [-0.2, 0];
%cfg.xlim = [-0. : 0.1 : 0.9]; 
cfg.ylim = [-0.03, 0.04]; 
%cfg.zlim = [-0.5e-13 0.5e-13];      % Set the 'color' limits.
ft_singleplotER(cfg, xv_vchannel, x_vchannel );






    