function Source_dec2vf_subj_prep(cfg_fun)
% this script us used to decode the stimulus at different ROIs. It will
% generate to output files (one for diagonal ms at 100 hz and 1 for TGM at
% 30 hz)

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

nperms = 15;

%% Starting source reconstruction
exp_main_path = cfg_fun.exp_main_path;
analysis_Path = [exp_main_path filesep 'Matlab' filesep  'MEG_scripts' filesep 'source_reconstruction' filesep]; % subjecf MEG dat

% loading rois
%load([analysis_Path  'roi_gridpoints.mat'])
%% Loading subject data

%subject = 49

SUBJS = {['0',num2str(cfg_fun.subject)]}
curSubj = SUBJS{1}

G = [];
G.subjPath = [exp_main_path filesep 'MEG'  filesep curSubj filesep]; % main subject folder
G.inputPath = [G.subjPath ,'func', filesep, 'source_based']; % subject saving artifact and preprocessing matrixes
G.outputPath = [G.subjPath ,'func', filesep, 'source_based']; % subject saving artifact and preprocessing matrixes

% loading data project along the virtual channels for each ROI
load([G.inputPath , filesep, 'vchan_ROIS_100.mat'])

% loading roi and power difference bewtween sp and sa at different temporal
% rois.

%load([G.inputPath, filesep, 'source_tempERFwin_ROIS.mat'])
load([G.inputPath , filesep, 'vchan_template_fun_ROIS.mat'])

% looping over ROIs and performing the selected analyses and saving the
% data

nrois = length(ROI_hf)


%% not really using this for decoding
%roi_sizes = [1, 10, 25, 50, 100]; % repweat the analyses for different roi sizes

conds = unique(ROI_hf(1).vchan_data.trialinfo(:,4)); % noise visual
modals = unique(ROI_hf(1).vchan_data.trialinfo(:,5));
vfields = unique(ROI_hf(1).vchan_data.trialinfo(:,6)); % center periphery


%% LETs use a decoder in order to get the data


for i_roi = 2 % only visual ROI
    for i_vf = 1 : 2 % only visual ROI
        roi_name = ROI_hf(i_roi).name;
        % what combinations I want to decode
        % combination_decoded = [0, 1] %; 0, 2; 1, 2]; % no signal vs low contrast, no signal vs high contrast, low contrast vs high contras
        combination_names = {'N vs S',} %'Ns vs HC','LC vs HC'};
        
        
        cfg = [];
        xdata = ROI_hf(i_roi).vchan_data;
        trialinfo = xdata.trialinfo;
        
        % Here I apply the baseline correction
        % cfg                 = [];  % to keep the data in raw format, use ft_preproc instead of baseline correction
        % cfg.demean = 'yes';
        % xdata = ft_preprocessing(cfg, xdata);
        
        %% Timelockbaseline correction to remove slow drift modulations but without affecting the neural activity right befgore the response
        cfg = [];
        cfg.baseline = [-2.0, 0.0];
        dat = ft_timelockbaseline(cfg, xdata)
        
        % Reformatting data to timelock structure
        cfg                = [];
        cfg.trials =  trialinfo(:, 6) == vfields(i_vf);
        cfg.keeptrials         = 'yes'
        alldat(i_roi).vfield(i_vf).vchan_decode_data      = ft_timelockanalysis(cfg, dat );
        
        trialinfo = alldat(i_roi).vfield(i_vf).vchan_decode_data.trialinfo;
        
        %% Zscoredata
        for imod = 1 : 2;
            idxs =  trialinfo(:,5) == modals(imod);
            for it = 1 : length(alldat(i_roi).vfield(i_vf).vchan_decode_data.time)
                %alldat(i_roi).vchan_decode_data
                Zdata = alldat(i_roi).vfield(i_vf).vchan_decode_data.trial(idxs, :, it);
                Zdata =  (Zdata - mean(Zdata(:))) /   std(Zdata(:));
                
                alldat(i_roi).vfield(i_vf).vchan_decode_data.trial(idxs, :, it) = Zdata
            end
        end
        
        %      AV = alldat(i_roi).vchan_decode_data.trial( trialinfo(:,5) == 1, :, :);
        %      V = alldat(i_roi).vchan_decode_data.trial( trialinfo(:,5) == 0, :, :);
        %
        %     figure()
        %      hist(mean(alldat(i_roi).vchan_decode_data.trial( trialinfo(:,5) == 1, :, 132),2));
        %      figure()
        %      hist(mean(alldat(i_roi).vchan_decode_data.trial( trialinfo(:,5) == 0, :, 132),2));
        %
        %      figure;
        %      imagesc(x)
        %         figure
        %         plot(1:251,mean(squeeze(mean(AV,2)),1),1:251, mean(squeeze(mean(V,2)),1))
        %
    end
end



for i_roi = 2
    roi_name = ROI_hf(i_roi).name;
    for i_vf = 1 : 2
        Zhat = [];
        dat = alldat(i_roi).vfield(i_vf).vchan_decode_data
        for iperm = 1 : nperms;
            %% Lets decode again, but now using a lower temporal resolution (to save space)
            % Reformatting data to timelock structure
            % cfg                = [];
            %cfg.resamplefs         = 30
            vchan_decode_data      = dat;% ft_resampledata(cfg,  dat );
            
            
            Y = vchan_decode_data.trial; %ROIdata.trial( indexes_selected(iper,:),:,:);     % Y: features x time x trials -> PIMS
            Y = permute(Y,[2 3 1]); % rearrange matrix as PIM %  Y: trials x features x  time -> PIMS
            
            % doing this to transform labels to 0 and 1 when labels are 1 and
            % 2.
            vchan_decode_data.trialinfo(:,4) = vchan_decode_data.trialinfo(:,4)  - min(vchan_decode_data.trialinfo(:,4))
            trialinfo = vchan_decode_data.trialinfo;
            
            labels = vchan_decode_data.trialinfo(:,4) > 0;
            X = labels';
            
            time = vchan_decode_data.time;
            label = vchan_decode_data.label;
            
            numF = size(Y, 1);
            numT = size(Y, 2);
            numN = size(Y, 3);
            
            
            % Create folds
            cfg = [];
            cfg.nFold = 5;
            folds = createFolds(cfg, X);
            
            cfg = [];
            cfg.folds = folds;
            cfg.feedback = 'yes';
            cfg.trainfun = 'train_array';
            cfg.traincfg.feedback = 'yes';
            cfg.traincfg.trainfun = 'train_LDA';
            cfg.traincfg.traincfg.gamma = 0.05;
            cfg.decodefun = 'decode_arrayGeneralization';
            cfg.decodecfg.feedback = 'yes';
            cfg.decodecfg.decodefun = 'decode_LDA';
            
            Zhat(iperm,:,:,:) =  decodeCrossValidation(cfg, X, Y);
            
        end
        
        Xhat = mean(Zhat, 1);
        predictions = (Xhat > 0);
        classes = squeeze(repmat(permute(X, [1 3 4 2]), [1, numT, numT, 1]));
        
        
        predictions = squeeze(predictions);
        
        % creating erf structure to save the diagonal
        
        struct_erf = {};
        struct_erf.time = vchan_decode_data.time; % update time information
        
        struct_erf.label = {[roi_name]};
        struct_erf.dimord = 'chan_time';
        ROI_decod(i_roi).vfield(i_vf).Xhat = Xhat;
        ROI_decod(i_roi).vfield(i_vf).classes = classes;
        ROI_decod(i_roi).vfield(i_vf).trialinfo = trialinfo;
        ROI_decod(i_roi).vfield(i_vf).name = roi_name;
        ROI_decod(i_roi).vfield(i_vf).temp_struct = struct_erf;
        
    end
end

save([G.outputPath, filesep, 'decod_2vf_visualROIS'], 'ROI_decod', '-v7.3');
%save([G.outputPath, filesep, 'decod_m_ROIS'], 'ROIm_decod', '-v7.3');


end
