text:
function tutorial_visual_single(bids_dir, reports_dir)
% TUTORIAL_VISUAL_SINGLE: Runs the Brainstorm/SPM group analysis pipeline (single subject, BIDS version).
%
% ONLINE TUTORIALS: https://neuroimage.usc.edu/brainstorm/Tutorials/VisualSingle
%
% INPUTS:
%    - bids_dir: Path to folder ds000117 (https://openneuro.org/datasets/ds000117)
%       |- derivatives/freesurfer/sub-XX                               : Segmentation folders generated with FreeSurfer
%       |- derivatives/meg_derivatives/sub-XX/ses-meg/meg/*.fif        : MEG+EEG recordings (processed with MaxFilter's SSS)
%       |- derivatives/meg_derivatives/sub-emptyroom/ses-meg/meg/*.fif : Empty room measurements
%    - reports_dir: If defined, exports all the reports as HTML to this folder

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Francois Tadel, Elizabeth Bock, 2016-2018

%% ===== SCRIPT VARIABLES =====
% Full list of subjects to process
SubjectNames = {'sub-01', 'sub-02', 'sub-03', 'sub-04', 'sub-05', 'sub-06', 'sub-07', 'sub-08', ...
                'sub-09', 'sub-10', 'sub-11', 'sub-12', 'sub-13', 'sub-14', 'sub-15', 'sub-16'};
% Empty-room dates for each subject (so that we can match automatically recordings with empty-room)
EmptyRoomSubj = 'sub-emptyroom';
AcquisitionDates = {'09-Apr-2009', '06-May-2009', '11-May-2009', '18-May-2009', '15-May-2009', '15-May-2009', '15-May-2009', '15-May-2009', ...
                    '15-May-2009', '15-May-2009', '01-Jun-2009', '01-Jun-2009', '01-Jun-2009', '26-Nov-2009', '08-Dec-2009', '08-Dec-2009'};
% Bad channels {iSubj} = {Run01, Run02, Run03, Run04, Run05, Run06}
BadChannels{1}  = {'EEG016', 'EEG070', 'EEG050',{'EEG008','EEG050'}, [], []};
BadChannels{2}  = {{'EEG027', 'EEG030', 'EEG038'}, 'EEG010', 'EEG010', 'EEG010', 'EEG010', 'EEG010'};
BadChannels{3}  = {{'EEG008','EEG017'}, {'EEG008','EEG017'}, {'EEG008','EEG017'}, {'EEG008','EEG017'}, {'EEG008','EEG017','EEG001'}, {'EEG008','EEG017','EEG020'}};
BadChannels{4}  = {{'EEG038'}, {'EEG038','EEG001','EEG016'}, {'EEG038','EEG001','EEG016'}, {'EEG038','EEG001'}, {'EEG038','EEG001','EEG016'}, {'EEG038','EEG001','EEG016'}};
BadChannels{5}  = {'EEG001', 'EEG001', [], [], [], []};
BadChannels{6}  = {'EEG068', [], 'EEG004', [], [], []};
BadChannels{7}  = {[], [], {'EEG004','EEG008'}, {'EEG004','EEG008'},{'EEG004','EEG008','EEG043','EEG045','EEG047'}, {'EEG004','EEG008'}};
BadChannels{8}  = {[], [], [], [], [], []};
BadChannels{9}  = {[], 'EEG004', 'EEG004', [], 'EEG004', 'EEG004'};
BadChannels{10} = {[], [], [], [], [], []};
BadChannels{11} = {{'EEG010','EEG050'}, 'EEG050', 'EEG050', 'EEG050', 'EEG050', 'EEG050'};
BadChannels{12} = {{'EEG024','EEG057'}, {'EEG024','EEG057'}, {'EEG024','EEG057'}, {'EEG024','EEG057','EEG070'}, {'EEG024','EEG057'}, {'EEG024','EEG057','EEG070'}};
BadChannels{13} = {'EEG009', 'EEG009', {'EEG009','EEG057','EEG69'}, 'EEG009', {'EEG009','EEG044'}, {'EEG009','EEG044'}};
BadChannels{14} = {'EEG029', 'EEG029', 'EEG029', {'EEG004','EEG008','EEG016','EEG029'}, {'EEG004','EEG008','EEG016','EEG029'}, {'EEG004','EEG008','EEG016','EEG029'}};
BadChannels{15} = {'EEG038', 'EEG038', 'EEG038', 'EEG038', {'EEG054','EEG038'}, 'EEG038'};
BadChannels{16} = {'EEG008', 'EEG008', 'EEG008', 'EEG008', 'EEG008', 'EEG008'};
% SSP components to remove {iSubj} = {sRun01, sRun02, sRun03, sRun03, sRun04, sRun05, sRun06},   sRun0X={ECG_GRAD,ECG_MAG}
SspSelect{1}  = {{1,1},         {1,1},   {1,1},   {1,1},   {1,1},   {1,1}};
SspSelect{2}  = {{1,1},         {1,1},   {1,1},   {1,1},   {3,1},   {1,1}};
SspSelect{3}  = {{[],1},        {1,1},   {1,1},   {1,1},   {1,1},   {1,1}};
SspSelect{4}  = {{[],1},        {[],1},  {[],1},  {[],1},  {[],1},  {[],1}};
SspSelect{5}  = {{2,1},         {1,1},   {1,1},   {[],1},  {1,1},   {1,1}};
SspSelect{6}  = {{2,1},         {2,1},   {1,1},   {2,1},   {1,1},   {2,1}};
SspSelect{7}  = {{1,1},         {1,1},   {1,1},   {1,1},   {1,1},   {1,1}};
SspSelect{8}  = {{1,1},         {1,1},   {[],1},  {2,1},   {1,1},   {2,1}};
SspSelect{9}  = {{1,1},         {1,1},   {[],1},  {[],1},  {1,1},   {1,1}};
SspSelect{10} = {{1,1},         {1,1},   {1,1},   {1,1},   {1,1},   {1,1}};
SspSelect{11} = {{[],1},        {[],1},  {[],1},  {[],1},  {[],[]}, {[],[]}};
SspSelect{12} = {{[1,2],[1,2]}, {1,1},   {1,1},   {1,1},   {1,1},   {1,1}};
SspSelect{13} = {{[],[]},       {[],[]}, {[],[]}, {[],[]}, {[],[]}, {[],[]}};
SspSelect{14} = {{1,1},         {1,1},   {1,1},   {1,1},   {1,1},   {1,1}};
SspSelect{15} = {{1,1},         {1,1},   {1,1},   {1,1},   {1,1},   {1,1}};
SspSelect{16} = {{1,1},         {1,1},   {1,1},   {2,1},   {1,1},   {1,1}};

%% ===== CREATE PROTOCOL =====
% Start brainstorm without the GUI
if ~brainstorm('status')
    brainstorm nogui
end
% Output folder for reports
if (nargin < 2) || isempty(reports_dir) || ~isdir(reports_dir)
    reports_dir = [];
end
% You have to specify the folder in which the tutorial dataset is unzipped
if (nargin < 1) || isempty(bids_dir) || ~file_exist(bids_dir) || ~file_exist(bst_fullfile(bids_dir, 'derivatives')) || ~file_exist(bst_fullfile(bids_dir, 'dataset_description.json'))
    error('The first argument must be the full path to the tutorial folder.');
end
% The protocol name has to be a valid folder name (no spaces, no weird characters...)
ProtocolName = 'TutorialVisual';
% Delete existing protocol
gui_brainstorm('DeleteProtocol', ProtocolName);
% Create new protocol
gui_brainstorm('CreateProtocol', ProtocolName, 0, 0);
% Set visualization filters: 40Hz low-pass, no high-pass
panel_filter('SetFilters', 1, 40, 0, [], 0, [], 0, 0);
% Set colormap: local color scale
bst_colormaps('SetMaxMode', 'meg', 'local');
bst_colormaps('SetMaxMode', 'eeg', 'local');

%% ===== PRE-PROCESS AND IMPORT =====
for iSubj = 1:16
    % Start a new report (one report per subject)
    bst_report('Start');
    disp(sprintf('\n===== IMPORT: SUBJECT #%d =====\n', iSubj));

    % If subject already exists: delete it
    [sSubject, iSubject] = bst_get('Subject', SubjectNames{iSubj});
    if ~isempty(sSubject)
        db_delete_subjects(iSubject);
    end

    % ===== FILES TO IMPORT =====
    % Build the path of the files to import
    AnatDir    = fullfile(bids_dir, 'derivatives', 'freesurfer', SubjectNames{iSubj}, 'ses-mri', 'anat');
    DataDir    = fullfile(bids_dir, 'derivatives',  'meg_derivatives', SubjectNames{iSubj}, 'ses-meg', 'meg');
    % Check if the folder contains the required files
    if ~file_exist(AnatDir)
        error(['The folder "' AnatDir '" does not exist.']);
    end
    if ~file_exist(DataDir)
        error(['The folder "' DataDir '" does not exist.']);
    end

    % ===== ANATOMY =====
    % Process: Import anatomy folder
    bst_process('CallProcess', 'process_import_anatomy', [], [], ...
        'subjectname', SubjectNames{iSubj}, ...
        'mrifile',     {AnatDir, 'FreeSurfer'}, ...
        'nvertices',   15000);

    % ===== PROCESS EACH RUN =====
    for iRun = 1:6
        % Files to import
        FifFile = bst_fullfile(DataDir, sprintf('%s_ses-meg_task-facerecognition_run-%02d_proc-sss_meg.fif', SubjectNames{iSubj}, iRun));

        % ===== LINK CONTINUOUS FILE =====
        % Process: Create link to raw file
        sFileRaw = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
            'subjectname',    SubjectNames{iSubj}, ...
            'datafile',       {FifFile, 'FIF'}, ...
            'channelreplace', 1, ...
            'channelalign',   0);
        % Set acquisition date
        panel_record('SetAcquisitionDate', sFileRaw.iStudy, AcquisitionDates{iSubj});

        % ===== PREPARE CHANNEL FILE =====
        % Process: Set channels type
        bst_process('CallProcess', 'process_channel_settype', sFileRaw, [], ...
            'sensortypes', 'EEG061, EEG064', ...
            'newtype',     'NOSIG');
        bst_process('CallProcess', 'process_channel_settype', sFileRaw, [], ...
            'sensortypes', 'EEG062', ...
            'newtype',     'EOG');
        bst_process('CallProcess', 'process_channel_settype', sFileRaw, [], ...
            'sensortypes', 'EEG063', ...
            'newtype',     'ECG');

        % Process: Remove head points
        sFileRaw = bst_process('CallProcess', 'process_headpoints_remove', sFileRaw, [], ...
            'zlimit', 0);
        % Process: Refine registration
        sFileRaw = bst_process('CallProcess', 'process_headpoints_refine', sFileRaw, []);
        % Process: Project electrodes on scalp
        sFileRaw = bst_process('CallProcess', 'process_channel_project', sFileRaw, []);

        % Process: Snapshot: Sensors/MRI registration
        bst_process('CallProcess', 'process_snapshot', sFileRaw, [], ...
            'target',   1, ...  % Sensors/MRI registration
            'modality', 1, ...  % MEG (All)
            'orient',   1, ...  % left
            'Comment',  sprintf('MEG/MRI Registration: Subject #%d, Run #%d', iSubj, iRun));
        bst_process('CallProcess', 'process_snapshot', sFileRaw, [], ...
            'target',   1, ...  % Sensors/MRI registration
            'modality', 4, ...  % EEG
            'orient',   1, ...  % left
            'Comment',  sprintf('EEG/MRI Registration: Subject #%d, Run #%d', iSubj, iRun));

        % ===== IMPORT TRIGGERS =====
        % Process: Read from channel
        bst_process('CallProcess', 'process_evt_read', sFileRaw, [], ...
            'stimchan',  'STI101', ...
            'trackmode', 2, ...  % Bit: detect the changes for each bit independently
            'zero',      0);
        % Process: Group by name
        bst_process('CallProcess', 'process_evt_groupname', sFileRaw, [], ...
            'combine', 'Unfamiliar=3,4', ...
            'dt',      0, ...
            'delete',  1);
        % Process: Rename event
        bst_process('CallProcess', 'process_evt_rename', sFileRaw, [], ...
            'src',  '3', ...
            'dest', 'Famous');
        % Process: Rename event
        bst_process('CallProcess', 'process_evt_rename', sFileRaw, [], ...
            'src',  '5', ...
            'dest', 'Scrambled');
        % Process: Add time offset
        bst_process('CallProcess', 'process_evt_timeoffset', sFileRaw, [], ...
            'info',      [], ...
            'eventname', 'Famous, Unfamiliar, Scrambled', ...
            'offset',    0.0345);
        % Process: Delete events
        bst_process('CallProcess', 'process_evt_delete', sFileRaw, [], ...
            'eventname', '1,2,6,7,8,9,10,11,12,13,14,15,16');
        % Process: Detect cHPI activity (Elekta):STI201
        bst_process('CallProcess', 'process_evt_detect_chpi', sFileRaw, [], ...
            'eventname',   'chpi_bad', ...
            'channelname', 'STI201', ...
            'method',      'off');  % Mark as bad when the HPI coils are OFF

        % ===== FREQUENCY FILTERS =====
        % Process: Notch filter: 50Hz 100Hz 150Hz 200Hz
        sFileClean = bst_process('CallProcess', 'process_notch', sFileRaw, [], ...
            'freqlist',    [50, 100, 150, 200], ...
            'sensortypes', 'MEG, EEG', ...
            'read_all',    0);
        % Process: Power spectrum density (Welch)
        sFilesPsd = bst_process('CallProcess', 'process_psd', [sFileRaw, sFileClean], [], ...
            'timewindow',  [], ...
            'win_length',  4, ...
            'win_overlap', 50, ...
            'sensortypes', 'MEG, EEG', ...
            'edit', struct(...
                 'Comment',         'Power', ...
                 'TimeBands',       [], ...
                 'Freqs',           [], ...
                 'ClusterFuncTime', 'none', ...
                 'Measure',         'power', ...
                 'Output',          'all', ...
                 'SaveKernel',      0));
        % Process: Snapshot: Frequency spectrum
        bst_process('CallProcess', 'process_snapshot', sFilesPsd, [], ...
            'target',   10, ...  % Frequency spectrum
            'Comment',  sprintf('Power spctrum: Subject #%d, Run #%d', iSubj, iRun));

        % ===== BAD CHANNELS =====
        if ~isempty(BadChannels{iSubj}{iRun})
            % Process: Set bad channels
            bst_process('CallProcess', 'process_channel_setbad', sFileClean, [], ...
                'sensortypes', BadChannels{iSubj}{iRun});
        end

        % ===== EEG REFERENCE =====
        % Process: Re-reference EEG
        bst_process('CallProcess', 'process_eegref', sFileClean, [], ...
            'eegref',      'AVERAGE', ...
            'sensortypes', 'EEG');

        % ===== DETECT ARTIFACTS ======
        % Process: Detect heartbeats
        bst_process('CallProcess', 'process_evt_detect_ecg', sFileClean, [], ...
            'channelname', 'EEG063', ...
            'timewindow',  [], ...
            'eventname',   'cardiac');
        % Different amplitude thresholds for different subjects
        if strcmpi(SubjectNames{iSubj}, 'sub-05')
            thresholdMAX = 50;
        else
            thresholdMAX = 100;
        end
        % Process: Detect: blink_BAD    -   Detects all events where the amplitude exceeds 100uV
        bst_process('CallProcess', 'process_evt_detect_threshold', sFileClean, [], ...
            'eventname',    'blink_BAD', ...
            'channelname',  'EEG062', ...
            'timewindow',   [], ...
            'thresholdMAX', thresholdMAX, ...
            'units',        3, ...  % uV (10^-6)
            'bandpass',     [0.3, 20], ...
            'isAbsolute',   1, ...
            'isDCremove',   0);

        % ===== SSP COMPUTATION =====
        % Process: SSP ECG: cardiac
        bst_process('CallProcess', 'process_ssp_ecg', sFileClean, [], ...
            'eventname',   'cardiac', ...
            'sensortypes', 'MEG GRAD', ...
            'usessp',      1, ...
            'select',      SspSelect{iSubj}{iRun}{1});
        bst_process('CallProcess', 'process_ssp_ecg', sFileClean, [], ...
            'eventname',   'cardiac', ...
            'sensortypes', 'MEG MAG', ...
            'usessp',      1, ...
            'select',      SspSelect{iSubj}{iRun}{2});
        % Process: Snapshot: SSP projectors
        bst_process('CallProcess', 'process_snapshot', sFileClean, [], ...
            'target',   2, ...
            'Comment',  sprintf('Subject #%d, Run #%d', iSubj, iRun));   % SSP projectors

        % ===== IMPORT BAD EVENTS =====
        % Get bad segments: this is typically done manually, not from a script
        BadSegments = GetBadSegments(iSubj, iRun);
        % Process: Import from file
        bst_process('CallProcess', 'process_evt_import', sFileClean, [], ...
            'evtfile', {BadSegments, 'ARRAY-TIMES'}, ...
            'evtname', 'BAD');

        % ===== IMPORT TRIALS =====
        % Process: Import MEG/EEG: Events
        sFilesEpochs = bst_process('CallProcess', 'process_import_data_event', sFileClean, [], ...
            'subjectname', SubjectNames{iSubj}, ...
            'condition',   '', ...
            'eventname',   'Famous, Scrambled, Unfamiliar', ...
            'timewindow',  [], ...
            'epochtime',   [-0.5, 1.2], ...
            'createcond',  0, ...
            'ignoreshort', 1, ...
            'usectfcomp',  1, ...
            'usessp',      1, ...
            'freq',        [], ...
            'baseline',    [-0.5, -0.0009]);

        % ===== AVERAGE: RUN =====
        % Process: Average: By trial group (folder average)
        sFilesAvg = bst_process('CallProcess', 'process_average', sFilesEpochs, [], ...
            'avgtype',    5, ...  % By trial group (folder average)
            'avg_func',   1, ...  % Arithmetic average:  mean(x)
            'weighted',   0, ...
            'keepevents', 0);
        % Process: Snapshot: Recordings time series
        bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
            'target',   5, ...  % Recordings time series
            'modality', 4, ...  % EEG
            'time',     0.11, ...
            'Comment',  sprintf('Subject #%d, Run #%d', iSubj, iRun));
        % Process: Snapshot: Recordings topography
        bst_process('CallProcess', 'process_snapshot', sFilesAvg, [], ...
            'target',   6, ...  % Recordings topography (one time)
            'modality', 4, ...  % EEG
            'time',     0.11, ...
            'Comment',  sprintf('Subject #%d, Run #%d', iSubj, iRun));

        % ===== COMPUTE NOISECOV: EEG =====
        % Process: Compute covariance (noise or data)
        bst_process('CallProcess', 'process_noisecov', sFilesEpochs, [], ...
            'baseline',       [-0.5, -0.0009], ...
            'sensortypes',    'EEG', ...
            'target',         1, ...  % Noise covariance     (covariance over baseline time window)
            'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
            'identity',       0, ...
            'copycond',       0, ...
            'copysubj',       0, ...
            'replacefile',    1);  % Replace
    end

    % Save report
    ReportFile = bst_report('Save', []);
    if ~isempty(reports_dir) && ~isempty(ReportFile)
        bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_' SubjectNames{iSubj} '.html']));
    end
end

%% ===== EMPTY ROOM RECORDINGS =====
disp(sprintf('\n===== IMPORT: EMPTY-ROOM =====\n'));
% Loop on all the noise sessions
NoiseFiles = {};
for ses = {'20090409', '20090506', '20090511', '20090515', '20090518', '20090601', '20091126', '20091208'}
    NoiseFiles{end+1} = fullfile(bids_dir, 'derivatives', 'meg_derivatives', EmptyRoomSubj, ['ses-' ses{1}], 'meg', ['sub-emptyroom_ses-' ses{1} '_task-noise_proc-sss_meg.fif']);
end
% Process: Create link to raw file
sFilesNoise = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
    'subjectname',    EmptyRoomSubj, ...
    'datafile',       {NoiseFiles, 'FIF'}, ...
    'channelreplace', 1, ...
    'channelalign',   0);
% Process: Notch filter: 50Hz 100Hz 150Hz 200Hz
sFileNoiseClean = bst_process('CallProcess', 'process_notch', sFilesNoise, [], ...
        'freqlist',    [50, 100, 150, 200], ...
        'sensortypes', 'MEG, EEG', ...
        'read_all',    0);
% Process: Compute noise covariance
bst_process('CallProcess', 'process_noisecov', sFileNoiseClean, [], ...
    'baseline',    [], ...
    'sensortypes', 'MEG', ...
    'target',      1, ...  % Noise covariance     (covariance over baseline time window)
    'dcoffset',    1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',    0, ...
    'copycond',    1, ...
    'copysubj',    1, ...
    'copymatch',   1, ...
    'replacefile', 2);  % Merge

%% ===== SOURCE ESTIMATION =====
% Start a new report (one report for the source estimation of all the subjects)
bst_report('Start');
% Loop on the subjects: This loop is separated from the previous one, because we should
% compute the BEM surfaces after importing all the runs, so that the registration is done
% using the high resolution head surface, instead of the smooth scalp BEM layer.
for iSubj = 1:length(SubjectNames)
    disp(sprintf('\n===== SOURCES: SUBJECT #%d =====\n', iSubj));

    % ===== BEM SURFACES =====
    % Process: Generate BEM surfaces
    bst_process('CallProcess', 'process_generate_bem', [], [], ...
        'subjectname', SubjectNames{iSubj}, ...
        'nscalp',      1082, ...
        'nouter',      642, ...
        'ninner',      642, ...
        'thickness',   4, ...
        'method',      'brainstorm');

    % ===== SELECT ALL AVERAGES =====
    % Process: Select data files in: */*
    sFilesAvg = bst_process('CallProcess', 'process_select_files_data', [], [], ...
        'subjectname',   SubjectNames{iSubj});
    % Process: Select file comments with tag: Avg
    sFilesAvg = bst_process('CallProcess', 'process_select_tag', sFilesAvg, [], ...
        'tag',    'Avg');  % Select only the files with the tag

    % ===== COMPUTE HEAD MODELS =====
    % Process: Compute head model (only for the first run of the subject)
    bst_process('CallProcess', 'process_headmodel', sFilesAvg(1), [], ...
        'sourcespace', 1, ...  % Cortex surface
        'meg',         3, ...  % Overlapping spheres
        'eeg',         3, ...  % OpenMEEG BEM
        'ecog',        1, ...  %
        'seeg',        1, ...  %
        'openmeeg',    struct(...
             'BemSelect',    [1, 1, 1], ...
             'BemCond',      [1, 0.0125, 1], ...
             'BemNames',     {{'Scalp', 'Skull', 'Brain'}}, ...
             'BemFiles',     {{}}, ...
             'isAdjoint',    0, ...
             'isAdaptative', 1, ...
             'isSplit',      0, ...
             'SplitLength',  4000));
    % Get all the runs for this subject (ie the list of the study indices)
    iStudyOther = setdiff(unique([sFilesAvg.iStudy]), sFilesAvg(1).iStudy);
    % Copy the forward model file to the other runs
    sHeadmodel = bst_get('HeadModelForStudy', sFilesAvg(1).iStudy);
    for iStudy = iStudyOther
        db_add(iStudy, sHeadmodel.FileName);
    end

    % ===== COMPUTE SOURCES: MEG =====
    % Process: Compute sources [2018]
    sAvgSrcMeg = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
        'output',  1, ...  % Kernel only: shared
        'inverse', struct(...
             'Comment',        'MN: MEG ALL', ...
             'InverseMethod',  'minnorm', ...
             'InverseMeasure', 'amplitude', ...
             'SourceOrient',   {{'fixed'}}, ...
             'Loose',          0.2, ...
             'UseDepth',       1, ...
             'WeightExp',      0.5, ...
             'WeightLimit',    10, ...
             'NoiseMethod',    'reg', ...
             'NoiseReg',       0.1, ...
             'SnrMethod',      'fixed', ...
             'SnrRms',         1e-06, ...
             'SnrFixed',       3, ...
             'ComputeKernel',  1, ...
             'DataTypes',      {{'MEG GRAD', 'MEG MAG'}}));
    % Process: Snapshot: Sources (one time)   -   Loop only to get a correct comment for the report
    for i = 1:length(sAvgSrcMeg)
        bst_process('CallProcess', 'process_snapshot', sAvgSrcMeg(i), [], ...
            'target',         8, ...  % Sources (one time)
            'orient',         4, ...  % bottom
            'time',           0.11, ...
            'threshold',      20, ...
            'Comment',        ['MEG sources: ' sFilesAvg(i).FileName]);
    end

    % ===== COMPUTE SOURCES: EEG =====
    % Process: Compute sources [2018]
    sAvgSrcEeg = bst_process('CallProcess', 'process_inverse_2018', sFilesAvg, [], ...
        'output',  1, ...  % Kernel only: shared
        'inverse', struct(...
             'Comment',        'MN: EEG', ...
             'InverseMethod',  'minnorm', ...
             'InverseMeasure', 'amplitude', ...
             'SourceOrient',   {{'fixed'}}, ...
             'Loose',          0.2, ...
             'UseDepth',       1, ...
             'WeightExp',      0.5, ...
             'WeightLimit',    10, ...
             'NoiseMethod',    'reg', ...
             'NoiseReg',       0.1, ...
             'SnrMethod',      'fixed', ...
             'SnrRms',         1e-06, ...
             'SnrFixed',       3, ...
             'ComputeKernel',  1, ...
             'DataTypes',      {{'EEG'}}));
    % Process: Snapshot: Sources (one time)   -   Loop only to get a correct comment for the report
    for i = 1:length(sAvgSrcEeg)
        bst_process('CallProcess', 'process_snapshot', sAvgSrcEeg(i), [], ...
            'target',         8, ...  % Sources (one time)
            'orient',         4, ...  % bottom
            'time',           0.11, ...
            'threshold',      10, ...
            'Comment',        ['EEG sources: ' sFilesAvg(i).FileName]);
    end
end
% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_sources.html']));
end

%% ===== TIME-FREQUENCY =====
% Start a new report (one report for the time-frequency of all the subjects)
bst_report('Start');
% List of conditions to process separately
AllConditions = {'Famous', 'Scrambled', 'Unfamiliar'};
% Channels to display in the screen capture, by order of preference (if the first channel is bad, use the following)
SelChannel = {'EEG070','EEG060','EEG065','EEG050','EEG003'};
% Compute one separate time-frequency average for each subject/run/condition
for iSubj = 1:length(SubjectNames)
    disp(sprintf('\n===== TIME-FREQUENCY: SUBJECT #%d =====\n', iSubj));
    for iRun = 1:6
        % Process: Select data files in: Subject/Run
        sTrialsAll = bst_process('CallProcess', 'process_select_files_data', [], [], ...
            'subjectname',   SubjectNames{iSubj}, ...
            'condition',     sprintf('sub-%02d_ses-meg_task-facerecognition_run-%02d_proc-sss_meg_notch', iSubj, iRun));
        % Loop on the conditions
        for iCond = 1:length(AllConditions)
            % Comment describing this average
            strComment = [SubjectNames{iSubj}, ' / ', sprintf('run_%02d', iRun), ' / ', AllConditions{iCond}];
            disp(['BST> ' strComment]);
            % Find the first good channel in the display list
            if isempty(BadChannels{iSubj}{iRun})
                iSel = 1;
            else
                iSel = find(~ismember(SelChannel,BadChannels{iSubj}{iRun}), 1);
            end
            % Process: Select file comments with tag: Avg
            sTrialsCond = bst_process('CallProcess', 'process_select_tag', sTrialsAll, [], ...
                'tag',    [AllConditions{iCond}, '_trial'], ...
                'search', 1, ...  % Search the file names
                'select', 1);  % Select only the files with the tag
            % Process: Time-frequency (Morlet wavelets), averaged across trials
            sTimefreq = bst_process('CallProcess', 'process_timefreq', sTrialsCond, [], ...
                'sensortypes', 'MEG MAG, EEG', ...
                'edit',        struct(...
                     'Comment',         ['Avg: ' AllConditions{iCond} ', Power, 6-60Hz'], ...
                     'TimeBands',       [], ...
                     'Freqs',           [6, 6.8, 7.6, 8.6, 9.7, 11, 12.4, 14, 15.8, 17.9, 20.2, 22.8, 25.7, 29, 32.7, 37, 41.7, 47.1, 53.2, 60], ...
                     'MorletFc',        1, ...
                     'MorletFwhmTc',    3, ...
                     'ClusterFuncTime', 'none', ...
                     'Measure',         'power', ...
                     'Output',          'average', ...
                     'RemoveEvoked',    0, ...
                     'SaveKernel',      0), ...
                'normalize',   'none');  % None: Save non-standardized time-frequency maps
            % Process: Extract time: [-200ms,900ms]
            sTimefreq = bst_process('CallProcess', 'process_extract_time', sTimefreq, [], ...
                'timewindow', [-0.2, 0.9], ...
                'overwrite',  1);
            % Screen capture of one sensor
            hFigTf = view_timefreq(sTimefreq.FileName, 'SingleSensor', SelChannel{iSel});
            bst_report('Snapshot', hFigTf, strComment, 'Time-frequency', [200, 200, 400, 250]);
            close(hFigTf);
        end
    end
end
% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_timefreq.html']));
end
end

%% ===== SUPPORT FUNCTIONS =====
function BadSeg = GetBadSegments(iSubj, iRun)
     BadSegments{1} = {...
         [247.867 248.185; 598.999 598.999; 598.999 598.999; 611.999 611.999; 612.999 612.999; 613.999 613.999; 616.999 616.999; 617.999 617.999; 623.999 623.999; 715.209 715.467], ...
         [84.791 85.166], ...
         [79.183 80.167], ...
         [64.309 65.185], ...
         [90.958 91.167; 178.005 178.355; 293.282 295.919; 312.298 316.479; 353.835 357.716], ...
         [60.292 66.802; 69.975 71.210; 105.233 107.586; 108.822 109.506; 376.225 376.325]};
     BadSegments{2} = {...
         [223.806 224.199; 279.772 279.895; 453.241 455.108; 692.423 692.593], ...
         [65.298 66.194; 304.727 306.178; 399.165 400.732], ...
         [203.141 205.085; 281.579 287.883; 420.395 421.128], ...
         [387.118 388.229; 440.318 441.900; 554.825 558.744], ...
         [71.000 80.999; 82.750 87.367; 149.528 149.667; 264.747 267.995; 368.415 371.973; 376.263 378.763; 398.334 401.551; 537.410 541.645], ...
         [38.000 47.999; 47.825 50.046; 61.298 61.384; 249.653 253.379; 282.917 283.820; 286.135 287.616; 298.167 300.196; 328.254 329.511; 335.957 337.817; 478.277 480.707]};
     BadSegments{3} = {...
         [406.312 407.207; 727.055 728.714], ...
         [84.894 85.156; 152.028 152.946; 297.835 298.915; 418.272 421.845; 554.084 554.794], ...
         [73.758 74.159; 378.212 378.536; 406.065 407.099; 470.541 471.698; 488.900 491.168; 529.596 530.453], ...
         [94.874 95.152; 317.385 321.374; 325.696 327.055; 439.220 439.829; 454.473 455.175; 486.196 486.829; 518.660 522.015; 524.400 525.249; 562.417 570.325], ...
         [96.208 97.181; 98.942 99.096; 135.005 135.754; 143.990 144.599; 250.139 250.247; 300.459 300.559; 338.265 339.322; 545.913 546.067], ...
         [91.415 92.156; 284.843 286.525; 297.886 298.404; 317.046 317.163; 332.698 332.791; 358.946 359.402; 428.405 428.775; 478.374 478.690; 549.866 550.128]};
     BadSegments{4} = {...
         [22.967 22.967; 50.036 50.098; 52.058 52.058; 156.653 156.653; 171.565 173.386; 239.544 242.105; 268.162 270.175; 268.992 268.992; 316.032 316.032; 338.283 339.000; 357.959 361.909; 370.871 370.871; 381.579 383.677; 437.731 437.731; 463.482 468.505; 476.135 479.838; 486.652 488.272; 504.860 508.999], ...
         [309.493 311.707; 342.681 344.525; 354.019 357.321; 390.023 391.225; 393.926 395.855; 404.221 405.069; 432.522 435.932; 459.048 460.715; 471.763 478.529; 549.387 551.999; 591.087 594.143; 608.541 611.079; 624.847 626.615; 649.648 651.570], ...
         [57.411 58.198; 88.346 88.955; 200.761 202.335; 227.016 227.688; 257.726 258.054; 356.798 359.005; 404.260 411.003], ...
         [46.000 54.823; 61.000 70.332; 203.005 207.125; 275.875 278.121; 313.500 314.824; 337.973 338.636; 422.505 426.239], ...
         [58.000 62.479; 78.250 85.166; 89.955 91.360; 116.322 117.888; 130.013 131.987; 149.509 150.489; 174.650 175.823; 182.030 183.334; 196.758 197.384; 204.458 204.697; 205.236 208.663; 311.028 316.383; 320.700 327.181; 332.437 335.354; 344.205 346.133; 374.208 374.865; 385.519 386.214; 441.942 444.241; 453.957 456.997; 486.039 487.004; 501.238 504.185; 512.962 514.675; 553.398 556.215], ...
         [41.406 45.743; 58.681 59.144; 108.086 108.896; 140.633 143.750; 196.110 199.474; 210.778 210.971; 234.649 235.143; 258.081 259.632; 339.101 340.805; 390.277 390.609; 438.935 442.122; 528.221 534.031]};
     BadSegments{5} = {...
         [265.539 265.778; 266.334 266.495; 268.479 268.965; 367.428 367.636; 439.655 442.779; 453.497 453.853; 504.997 505.329; 519.513 519.683; 595.674 595.982; 602.000 602.463], ...
         [121.113 121.499; 124.971 126.213; 253.735 254.075; 272.232 272.464; 272.895 273.104; 346.368 346.645; 368.812 369.052; 406.382 406.605; 452.920 453.113; 454.903 455.112; 507.655 507.840; 508.766 509.013; 584.853 585.030; 594.831 595.656; 597.261 602.249], ...
         [37.251 37.497; 38.825 39.056; 40.615 41.849; 43.624 44.758; 53.333 53.641; 54.698 55.076; 57.668 59.196; 79.129 79.360; 81.475 81.714; 122.658 123.375; 284.787 285.296; 288.754 288.993; 345.790 346.022; 421.212 421.459; 481.428 482.207; 503.408 503.685; 504.272 504.449; 524.451 524.714; 526.913 531.499], ...
         [87.322 88.178; 91.085 91.325; 95.121 95.491; 114.174 114.397; 129.874 130.113; 151.220 151.544; 281.689 281.959; 532.966 533.345], ...
         [59.176 60.218; 74.854 75.317; 308.180 309.877; 380.705 381.059], ...
         [182.382 183.245; 196.220 196.736; 276.018 276.327; 292.490 294.086; 370.755 370.847; 435.644 436.624; 467.535 468.460; 522.838 525.847]};
     BadSegments{6} = {...
         [141.690 142.424; 157.070 157.417; 355.138 356.025; 423.999 423.999; 424.999 424.999; 426.999 426.999; 427.999 427.999; 486.430 488.151; 493.999 493.999; 501.511 501.619; 501.549 501.549; 501.585 501.585; 502.999 502.999; 503.999 503.999; 540.999 540.999; 541.999 541.999; 555.999 555.999; 556.999 556.999; 561.999 561.999; 563.999 563.999; 564.999 564.999; 565.999 565.999; 567.999 567.999], ...
         [64.898 65.161; 71.700 72.718; 226.185 226.740; 324.124 324.425; 329.062 329.301; 486.143 486.975], ...
         [62.300 63.148; 266.254 266.639; 409.920 410.221], ...
         [54.048 55.214; 330.893 331.255], ...
         [185.616 186.495; 331.411 331.796; 386.843 387.028; 387.999 387.999; 389.999 389.999; 434.575 434.875; 519.802 519.995], ...
         [44.720 45.167; 211.446 211.964; 368.955 369.172]};
     BadSegments{7} = {...
         [154.966 155.144; 551.639 551.855], ...
         [88.774 89.167; 107.999 107.999; 109.999 109.999; 110.999 110.999; 112.999 112.999; 113.999 113.999; 114.999 114.999; 119.999 119.999; 121.999 121.999; 124.223 125.465; 137.011 138.871], ...
         [81.627 82.136; 377.953 381.046], ...
         [241.136 242.171; 543.849 544.196; 596.639 598.553; 600.227 601.075; 603.999 603.999; 605.999 605.999; 609.999 609.999; 611.305 612.809; 614.999 614.999; 615.803 616.937; 623.694 625.877; 653.999 653.999; 655.055 655.756; 663.999 663.999], ...
         [68.852 69.481; 74.034 75.200; 78.426 78.600; 104.497 105.963; 253.951 254.034; 256.964 257.038; 257.915 258.048; 323.254 324.156; 365.880 368.131; 369.952 370.060; 371.728 372.999; 430.931 431.965; 535.521 542.999], ...
         [70.271 71.205; 94.441 96.445; 98.613 99.126; 112.318 112.749; 131.686 132.265; 148.935 150.615; 161.120 161.600; 205.325 208.214; 215.035 215.863; 217.403 218.818; 286.178 287.171; 399.075 404.853]};
     BadSegments{8} = {...
         [238.505 238.546; 256.354 257.224; 316.116 316.167; 341.519 341.558; 356.493 356.566; 380.095 380.170; 391.906 392.048; 448.457 448.562; 469.931 470.028; 530.488 530.575; 555.180 558.136; 562.152 562.245; 588.403 588.502; 625.205 625.662; 638.019 638.438; 649.982 650.008; 650.691 651.357; 651.925 652.061; 665.445 665.472; 695.923 696.015; 706.528 706.720; 729.706 732.534], ...
         [99.546 106.114; 113.245 114.199; 150.885 154.989; 277.486 278.075; 333.645 335.574; 339.358 340.646; 371.860 375.394; 497.459 499.156], ...
         [49.008 50.205; 231.446 233.406; 329.590 329.659; 355.101 356.019; 360.733 360.973; 372.955 374.891; 389.283 392.068; 453.610 455.431; 464.632 465.265; 489.209 489.996; 514.777 515.295], ...
         [178.368 179.232; 293.865 294.521; 418.252 418.314; 450.124 450.209; 480.236 480.445; 492.725 493.975; 495.170 497.624; 500.382 500.459; 504.247 504.402; 628.017 628.079; 628.827 628.905; 630.209 630.332], ...
         [100.616 101.172; 107.691 108.539; 188.814 188.875; 193.119 193.281; 207.964 208.033; 432.254 432.385; 489.834 489.911; 517.960 518.037; 518.955 519.040; 520.938 521.062; 521.945 522.037; 523.959 524.052; 525.942 526.057; 526.945 527.015; 528.958 529.059; 531.138 531.192; 531.979 532.048; 532.951 533.028; 533.962 534.024], ...
         [135.997 137.232; 383.565 383.657; 418.763 418.955]};
     BadSegments{9} = {...
         [215.107 216.187; 262.388 262.388; 287.519 287.635; 289.895 290.135; 311.999 311.999; 350.161 351.179; 526.154 527.033; 564.999 564.999; 584.935 585.059; 587.999 587.999; 601.999 601.999; 603.999 603.999; 608.999 608.999; 612.999 612.999], ...
         [47.667 47.775; 49.172 54.681; 67.195 70.805; 78.988 80.477; 138.701 138.948; 138.727 138.727; 138.728 138.728; 138.728 138.728; 138.734 138.734; 138.734 138.734; 138.881 138.881; 138.881 138.881; 138.907 138.907; 140.993 141.093; 141.037 141.037; 141.045 141.045; 155.795 155.864; 155.822 155.822; 155.849 155.849; 168.999 168.999; 206.999 206.999; 219.999 219.999; 225.999 225.999; 226.999 226.999; 228.999 228.999; 236.999 236.999; 242.463 242.463; 242.463 242.463; 242.487 242.487; 247.999 247.999; 251.999 251.999; 252.483 252.483; 253.315 254.163; 265.999 265.999; 267.999 267.999; 272.358 272.358; 272.410 272.410; 298.999 298.999; 300.999 300.999; 314.037 314.037; 314.047 314.047; 321.813 321.813; 321.813 321.813; 321.833 321.833; 329.117 329.117; 329.117 329.117; 329.156 329.156; 346.999 346.999; 347.999 347.999; 349.320 349.320; 349.329 349.329; 352.528 355.869; 364.040 364.040; 396.278 397.420; 404.865 404.865; 407.905 407.905; 407.905 407.905; 407.954 407.954; 418.454 418.454; 418.454 418.454; 418.486 418.486; 441.999 441.999; 444.999 444.999; 447.999 447.999; 453.550 453.650; 454.931 455.055; 457.999 457.999; 479.964 480.079; 481.999 481.999; 482.949 483.073; 488.948 489.064; 511.999 511.999; 520.999 520.999; 523.999 523.999; 526.954 527.069; 533.999 533.999; 537.999 537.999], ...
         [82.889 83.198; 206.143 207.424; 403.873 404.096; 406.790 407.485; 413.936 414.075; 415.912 416.112; 536.621 537.532], ...
         [182.999 182.999; 182.999 182.999; 183.999 183.999; 195.999 195.999; 208.999 208.999; 209.999 209.999; 278.601 278.955; 413.913 414.067; 415.250 417.792; 419.940 420.087; 420.999 420.999; 512.999 512.999; 514.363 516.254; 521.917 522.134; 522.999 522.999; 523.915 524.101; 533.999 533.999; 538.999 538.999; 539.892 540.062; 543.999 543.999; 548.874 549.089; 549.900 550.062; 552.755 554.075; 585.999 585.999; 587.957 588.073; 588.999 588.999; 594.999 594.999; 603.971 604.056; 604.928 605.097; 615.986 618.495; 623.999 623.999], ...
         [53.215 54.203; 164.227 168.965; 201.433 202.775; 292.973 295.889; 303.961 304.817; 309.354 311.236; 313.123 313.639; 322.547 323.017; 331.141 334.852; 356.637 356.985; 367.855 368.079; 369.995 373.459; 377.849 378.142; 406.072 407.306; 436.164 436.665; 458.033 458.905; 516.889 518.656; 517.812 518.684], ...
         [113.945 115.225; 198.570 198.863; 264.162 264.795; 383.705 385.641; 396.477 397.064; 399.706 406.457; 452.261 453.179; 486.338 487.102; 498.869 499.579; 507.968 508.925; 546.266 547.285; 558.825 560.128]};
     BadSegments{10} = {...
         [235.104 236.184; 324.272 325.028; 330.799 331.401; 541.062 541.826; 564.747 565.072], ...
         [100.581 101.229; 285.644 285.644; 285.644 285.644; 297.979 298.033; 298.999 298.999; 300.949 301.026; 301.999 301.999; 303.999 303.999; 304.999 304.999; 306.999 306.999; 307.999 307.999; 310.999 310.999; 311.999 311.999; 313.971 314.025; 314.999 314.999; 316.995 317.035; 317.999 317.999; 319.999 319.999; 320.923 321.046; 326.971 327.048; 327.999 327.999; 329.999 329.999; 330.999 330.999; 332.974 333.028; 333.999 333.999; 335.971 336.025; 336.999 336.999; 338.973 339.035; 339.999 339.999; 341.951 342.044; 343.999 343.999; 344.964 345.033; 346.999 346.999; 347.999 347.999; 452.015 453.411; 458.735 459.776; 467.974 468.089], ...
         [41.266 41.675; 53.905 55.240; 159.159 159.345; 220.255 220.394], ...
         [66.751 67.190; 201.674 201.812; 294.119 295.168; 303.188 303.528; 316.992 317.037; 317.999 317.999; 317.999 317.999; 319.999 319.999; 320.980 321.042; 322.999 322.999; 325.999 325.999; 326.955 327.048; 328.999 328.999; 333.999 333.999; 334.999 334.999; 342.999 342.999; 343.999 343.999; 351.915 352.038; 352.999 352.999; 407.979 408.056; 409.999 409.999; 411.999 411.999; 415.405 416.154; 435.999 435.999; 436.962 437.046; 468.885 469.100; 469.999 469.999; 470.999 470.999; 516.210 516.357], ...
         [105.369 106.172; 141.468 142.178; 199.008 199.818; 201.269 201.716; 352.851 353.083; 449.865 452.041; 459.508 459.802], ...
         [229.948 230.033; 230.999 230.999; 230.999 230.999; 259.588 259.990; 343.970 345.798; 361.999 361.999; 362.970 363.046; 364.999 364.999; 366.999 366.999; 367.962 368.031; 372.521 374.179; 405.999 405.999; 409.999 409.999; 411.999 411.999; 412.999 412.999; 434.999 434.999; 435.970 436.031; 516.999 516.999; 519.999 519.999]};
     BadSegments{11} = {...
         [179.045 180.225; 323.999 323.999; 323.999 323.999; 324.999 324.999; 327.978 328.025; 328.966 329.028; 330.999 330.999; 365.999 365.999; 367.999 367.999; 370.999 370.999; 371.999 371.999; 373.999 373.999; 375.965 376.042; 377.999 377.999], ...
         [52.107 53.156; 521.077 521.155], ...
         [65.331 66.156], ...
         [81.579 82.143; 108.565 108.565; 108.566 108.566; 112.176 112.176; 122.377 122.377; 122.565 122.565; 213.882 213.882; 215.305 215.305; 224.851 224.851; 224.919 224.919; 255.815 255.815; 257.893 257.893; 359.952 359.952; 361.630 361.754; 370.285 370.285; 376.853 376.853; 511.600 511.600; 513.631 513.631; 513.988 513.988; 518.215 518.215], ...
         [63.205 64.154; 170.951 171.036; 176.841 176.841; 176.842 176.842; 177.282 177.282; 223.886 223.971; 259.963 259.963; 261.547 261.547; 341.185 341.302; 368.999 368.999; 370.999 370.999; 374.999 374.999; 382.971 383.033; 383.999 383.999; 386.999 386.999; 388.958 389.035; 390.999 390.999; 391.955 392.031; 394.955 395.040; 396.999 396.999; 398.999 398.999; 400.999 400.999; 402.999 402.999; 404.963 405.033; 406.999 406.999; 429.829 429.829; 430.346 430.346; 465.184 465.184; 470.290 470.290], ...
         [50.195 51.145; 158.850 159.012]};
     BadSegments{12} = {...
         [193.435 194.175; 500.000 501.000; 08.988 509.868; 528.999 528.999; 528.999 528.999; 529.999 529.999; 531.999 531.999; 547.999 547.999], ...
         [133.484 134.185; 134.911 135.065; 553.999 553.999; 553.999 553.999; 554.999 554.999; 557.999 557.999; 558.999 558.999; 564.999 564.999; 565.977 566.046; 568.999 568.999; 569.951 570.028; 571.999 571.999; 579.959 580.028; 580.977 581.055; 583.971 584.033; 583.999 583.999; 585.999 585.999; 586.999 586.999; 588.980 589.026; 590.999 590.999; 591.958 592.044; 595.999 595.999], ...
         [46.028 47.216; 129.557 130.236; 200.999 200.999; 200.999 200.999; 201.975 202.075; 203.999 203.999; 204.967 205.053; 213.627 214.515; 218.958 219.044; 222.959 228.252; 309.999 309.999; 311.897 317.537; 311.999 311.999; 351.968 353.234; 399.999 399.999; 446.508 451.454; 487.999 487.999; 512.278 512.926], ...
         [82.535 84.062; 102.247 102.355; 134.999 134.999; 134.999 134.999; 136.999 136.999; 138.973 139.042; 147.999 147.999; 148.999 148.999; 149.999 149.999; 193.999 193.999; 194.999 194.999; 224.984 226.095; 241.954 242.077; 243.948 244.087; 280.999 280.999; 281.999 281.999; 305.000 305.451; 310.945 311.068; 328.999 328.999; 352.999 352.999; 353.999 353.999; 397.025 397.411; 415.944 418.096; 415.999 415.999; 470.999 470.999; 477.971 478.041; 483.997 484.737; 492.958 493.089; 493.999 493.999; 494.999 494.999; 522.999 522.999; 523.999 523.999; 524.955 525.039], ...
         [175.183 176.155; 246.112 246.991; 412.005 413.340; 483.915 484.061; 485.959 486.028; 497.360 498.880; 594.457 594.944; 596.641 598.006; 615.277 618.710], ...
         [66.900 72.232; 75.510 78.688; 543.495 545.999]};
     BadSegments{13} = {...
         [307.246 308.179; 627.881 628.089; 629.941 630.049], ...
         [171.506 172.301; 172.999 172.999; 430.999 430.999; 432.999 432.999; 434.999 434.999; 448.999 448.999; 479.999 479.999; 489.999 489.999; 490.999 490.999; 491.999 491.999], ...
         [92.059 93.209], ...
         [52.791 53.231; 54.999 54.999; 65.985 68.138; 91.240 91.386; 92.137 92.207; 105.346 105.493; 121.120 121.398; 146.546 146.955; 194.025 197.652; 242.985 243.571; 247.101 247.556; 269.889 270.083; 270.946 271.070; 274.408 275.866; 294.565 295.267; 319.370 319.879; 365.999 365.999; 375.913 376.044; 376.869 377.055; 377.999 377.999; 390.130 393.225; 403.965 404.035; 404.953 405.075; 419.265 419.565; 427.015 427.154; 427.879 428.118; 427.999 427.999; 428.913 429.067; 480.945 482.095; 484.285 484.571; 516.929 517.099; 517.963 518.079], ...
         [115.921 116.060; 117.999 117.999; 120.999 120.999; 126.896 127.058; 130.955 131.039; 131.957 132.043; 132.961 133.015; 134.999 134.999; 135.950 136.043; 139.999 139.999; 143.999 143.999; 144.850 145.066; 145.999 145.999; 146.999 146.999; 160.929 161.045; 161.999 161.999; 189.086 189.310; 193.939 194.063; 195.975 196.045; 199.961 200.045; 201.952 202.059; 202.939 205.084; 210.999 210.999; 211.738 213.629; 217.954 218.061; 218.987 219.049; 219.975 220.028; 258.953 259.037; 259.952 260.045; 260.963 261.032; 261.958 262.043; 262.954 263.054; 262.999 262.999; 288.999 288.999; 295.477 295.655; 298.976 299.054; 299.979 300.018; 300.975 301.005; 336.952 337.044; 337.977 338.032; 338.988 339.027; 338.999 338.999; 339.991 340.045; 340.975 341.029; 341.978 342.016; 343.976 344.054; 361.961 362.045; 365.942 366.050; 370.999 370.999; 384.961 385.045; 409.977 410.047; 415.984 416.030; 426.944 427.059; 430.999 430.999; 436.999 436.999; 448.999 448.999; 449.500 452.131; 481.962 482.055; 482.973 483.027; 484.975 485.060; 503.999 503.999; 504.999 504.999; 507.967 508.028; 508.955 509.032; 510.999 510.999; 512.999 512.999; 518.865 520.292; 522.958 523.043; 523.999 523.999; 524.948 525.034; 525.959 526.044; 526.946 527.055; 527.942 528.027; 528.961 529.053; 529.945 530.022; 535.999 535.999; 536.957 537.058; 538.967 539.036; 539.963 540.040; 541.976 542.045; 543.999 543.999; 551.956 552.034; 552.990 553.028; 553.977 554.024; 555.968 556.022; 556.204 557.454; 557.963 558.032; 559.969 560.023; 563.973 564.035; 563.999 563.999], ...
         [87.308 87.639; 107.039 108.189; 427.943 428.035; 437.965 438.027; 439.973 440.026; 491.275 492.571]};
     BadSegments{14} = {...
         [365.982 367.194; 368.999 368.999; 368.999 368.999; 369.999 369.999; 371.999 371.999; 373.999 373.999; 375.942 376.050; 376.999 376.999; 378.999 378.999; 380.999 380.999; 382.999 382.999; 383.999 383.999; 386.999 386.999; 387.999 387.999; 418.999 418.999; 444.999 444.999; 690.178 693.218], ...
         [101.796 102.183; 218.999 218.999; 219.999 219.999; 221.999 221.999; 222.999 222.999; 227.999 227.999; 229.999 229.999; 230.999 230.999; 232.999 232.999; 233.999 233.999; 235.999 235.999; 237.999 237.999; 238.985 239.031; 240.999 240.999; 242.999 242.999; 243.999 243.999; 246.999 246.999; 247.999 247.999; 253.999 253.999; 258.999 258.999; 259.999 259.999; 260.227 261.123; 264.999 264.999; 265.999 265.999; 268.999 268.999; 269.999 269.999; 274.999 274.999; 276.999 276.999; 277.999 277.999; 280.999 280.999; 282.966 283.021; 284.999 284.999; 285.999 285.999; 360.999 360.999; 362.973 363.026; 364.999 364.999; 365.999 365.999; 367.999 367.999; 368.985 369.024; 370.999 370.999; 371.999 371.999; 378.999 378.999; 379.999 379.999; 381.999 381.999; 383.999 383.999; 385.999 385.999; 387.999 387.999; 388.999 388.999; 398.999 398.999; 415.999 415.999; 417.999 417.999; 418.999 418.999; 421.999 421.999; 422.999 422.999], ...
         [89.285 90.187], ...
         [85.991 87.179; 88.999 88.999; 89.999 89.999; 92.966 93.027; 93.999 93.999; 107.999 107.999; 110.999 110.999; 111.971 112.025; 118.999 118.999; 121.999 121.999; 122.999 122.999; 166.999 166.999; 168.957 169.035; 173.999 173.999; 174.999 174.999; 175.957 176.035; 198.999 198.999; 199.999 199.999; 202.999 202.999; 203.999 203.999; 205.999 205.999; 207.999 207.999; 208.999 208.999; 210.974 211.044; 211.999 211.999; 212.999 212.999; 232.984 233.023; 234.999 234.999; 235.999 235.999; 236.999 236.999; 239.958 240.044; 240.999 240.999], ...
         [71.132 72.227], ...
         [98.134 98.827; 110.497 110.814; 114.000 117.272; 279.999 279.999; 279.999 279.999; 280.999 280.999; 282.999 282.999; 283.999 283.999; 284.999 284.999; 286.999 286.999; 287.999 287.999; 288.999 288.999; 289.999 289.999; 291.966 292.044; 292.999 292.999; 346.999 346.999; 437.608 437.646; 518.988 519.104; 522.969 523.015; 534.999 534.999; 536.967 537.029; 562.725 563.426]};
     BadSegments{15} = {...
         [66.755 67.172; 257.988 258.042; 258.714 259.385; 260.958 261.044; 265.931 266.062; 267.968 268.022; 268.978 269.033; 270.992 271.031; 272.964 273.017; 273.968 274.030; 275.975 276.059; 276.970 277.046; 288.976 289.038; 289.984 290.015; 295.994 296.025; 296.982 297.013; 300.966 301.152; 364.730 364.800], ...
         [46.531 46.555; 65.156 66.148; 87.535 87.642; 183.896 183.934; 294.342 294.375; 330.004 330.528; 333.785 333.815; 345.177 345.203; 399.544 399.573; 480.964 480.980], ...
         [53.557 54.175; 148.851 149.136], ...
         [43.681 44.214; 361.255 361.587; 409.462 411.021], ...
         [90.868 91.130; 410.935 411.044], ...
         [64.786 65.141; 132.385 132.786; 177.735 178.252; 278.747 278.902; 313.959 314.028; 314.970 315.031]};
     BadSegments{16} = {...
         [38.958 44.335; 305.280 312.826; 324.778 326.067; 393.183 398.367; 403.410 422.854], ...
         [49.726 50.251; 58.675 65.465; 263.918 264.088; 321.795 322.195; 522.881 523.081], ...
         [61.245 62.325; 65.628 65.905; 315.686 317.275; 335.184 336.449; 387.870 388.117; 389.892 390.062], ...
         [56.020 57.255; 60.711 61.096; 64.353 64.924; 73.202 73.757; 76.134 76.844; 95.625 96.859; 171.436 172.625; 178.359 178.715; 212.828 213.415; 352.425 352.617; 367.755 369.097; 488.400 489.426], ...
         [71.665 72.513; 76.347 77.181; 127.134 128.862], ...
         [46.190 47.209; 53.461 55.752; 194.317 194.510; 213.101 213.240; 216.356 216.442; 225.962 226.286; 245.163 245.464; 252.041 253.422; 254.648 254.856; 292.966 294.340; 340.928 341.067; 346.298 346.452]};
    BadSeg = BadSegments{iSubj}{iRun}';
end

author: brainstorm-tools


text:
%Speed comparison between isStopWord functions implemented with:
    %1)cellArray
    %2)containers.map

n = 1000;

text = fileread('C:\Users\FLOW\Documents\Internships_Research\faridani-MatlabNLP-2e75adc\isStopWord\test.txt');
text = sanitizeComment(text);   %remove all but letters
text = lower(text);
text = regexp(text,' ','split'); %text now in cell array
%text = text(1:n);   %determine size of test text based on given n

testMap = containers.Map;
for i = 1:size(text,2)
    word = porterStemmer(cell2mat(text(i)));
    if isKey(testMap,word)
    else
        testMap(word) = word;
    end
end

testCellArray = keys(testMap);

tic;

mapOut = isStopWordMap('aura',testMap);

mapTime = toc;

tic;

caOut = isStopWord('aura',testCellArray);

caTime = toc;

fprintf('map implementation of isStopWord returns %d in %f\n',mapOut,mapTime)
fprintf('cell array implementation of isStopWord returns %d in %f\n',caOut,caTime)

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function rgbVector = str2rgb(colorString)
%STR2RGB   Converts a string representation of a color to an RGB triple.
%   X = STR2RGB(STR) converts the string STR into a three-element row
%   vector X (an RGB triple). STR can be one of three string
%   representations of colors that MATLAB uses (see ColorSpec help): a full
%   color name ('yellow'), a single character ('y'), or a string of numbers
%   within the range of 0 and 1 ('[1 1 0]' or '1 1 0').
%
%   If the string STR does not represent a valid string representation of a
%   color, STR2RGB(STR) returns NaN.
%
%   NOTE: STR2RGB does not use eval.
%
%   Examples:
%      str2rgb('yellow')        returns [1 1 0]
%      str2rgb('y')             returns [1 1 0]
%      str2rgb('[1 1 0]')       returns [1 1 0]
%      str2rgb('1 1 0')         returns [1 1 0]
%      str2rgb('[1; 1; 0]')     returns [1 1 0]
%      str2rgb('[0 0.5 0.91]')  returns [0 0.5000 0.9100]
%      str2rgb('purple')        returns NaN
%      str2rgb('[1 2]')         returns NaN

% Author: Ken Eaton
% Last modified: 4/2/08
%--------------------------------------------------------------------------

  if (nargin == 0),
    error('str2rgb:notEnoughInputs','Not enough input arguments.');
  end
  if (~ischar(colorString)),
    error('str2rgb:badArgumentType',...
          'Input argument should be of type char.');
  end
  expression = {'^\s*yellow\s*$','^\s*magenta\s*$','^\s*cyan\s*$',...
                '^\s*red\s*$','^\s*green\s*$','^\s*blue\s*$',...
                '^\s*white\s*$','^\s*black\s*$','^\s*y\s*$','^\s*m\s*$',...
                '^\s*c\s*$','^\s*r\s*$','^\s*g\s*$','^\s*b\s*$',...
                '^\s*w\s*$','^\s*k\s*$','[\[\]\;\,]'};
  replace = {'[1 1 0]','[1 0 1]','[0 1 1]','[1 0 0]','[0 1 0]',...
             '[0 0 1]','[1 1 1]','[0 0 0]','[1 1 0]','[1 0 1]',...
             '[0 1 1]','[1 0 0]','[0 1 0]','[0 0 1]','[1 1 1]',...
             '[0 0 0]',' '};
  rgbVector = sscanf(regexprep(colorString,expression,replace),'%f').';
  if ((numel(rgbVector) ~= 3) || any((rgbVector < 0) | (rgbVector > 1))),
    rgbVector = nan;
  end

end

author: Sable

repo: mcbench-benchmarks


text:
function gr_test07 ( node_num, node_coordinates, edge_num, edge_nodes );

%*****************************************************************************80
%
%% GR_TEST07 tests GR_INCIDENCE_MATRIX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GR_TEST07:\n' );
  fprintf ( 1, '  GR_INCIDENCE_MATRIX computes the incidence matrix.\n' );

  incidence_matrix = gr_incidence_matrix ( node_num, node_coordinates, ...
    edge_num, edge_nodes );

  i4mat_print ( edge_num, node_num, incidence_matrix, ...
    '  The incidence matrix:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Code to perform Raycasting on a Supplied Map to return simulated Lidar
% Range values.
% % Author- Shikhar Shrestha, IIT Bhubaneswar

% Input parameters are x-position, y-position, Binary image of map, number
% of rays to cast, Lidar range in cm.
% % Returned Range values are obtained on scanning clockwise.
function [range]=castrays(xc,yc,map,n,lidarrange)
% % Conversion factor from pixels to centimeters.
pixeltocm=1;
if map(ceil(yc),ceil(xc))==1
    range=zeros(n,1);
    return;
end
figure(2);
hold off;
imshow(map);
hold on;
plot(xc,yc,'b*');

if nargin==3
    n=20;
    lidarrange=500;
end

if nargin==4
lidarrange=500;
end

range=zeros(n,1);
thetastep=360/n;

for i=1:n
theta=thetastep*i;
r=linspace(0,lidarrange,1000);
x=xc+(r*cosd(theta));
y=yc+(r*sind(theta));
% Removing points out of map
temp=[];
for k=1:numel(x)
    if x(k)>size(map,2) || y(k)>size(map,1) || x(k)<=0 || y(k)<=0
        temp=[temp;k];
    end
end
x(temp)=[];
y(temp)=[];
figure(2);
plot(x,y,'r');
% Computing Intersections
xint=round(x);
yint=round(y);
% Correcting zero map indexing
for l=1:numel(xint)
    if xint(l)==0
        xint(l)=1;
    end
    if yint(l)==0
        yint(l)=1;
    end
end

b=[];
for j=1:numel(xint)
b=[b;map(yint(j),xint(j))];
end
ind=find(b==1);

if ~isempty(ind)
xb=x(ind(1));
yb=y(ind(1));
figure(2);
plot(xb,yb,'g*');
dist=sqrt((xc-xb).^2 + (yc-yb).^2);
range(i)=dist;
end
pause(0.000001);
end
% Converting to mm from pixels.
range=range*10*pixeltocm;
% % % % % % % range=flipud(range);


text:
% needs ops.RegFile, ops.xrange, ops.yrange, ops.NavgFramesSVD
% Ly, Lx are the size of each frame. nimgbatch: the number of frames loaded per batch
% nt0 is the number of timepoints to bin over. If a sixth argument is
% present, it does not subtract the mean of each batch.
function mov = loadAndBin(ops, Ly, Lx, nimgbatch, nt0, clustModel)

ix = 0;
fid = fopen(ops.RegFile, 'r');
mov = zeros(numel(ops.yrange), numel(ops.xrange), ops.NavgFramesSVD, 'single');
ij = 0;
while 1
    % load frames
    data = fread(fid,  Ly*Lx*nimgbatch, '*int16');
    if isempty(data)
        break;
    end
    data = single(data);
    data = reshape(data, Ly, Lx, []);

    % ignore bad frames
    badi = ops.badframes(ix + [1:size(data,3)]);
%     data(:,:, badi) = [];

    % subtract off the mean of this batch
    if nargin<=5
        data = bsxfun(@minus, data, mean(data,3));
    end
    %     data = bsxfun(@minus, data, ops.mimg1);

    nSlices = nt0*floor(size(data,3)/nt0);
    if nSlices~=size(data,3)
        data = data(:,:, 1:nSlices);
    end

    % bin data
    data = reshape(data, Ly, Lx, nt0, []);
    davg = squeeze(mean(data,3));

    mov(:,:,ix + (1:size(davg,3))) = davg(ops.yrange, ops.xrange, :);

    ix = ix + size(davg,3);
    ij = ij + 1;
end
fclose(fid);

mov = mov(:, :, 1:ix);

author: cortex-lab


text:
classdef nnloss < nntest
  properties (TestParameter)
    loss = {...
      'classerror', 'log', 'softmaxlog', 'mhinge', 'mshinge', ...
      'binaryerror', 'binarylog', 'logistic', 'hinge'}
    weighed = {false, true}
  end

  properties
    x
  end

  methods
    function [x,c,dzdy,instanceWeights] = getx(test,loss)
      numClasses = 3 ;
      numAttributes = 5 ;
      numImages = 3 ;
      w = 5 ;
      h = 4 ;
      switch loss
        case {'log', 'softmaxlog', 'mhinge', 'mshinge', 'classerror'}
          % multiclass
          instanceWeights = test.rand(h,w, 'single') / test.range / (h*w) ;
          c = randi(numClasses, h,w,1,numImages) ;
          c = test.toDevice(c) ;
        otherwise
          % binary
          instanceWeights = test.rand(h,w, numAttributes, 'single') / test.range / (h*w*numAttributes) ;
          c = sign(test.randn(h,w,numAttributes, numImages)) ;
      end
      c = single(c) ;
      switch loss
        case {'log'}
          x = test.rand(h,w, numClasses, numImages, 'single') / test.range * .60 + .20 ;
          x = bsxfun(@rdivide, x, sum(x,3)) ;
        case {'binarylog'}
          x = test.rand(h,w, numAttributes, numImages, 'single') / test.range * .60 + .20 ;
        case {'softmaxlog'}
          x = test.randn(h,w, numClasses, numImages, 'single') / test.range ;
        case {'mhinge', 'mshinge', 'classerror'}
          x = test.randn(h,w, numClasses, numImages, 'single') / test.range ;
        case {'hinge', 'logistic', 'binaryerror'}
          x = test.randn(h,w, numAttributes, numImages, 'single') / test.range ;
      end
      dzdy = test.randn(1,1) / test.range ;
    end
  end

  methods (Test)
    function nullcategories(test, loss, weighed)
      [x,c,dzdy,instanceWeights] = test.getx(loss) ;
      % make a number of categories null
      c(:) = c(:) .* (test.randn(numel(c),1) > 0) ;
      opts = {'loss',loss} ;
      if weighed, opts = {opts{:}, 'instanceWeights', instanceWeights} ; end
      y = vl_nnloss(x,c,[],opts{:}) ;
      dzdx = vl_nnloss(x,c,dzdy,opts{:}) ;
      test.der(@(x) vl_nnloss(x,c,[],opts{:}), x, dzdy, dzdx, 0.001, -5e-1) ;
    end

    function convolutional(test, loss, weighed)
      [x,c,dzdy,instanceWeights] = test.getx(loss) ;
      opts = {'loss',loss} ;
      if weighed, opts = {opts{:}, 'instanceWeights', instanceWeights} ; end
      y = vl_nnloss(x,c,[],opts{:}) ;
      dzdx = vl_nnloss(x,c,dzdy,opts{:}) ;
      test.der(@(x) vl_nnloss(x,c,[],opts{:}), x, dzdy, dzdx, 0.001, -5e-1) ;
    end

  end
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function a = dirichlet_fit_s_simple(data,a)
% DIRICHLET_FIT_S_SIMPLE   Initial guess for Dirichlet precision.
%
% DIRICHLET_FIT_S_SIMPLE(data,a) returns an initial guess for the Dirichlet
% parameter vector A, by scaling the input A.

bar_p = mean(log(data));
m = a/sum(a);
bar_p = sum(m.*bar_p);
s = 1/(sum(m.*log(m)) - bar_p);
s = s*(length(m)-1)/2;
a = s*m;

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
% test for hessian tensor

n = 256;
name = 'barb';
M = load_image(name);
M = rescale(crop(M,n));

options.order = 2;
T = compute_hessian(M, options);

[e1,e2,l1,l2] = perform_tensor_decomp(T);
l1 = abs(l1); l2 = abs(l2);
T = perform_tensor_recomp(e1,e2,l2,l1);

sigma = 5;
T = perform_blurring(T,sigma);

options.sub = 8;
clf;
plot_tensor_field(T, M, options);
figure;
clf;
plot_tensor_field_nb(T, M, options);

author: gpeyre

repo: matlab-toolboxes


text:
function [nodes, edges, ages, error] = removeUnconnected(nodes, edges,ages,error) %#codegen

NumOfNodes = size(nodes,2);

i = 1;
while NumOfNodes >= i
    if any(edges(i,:)) == 0

        edges = [edges(1:i-1,:); edges(i+1:NumOfNodes,:);];
        edges = [edges(:,1:i-1)  edges(:,i+1:NumOfNodes);];

        ages = [ages(1:i-1,:); ages(i+1:NumOfNodes,:);];
        ages = [ages(:,1:i-1)  ages(:,i+1:NumOfNodes);];

        nodes = [nodes(:,1:i-1) nodes(:,i+1:NumOfNodes);];
        error = [error(1,1:i-1) error(1,i+1:NumOfNodes);];

        NumOfNodes = NumOfNodes - 1;

        i = i -1;
    end
    i = i+1;
end

% MEX-code generation

% mexcfg = coder.config('mex');
% mexcfg.DynamicMemoryAllocation = 'AllVariableSizeArrays';
% codegen -config mexcfg removeUnconnected -args {coder.typeof(nodes,[Inf Inf]),coder.typeof(edges,[Inf Inf]),coder.typeof(ages,[Inf Inf]),coder.typeof(error,[1 Inf])}


text:
function obj = checkToolboxes(obj)
	% Biafra Ahanonu
	% Started: 2021.03.25 [22:11:25] (branched from ciatah.m)

	license('inuse')
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function linplus_test03 ( )

%*****************************************************************************80
%
%% TEST03 tests R83_CR_FA, R83_CR_SL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03\n' );
  fprintf ( 1, '  For a real tridiagonal matrix,\n' );
  fprintf ( 1, '  using CYCLIC REDUCTION,\n' );
  fprintf ( 1, '  R83_CR_FA factors;\n' );
  fprintf ( 1, '  R83_CR_SL solves.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
  fprintf ( 1, '  The matrix is NOT symmetric.%d\n' );
  fprintf ( 1, '\n' );
%
%  Set the matrix values.
%
  a(1,1) = 0.0E+00;
  for j = 2 : n
    a(1,j) = j;
  end

  for j = 1 : n
    a(2,j) = 4.0E+00 * j;
  end

  for j = 1 : n - 1
    a(3,j) = j;
  end
  a(3,n) = 0.0E+00;

  r83_print ( n, a, '  The matrix:' );
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
  b = r83_mxv ( n, a, x );
  x(1:n) = 0.0E+00;
%
%  Factor the matrix.
%
  a_cr = r83_cr_fa ( n, a );
%
%  Solve the linear system.
%
  x = r83_cr_sl ( n, a_cr, b );

  r8vec_print_some ( n, x, 10, '  Solution:' );

  return
end


text:
function varargout=abs(varargin)
%ABS (overloaded)

switch class(varargin{1})
    case 'double'
        error('Overloaded SDPVAR/ABS CALLED WITH DOUBLE. Report error')

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.
        if isreal(varargin{1})
            varargout{1} = yalmip('define',mfilename,varargin{1});
            B = getbase(varargin{1});
            if all(all(fix(B)==B))
                CurrInt = yalmip('quantvariables');
                if all(ismember(getvariables(varargin{1}),CurrInt))
                    yalmip('setintvariables',[CurrInt getvariables(varargout{1})]);
                end
            end
        else
            % For complex args, abs(X) is defined [norm(X(i,j),2)] in MATLAB
            y = [];
            x = varargin{1};
            for i = 1:size(x,1)
                temp = [];
                for j = 1:size(x,2)
                    temp = [temp norm(extsubsref(x,i,j))];
                end
                y = [y;temp];
            end
            varargout{1} = y;
        end

    case 'char' % YALMIP send 'graph' when it wants the epigraph or hypograph
        switch varargin{1}
            case 'graph'
                % Description using epigraphs
                t = varargin{2};
                X = varargin{3};
                varargout{1} = [1 -1;-1 -1]*[X;t] <= [0;0];
                varargout{2} = struct('convexity','convex','monotonicity','none','definiteness','positive','model','graph');
                varargout{3} = X;

            case {'exact','integer','callback'}
                % Exact description using binary variables
                t = varargin{2};
                X = varargin{3};
                d = varargin{4};
                [M,m]=derivebounds(X);
                if m>=0
                    F = (t == X);
                elseif M<=0
                    F = (t == -X);
                else
                    maxABSX = max([abs(m) abs(M)],[],2);
                    F = [[0 1 0;
                     0 -1 0;
                     1 0 -M;
                     1 1 -2*maxABSX;
                     -1 -1 0;
                     -1 0 -m;
                     -1 1 2*maxABSX;1 -1 0]*[X;t;d] <= [maxABSX;0;0;0;0;-m;2*maxABSX;0]];
                end

                varargout{1} = F;
                varargout{2} = struct('convexity','convex','monotonicity','none','definiteness','positive','model','integer');
                varargout{3} = X;

            otherwise
                error('SDPVAR/ABS called with CHAR argument?');
        end
    otherwise
        error('Strange type on first argument in SDPVAR/ABS');
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function xMax = SmoothMax(x,alpha)
%
% xMax = SmoothMax(x, alpha)
%
% This function is a smooth version of the max(x) function, commonly known
% as the 'soft-max' function. In the limit as the smoothing parameter goes
% to zero, this function will return max(x). In the limit as the smoothing
% parameter goes to infinity, this will return sum(x).
%
%INPUTS:
% x = a vector or matrix of inputs
% alpha = a smoothing parameter. Asymtotes:  0->max(x), inf->sum(x)
%
%OUTPUTS:
% xMax = the smooth maximum of x
%
%
% Written by Matthew Kelly
% October 2013
% Cornell University
%

if nargin==1
    alpha = 1;
end

xMax = alpha*log(sum(exp(x/alpha)));

if isinf(xMax)
   %Then the scaling was too sharp
    xMax = max(x);
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
%   TTeMPS Toolbox.
%   Michael Steinlechner, 2013-2016
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt

function [X, residual, cost, times] = amen_fast( L, F, X, opts )

t_start = tic();
% set default opts
if ~exist( 'opts', 'var');        opts = struct();     end
if ~isfield( opts, 'nSweeps');    opts.nSweeps = 4;    end
if ~isfield( opts, 'maxrank');    opts.maxrank = 20;   end
if ~isfield( opts, 'maxrankRes'); opts.maxrankRes = 4; end
if ~isfield( opts, 'tolRes');     opts.tolRes = 1e-13;  end
if ~isfield( opts, 'tol');        opts.tol = 1e-13;     end
if ~isfield( opts, 'solver');     opts.solver = 'direct';   end
if ~isfield( opts, 'prec');       opts.prec = true;   end

d = X.order;
n = X.size;

normF = norm(F);
cost = cost_function( L, X, F );
residual = norm( apply(L, X) - F ) / normF;
times = toc(t_start);

for sweep = 1:opts.nSweeps
    X = orthogonalize(X, 1);
    for mu = 1:d-1
        disp( ['Current core: ', num2str(mu)] )

        % STEP 1: Solve mu-th core opimization
        F_mu = contract( X, F, mu );
        sz = [X.rank(mu), X.size(mu), X.rank(mu+1)];

        if strcmpi( opts.solver, 'direct' )
            % if system very small
            L_mu = contract( L, X, mu );
            U_mu = L_mu \ F_mu(:);
            X.U{mu} = reshape( U_mu, sz );
        elseif strcmpi( opts.solver, 'pcg' )
            [left, right] = Afun_prepare( L, X, mu );
            [B2, V, E] =  prepare_precond( L.L0, X, mu );

            U_mu = pcg( @(y) Afun( L, y, mu, sz, left, right), ...
                     F_mu(:), ...
                     1e-10, 1000, ...
                     @(y) apply_precond( B2, V, E, y, sz ), [],...
                     X.U{mu}(:) );
            X.U{mu} = reshape( U_mu, sz );
        else
            error( 'Unknown opts.solver type. Use either ''direct'' (default) or ''pcg''.' )
        end

        % STEP 2: Calculate current residual and cost function
        res =  F - apply(L, X);
        residual = [residual; norm( res ) / normF];
        cost = [cost; cost_function( L, X, F )];
        disp(['Rel. residual: ' num2str(residual(end)) ', Current rank: ' num2str(X.rank) ]);

        % STEP 3: Augment mu-th and (mu+1)-th core with (truncated) residual
        R = contract( X, res, [mu, mu+1] );
        R_combined = unfold(R{1},'left') * unfold(R{2},'right');
        if opts.prec
            R_combined = precond_residual( L.L0, X, R_combined, mu );
        end
        [uu,ss,~] = svd( R_combined, 'econ');
        s = find( diag(ss) > opts.tolRes*norm(diag(ss)), 1, 'last' );
		if opts.maxrankRes ~= 0
			s = min( s, opts.maxrankRes );
		end
        R{1} = reshape( uu(:,1:s)*ss(1:s,1:s), [X.rank(mu), n(mu), s]);
        %R{2} = reshape( vv(:,1:s)', [s, n(mu+1), X.rank(mu+2)]);

        left = cat(3, X.U{mu}, R{1});
        %right = cat(1, X.U{mu+1}, R{2});

        % STEP 4: Move orthogonality to (mu+1)-th core while performing rank truncation
        [U,S,~] = svd( unfold(left,'left'), 'econ' );
        t = find( diag(S) > opts.tol*norm(diag(S)), 1, 'last' );
        t = min( t, opts.maxrank );
        X.U{mu} = reshape( U(:,1:t), [X.rank(mu), n(mu), t] );
        X.U{mu+1} = rand( t, n(mu+1), X.rank(mu+2));

        times = [times; toc(t_start)];
    end
    for mu = d:-1:2
        disp( ['Current core: ', num2str(mu)] )

        % STEP 1: Solve mu-th core opimization
        F_mu = contract( X, F, mu );
        sz = [X.rank(mu), X.size(mu), X.rank(mu+1)];

        if strcmpi( opts.solver, 'direct' )
            L_mu = contract( L, X, mu );
            U_mu = L_mu \ F_mu(:);
            X.U{mu} = reshape( U_mu, size(X.U{mu}) );
        elseif strcmpi( opts.solver, 'pcg' )
            [left, right] = Afun_prepare( L, X, mu );
            [B2, V, E] =  prepare_precond( L.L0, X, mu );

            U_mu = pcg( @(y) Afun( L, y, mu, sz, left, right), ...
                     F_mu(:), ...
                     1e-10, 1000, ...
                     @(y) apply_precond( B2, V, E, y, sz ), [],...
                     X.U{mu}(:) );
            X.U{mu} = reshape( U_mu, sz );
        else
            error( 'Unknown opts.solver type. Use either ''direct'' (default) or ''diag''.' )
        end

        % STEP 2: Calculate current residual and cost function
        res =  F - apply(L, X);
        residual = [residual; norm( res ) / normF];
        disp(['Rel. residual: ' num2str(residual(end)) ', Current rank: ' num2str(X.rank) ]);
        cost = [cost; cost_function( L, X, F )];

        % STEP 3: Augment mu-th and (mu+1)-th core with (truncated) residual
        R = contract( X, res, [mu-1, mu] );
        R_combined = unfold(R{1},'left') * unfold(R{2},'right');
        if opts.prec
            R_combined = precond_residual( L.L0, X, R_combined, mu-1 );
        end
        [~,ss,vv] = svd( R_combined, 'econ');
        s = find( diag(ss) > opts.tolRes*norm(diag(ss)), 1, 'last' );
		if opts.maxrankRes ~= 0
			s = min( s, opts.maxrankRes );
		end
        R{2} = reshape( ss(1:s,1:s)*vv(:,1:s)', [s, n(mu), X.rank(mu+1)]);

        right = cat(1, X.U{mu}, R{2});

        % STEP 4: Move orthogonality to (mu+1)-th core while performing rank truncation
        [~,S,V] = svd( unfold(right,'right'), 'econ' );
        t = find( diag(S) > opts.tol*norm(diag(S)), 1, 'last' );
        t = min( t, opts.maxrank );
        X.U{mu} = reshape( V(:,1:t)', [t, n(mu), X.rank(mu+1)] );
        X.U{mu-1} = rand( X.rank(mu-1), n(mu-1), t);

        times = [times; toc(t_start)];
    end

end

end

function res = cost_function( L, X, F )
res = 0.5*innerprod( X, apply(L, X) ) - innerprod( X, F );
end

function res = euclid_grad( L, X, F )
res = apply(L, X) - F;
end

function res = precond_residual( L0, X, R_combined, idx )
    n = size(L0, 1);
    rl = X.rank(idx);
    rr = X.rank(idx+2);

    B1 = zeros( rl );
    % calculate B1 part:
    for i = 1:idx-1
        % apply L to the i'th core
        tmp = X;
        tmp.U{i} = tensorprod_ttemps( tmp.U{i}, L0, 2 );
        B1 = B1 + innerprod( X, tmp, 'LR', idx-1);
    end

    % calculate B2 part:
    B2 = kron( L0, speye(n) ) + kron( speye(n), L0 );

    B3 = zeros( rr );
    % calculate B3 part:
    for i = idx+2:X.order
        tmp = X;
        tmp.U{i} = tensorprod_ttemps( tmp.U{i}, L0, 2 );
        B3 = B3 + innerprod( X, tmp, 'RL', idx+2);
    end

    [V,E] = eig( kron( eye(rr), B1 ) + kron( B3, eye(rl) ) );
    E = diag(E);

    R_combined = reshape( R_combined, [rl, n*n, rr] );
    rhs = matricize( R_combined, 2 ) * V;
    Y = zeros(size(rhs));
    for i=1:length(E)
        Y(:,i) = (B2 + E(i)*speye(n*n)) \ rhs(:,i);
    end
    res = tensorize( Y*V', 2, [rl, n*n, rr] );
    res = reshape( res, [rl*n, n*rr] );
end

function [left, right] = Afun_prepare( A, x, idx )
    y = A.apply(x);
    if idx == 1
        right = innerprod( x, y, 'RL', idx+1 );
        left = [];
    elseif idx == x.order
        left = innerprod( x, y, 'LR', idx-1 );
        right = [];
    else
        left = innerprod( x, y, 'LR', idx-1 );
        right = innerprod( x, y, 'RL', idx+1 );
    end
end

function res = Afun( A, U, idx, sz, left, right )

    V = reshape( U, sz );
    V = A.apply( V, idx );

    if idx == 1
        tmp = tensorprod_ttemps( V, right, 3 );
    elseif idx == A.order
        tmp = tensorprod_ttemps( V, left, 1 );
    else
        tmp = tensorprod_ttemps( V, right, 3);
        tmp = tensorprod_ttemps( tmp, left, 1);
    end

    res = tmp(:);
end
function [B2, V, E] = prepare_precond( L0, X, idx )
    n = size(L0, 1);
    rl = X.rank(idx);
    rr = X.rank(idx+1);

    B1 = zeros( rl );
    % calculate B1 part:
    for i = 1:idx-1
        % apply L to the i'th core
        tmp = X;
        tmp.U{i} = tensorprod_ttemps( tmp.U{i}, L0, 2 );
        B1 = B1 + innerprod( X, tmp, 'LR', idx-1);
    end

    % calculate B2 part:
    B2 = L0;

    B3 = zeros( rr );
    % calculate B3 part:
    for i = idx+1:X.order
        tmp = X;
        tmp.U{i} = tensorprod_ttemps( tmp.U{i}, L0, 2 );
        B3 = B3 + innerprod( X, tmp, 'RL', idx+1);
    end

    [V,E] = eig( kron( eye(rr), B1 ) + kron( B3, eye(rl) ) );
    E = diag(E);
end
function res = apply_precond( B2, V, E, rhs, sz )
    n = size(B2, 1);
    rhs = reshape( rhs, sz );
    rhs = matricize( rhs, 2 ) * V;
    Y = zeros(size(rhs));
    for i=1:length(E)
        Y(:,i) = (B2 + E(i)*speye(n)) \ rhs(:,i);
    end
    res = tensorize( Y*V', 2, sz );
    res = res(:);
end

author: NicolasBoumal

repo: manopt


text:
%% Open-Loop Steering Simple
% This example simulates a simple vehicle with an open loop
% sinusoidal steering input.
%
% <<OpenLoopSteeringSimple.gif>>
%
%% Simulation models and parameters
% First, all classes of the package are imported with
%

clear ; close all ; clc

import VehicleDynamicsLateral.*

%% Tire and vehicle model
% Choosing tire and vehicle model. In this case, the parameters are defined
% by the user.
%

% Choosing tire
TireModel = TirePacejka();
% Choosing vehicle
System = VehicleSimpleNonlinear();

% Tire model and friction
System.tire = TireModel;
System.muy  = 1.0;

%% Simulation parameters
% Choosing simulation time span
%

T = 4;                      % Total simulation time [s]
resol = 50;                 % Resolution
TSPAN = 0:T/resol:T;        % Time span [s]

%% Open-loop steering input
% Steering angle
System.deltaf = 1*pi/180*sin(T^-1*2*pi*TSPAN);

%% Run simulation
% To define a simulation object (simulator) the arguments must be the
% vehicle object and the time span.

simulator = Simulator(System, TSPAN);

% Simulation
simulator.Simulate();

%% Results

% Retrieving states
XT      = simulator.XT;
YT      = simulator.YT;
PSI     = simulator.PSI;
VEL     = simulator.VEL;
ALPHAT  = simulator.ALPHAT;
dPSI    = simulator.dPSI;

figure(1)
plot(TSPAN,180/pi*System.deltaf)
xlabel('time [s]')
ylabel('Steering angle [deg]')

%%
% Frame and animation

g = Graphics(simulator);
g.TractorColor = 'r';

g.Frame();
g.Animation();
% g.Animation('html/OpenLoopSteeringSimple');       % Uncomment to save animation gif

%% See Also
%
% <../../../index.html Home> | <../../OpenLoopSteeringArticulated/html/OpenLoopSteeringArticulated.html Open-Loop Steering Articulated>
%

author: andresmendes

repo: Vehicle-Dynamics-Lateral


text:
function marginal = marginal_nodes(engine, nodes, t, add_ev)
% MARGINAL_NODES Compute the marginal on the specified query nodes (hmm)
% marginal = marginal_nodes(engine, nodes, t, add_ev)
%
% 'nodes' must be a single node.
% t is the time slice.

if nargin < 3, t = 1; end
if nargin < 4, add_ev = 0; end

assert(length(nodes)==1)
ss = engine.slice_size;

i = nodes(1);
bigT = engine.one_slice_marginal(:,t);
dom = i + (t-1)*ss;

ns = engine.eff_node_sizes(:);
bigdom = 1:ss;
marginal.T = marg_table(bigT, bigdom + (t-1)*ss, ns(bigdom), dom, engine.maximize);

marginal.domain = dom;
marginal.mu = [];
marginal.Sigma = [];

if add_ev
  marginal = add_ev_to_dmarginal(marginal, engine.evidence, engine.node_sizes);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function N1 = plus(N1, N2)
%PLUS   Plus for CHEBOP2.
%
% N = PLUS(N1, N2) is the same as N = N1 + N2.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Check domains are the same.
if ( all( N1.domain ~= N1.domain ) )
    error('CHEBFUN:CHEBOP2:plus:domain', ...
        'Domains of operators must be identical.');
end

% Plus the coefficient matrices together.
A = N1.coeffs;
B = N2.coeffs;

% PLUS the coefficient cellarrays or matrices together:
if ( iscell(A) && iscell(B) )
    C = cell(max(size(A, 1), size(B, 1)), max(size(A, 2), size(B, 2)));
    if size(B, 1) < size(C, 1)
        B{size(C, 1),1} = [];
    end
    if ( size(B, 2) < size(C, 2) )
        B{1,size(C, 2)} = [];
    end
    if ( size(A, 1) < size(C, 1) )
        A{size(C, 1),1} = [];
    end
    if ( size(A, 2) < size(C, 2) )
        A{1,size(C, 2)} = [];
    end
    for jj = 1:size(C, 1)
        for kk = 1:size(C, 2)
            if ( isempty(B{jj,kk}) )
                C{jj,kk} = A{jj,kk};
            elseif ( isempty(A{jj,kk}) )
                C{jj,kk} = B{jj,kk};
            else
                C{jj,kk} = A{jj,kk} + B{jj,kk};
            end
        end
    end

elseif ( iscell(A) )
    BB = cell(size(B, 2), size(B, 1));
    for jj = 1:size(B, 1)
        for kk = 1:size(B, 2)
            BB(kk,jj) = {B(jj,kk)};
        end
    end
    N2.coeffs = BB;
    N1 = plus(N1, N2);
    return

elseif ( iscell(B) )
    AA = cell(size(A, 2), size(A, 1));
    for jj = 1:size(A, 1)
        for kk = 1:size(A, 2)
            AA(kk,jj) = {A(jj,kk)};
        end
    end
    N1.coeffs = AA;
    N1 = plus(N1, N2);
    return

else
    C = zeros(max(size(A, 1), size(B, 1)), max(size(A, 2), size(B, 2)));
    if size(B, 1) < size(C, 1)
        B(size(B, 1)+1:size(C, 1),:) = 0;
    end
    if size(B, 2) < size(C, 2)
        B(:,size(B, 2)+1:size(C, 2)) = 0;
    end
    if size(A, 1) < size(C, 1)
        A(size(A, 1)+1:size(C, 1),:) = 0;
    end
    if size(A, 2) < size(C, 2)
        A(:,size(A, 2)+1:size(C, 2)) = 0;
    end
    C = A + B;
end

% Loop over to remove empty cells (replace with zeros).
for jj = 1:size(C, 1)
    for kk = 1:size(C, 2)
        if ( isempty(C{jj,kk}) )
            C{jj,kk} = 0;
        end
    end
end

% Do not add boundary conditions together. Ignore them.
if ( ~isempty(N1.lbc) || ~isempty(N1.rbc) ||...
        ~isempty(N1.ubc) || ~isempty(N1.dbc) )
    warning('CHEBFUN:CHEBOP2:plus:ignoredBCs', 'BCs were ignored.');
    N1.lbc =[];
    N1.rbc =[];
    N1.ubc =[];
    N1.dbc =[];
end
if ( ~isempty(N2.lbc) || ~isempty(N2.rbc) ||...
        ~isempty(N2.ubc) || ~isempty(N2.dbc) )
    warning('CHEBFUN:CHEBOP2:plus:ignoredBCs', 'BCs were ignored.');
    N2.lbc =[];
    N2.rbc =[];
    N2.ubc =[];
    N2.dbc =[];
end

% Assign to output:
N1.coeffs = C;
N1.xorder = max(N1.xorder, N2.xorder);
N1.yorder = max(N1.yorder, N2.yorder);
% Plus the ops together:
op1 = N1.op;
op2 = N2.op;
N1.op = @(u) op1(u) + op2(u);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% calculate error between an ellipse and a point
% expected input form : ellipse = [x0 y0 A B alpha], point = [x y]
function error = calcEllipseCost(ellipse, point)

% Initial setting
x = point(2);   y = point(1);
x0 = ellipse(1);    y0 = ellipse(2);
a = ellipse(4);     b = ellipse(3);
alpha = ellipse(5);
Q = [cos(alpha), sin(alpha); -sin(alpha), cos(alpha)];  % inverse of rotation matrix

% E(x',y') = abs(x'^2/a^2 + y'^2/b^2 - 1)
v = Q*[x-x0;y-y0];
v(1) = v(1)/a;
v(2) = v(2)/b;
error = abs(1-v'*v);

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
clear;
clc;
close all;

H = [-6 -6 -7 0 7 6 6 -3 -3 0 0 -6; -7 2 1 8 1 2 -7 -7 -2 -2 -7 -7];

[R, ~] = ch_rotation_2d(H, deg2rad(-90));
H = R*H;

x = H(1,:)'; y = H(2,:)';

axis('square');
axis equal
plot(x, y, 'o', x, y, '-');

author: yandld


text:
function [ax, data] = fgplvmFieldPlot(model, YLbls);

% FGPLVMFIELDPLOT 2-D field plot of the dynamics.

% FGPLVM

if isempty(YLbls)
  symbol = [];
else
  symbol = getSymbols(size(YLbls,2));
end

x1 = linspace(min(model.X(:, 1))*1.1, max(model.X(:, 1))*1.1, 150);
x2 = linspace(min(model.X(:, 2))*1.1, max(model.X(:, 2))*1.1, 150);
[X1, X2] = meshgrid(x1, x2);
XTest = [X1(:), X2(:)];

fhandle = str2func([model.type 'PosteriorMeanVar']);
if str2num(version('-release'))>13
  [mu, varsigma] = fhandle(model, XTest);
else
  [mu, varsigma] = feval(fhandle, model, XTest);
end
figure(1)
clf
% Create the plot for the data
clf
ax = axes('position', [0.05 0.05 0.9 0.9]);
hold on

C = log10(reshape(1./varsigma(:, 1), size(X1)));
C = C - min(min(C));
C = C/max(max(C));
C = round(C*63);
image(x1, x2, C);

% [c, h] = contourf(X1, X2, log10(reshape(1./varsigma(:, 1), size(X1))), 128);
% shading flat
colormap gray;
%colorbar
data = lvmtwoDPlot(model.X, YLbls, symbol);
xLim = [min(XTest(:, 1)) max(XTest(:, 1))];
yLim = [min(XTest(:, 2)) max(XTest(:, 2))];
set(ax, 'xLim', xLim);
set(ax, 'yLim', yLim);

set(ax, 'fontname', 'arial');
set(ax, 'fontsize', 20);

function returnVal = lvmtwoDPlot(X, label, symbol)

% LVMTWODPLOT Helper function for plotting the labels in 2-D.

returnVal = [];

if ~isempty(label)
  for i = 1:size(X, 1)
    labelNo = find(label(i, :));
    try
      returnVal = [returnVal; plot(X(i, 1), X(i, 2), symbol{labelNo})];
    catch
      if strcmp(lasterr, 'Index exceeds matrix dimensions.')
	error(['Only ' num2str(length(symbol)) ' labels supported (it''s easy to add more!)'])
      end
    end
  end
else
  returnVal = plot(X(:, 1), X(:, 2), 'rx');
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function Istitched = vertical_stitch(I,height)
% I must be within 0 and 1
Istitched = I;
while size(Istitched,1)<height
    window_size = round(size(I,1)/4);
    up = I((size(I,1)-window_size+1):size(I,1),:);
    down = I(1:window_size,:);
    aup = zeros(window_size, size(up,2));
    adown = zeros(window_size, size(up,2));
    for i=1:window_size
        aup(i,:) = up(i,:)*(1-i/window_size);
        adown(i,:) = down(i,:)*i/window_size;
    end
    Istitched = [Istitched(1:(size(Istitched,1)-window_size),:);aup+adown;Istitched((window_size+1):size(Istitched,1),:)];
end

Istitched = Istitched(1:height,:);

author: fumin

repo: pencil

save_path: github-repos/MATLAB/fumin-pencil


text:
function files = genfiles(path,name)
% recurisviely look for files in a path
% example usage:
% files = genfiles('~/Downloads','*.jpg');

folders = genpathsplit(path);

files = [];

for folderID=1:length(folders)
    cfiles = dir(fullfile(folders{folderID},name));
    for fileID=1:length(cfiles)
        files{end+1} = fullfile(folders{folderID},cfiles(fileID).name);
    end
end

author: jianxiongxiao

repo: ProfXkit


text:
function [mGal] = cmps22mGal(cmps2)
% Convert acceleration from centimeters per square centimeter to milligalileos
% Chad A. Greene 2012
mGal = cmps2*1e+3;
end

author: Sable


text:
function n = double(W, option)
   %DOUBLE returns a waveform's data as a double type
   %   N = double(waveform  [, option] )
   %   this is just a fancy way of saying get(W,'data')
   %   but can handle vectors of W. The data lengths of the vectors do not
   %   need to be the same.  By default it will zero-pad the end of all short
   %   vectors.  Using the OPTION, you can specify a NaN padding instead.
   %
   %   Input Arguments
   %       WAVEFORM: a waveform object   assumed as 1xN DIMENSIONAL
   %               (If not, it will be treated as though it were 1xN!)
   %       OPTION: Handle to any function that creates an array of a
   %               dictated size.  The chosen function should work like:
   %               A = (nRows,nCols); and return a double array of the
   %               appropriate size. ex. @ones, @zeros, @nan.
   %
   %   Output Arguments
   %       N: the columns of data
   %
   %   example:
   %       % let W(1) be a waveform with data [0 1 2 5 3]
   %       % let W(2) be a waveform with data [2 5]
   %
   %       Y = double(W);
   %       % Y is now [[ 0 2; 1 5; 2 0; 5 0; 3 0]
   %
   %       Z = double(W,@nan); pads with NaN instead of zeros
   %       % Z is now [ 0 2; 1 5; 2 NaN; 5 NaN; 3 NaN]
   %
   % See also WAVEFORM/GET, NAN

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   if numel(W) <= 1
      n = W.data; %  assumption --> data is already a column of double!
      return
   end

   if nargin == 1
      createDefaultArray = @zeros;
   else
      if ischar(option)
         createDefaultArray = str2func(option);
      elseif isa(option,'function_handle')
         createDefaultArray = option;
      else
         error('waveform:double:invalidOption',...
            '''option'' must be either a function handle or function name.');
      end
   end

   m = get(W,'data_length');
   if isempty(m)
      n = [];
      return
   end

   m=m(:);
   if all(m == m(1))
      n = [W(:).data]; % easy-peazy since all data is of same length.
   else
      n = createDefaultArray(max(m),numel(W));
      for x=1:numel(W)
         n(1:m(x),x) = W(x).data; % assumes that data is a column of double!
      end
   end
end

author: geoscience-community-codes

repo: GISMO


text:
function node_rhs = rhs ( node_num, node_xy )

%*****************************************************************************80
%
%% RHS gives the right-hand side of the differential equation.
%
%  Discussion:
%
%    This routine is set up for the L-shaped region, with exact solution
%    U = X^2 + Y^2.  Hence, the right hand side of the equation is
%    exactly -4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM),
%    the coordinates of the points.
%
%    Output, real NODE_RHS(NODE_NUM,1), the value of the
%    right hand side function at the points.
%
  node_rhs(1:node_num,1) = -4.0                     ...
                       + node_xy(1,1:node_num).^2 ...
                       + node_xy(2,1:node_num).^2;

  return
end

author: johannesgerer


text:
% [LEV,IND] = spyrBand(PYR,INDICES,LEVEL,BAND)
%
% Access a band from a steerable pyramid.
%
%   LEVEL indicates the scale (finest = 1, coarsest = spyrHt(INDICES)).
%
%   BAND (optional, default=1) indicates which subband
%     (1 = vertical, rest proceeding anti-clockwise).

% Eero Simoncelli, 6/96.

function res =  spyrBand(pyr,pind,level,band)

if (exist('level') ~= 1)
  level = 1;
end

if (exist('band') ~= 1)
  band = 1;
end

nbands = spyrNumBands(pind);
if ((band > nbands) | (band < 1))
  error(sprintf('Bad band number (%d) should be in range [1,%d].', band, nbands));
end

maxLev = spyrHt(pind);
if ((level > maxLev) | (level < 1))
  error(sprintf('Bad level number (%d), should be in range [1,%d].', level, maxLev));
end

firstband = 1 + band + nbands*(level-1);
res  = pyrBand(pyr, pind, firstband);

repo: SelfExSR


text:
%% HEATED_PLATE_FSU runs HEATED_PLATE_FUN using fsuClusterMatlab.
%
%  Discussion:
%
%    The arguments to fsuClusterMatlab mean:
%
%    [] allows us to specify an output directory;
%    [] allows us to specify switches to the MOAB queue handler, such as more time.
%    'p" indicates this is a "parfor" job;
%    'w' indicates that the MATLAB session should wait for completion;
%    4 is the number of processors requested;
%    @heated_plate_fun indicates the function we want to execute;
%    { tol, filename } contains input arguments.
%
%    The function "heated_plate_fun" must correspond to a MATLAB M-file
%    "heated_plate_fun.m" and that file must be in MATLAB's path.
%
%    The easiest way to ensure that the file is in MATLAB's path is to
%    create a subdirectory called "matlab" immediately below your home directory,
%    and place the function file there (as well as any other files needed by
%    the job).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 December 2010
%
%  Author:
%
%    John Burkardt
%
  tol = 0.001;
  filename = 'solution.txt';

  results = fsuClusterMatlab ( [], [], 'p', 'w', 4, @heated_plate_fun, { tol, filename } )

repo: jburkardt-m


text:
function [ff, pp, fi] = uiget1Dfile()
% open a gui for 1D file selection
% returns the same output as uigetfile

[ff, pp, fi] = uigetfile( {'*.1D;*.1D.dset' , ...
                         'AFNI 1D Files (*.1D, *.1D.dset)';},...
                         'Select 1D file' );
return

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function showweb(action)
    %  M file that starts a Browser displaying
    %  the HTML version of the ZMAP Users Guide
    %
    %  Stefan Wiemer   6/96

    report_this_filefun();
    disp('Attempting to open browser - please be patient...');

    switch(action)

        case 'new'
            file_to_manually_open = 'help/IntrotoZMAP6.htm';
            page_location = which('IntrotoZMAP6.htm');

        case 'stress'
            file_to_manually_open = 'help/stressinversions.pdf';
            page_location = which('stressinversions.htm');

        case 'data'
            %FIXME data no longer exists here.
            file_to_manually_open = 'help/onlinedata.htm';
            page_location = ' http://seismo.ethz.ch/staff/stefan/zmap6/help/onlinedata.htm ';

        case 'fractal'
            file_to_manually_open='help/FDH0.htm';
            page_location = which('FDH0.htm');

        case 'explproba'
            file_to_manually_open='help/explproba.htm';
            page_location = which('explproba.htm');

        case '3dgrids'
            file_to_manually_open='help/3dgrid.htm'; %was 3dgrids, but that doesn't existt
            page_location = which('3dgrid.htm');

        case 'topo'
            file_to_manually_open='help/plottopo.htm';
            page_location = which('plottopo.htm');
        otherwise
            file_to_manually_open=[];
            page_location = [];
    end

    if isempty(page_location)
        return
    end

    try
        % on mac, at least, page_location need not be preceded by 'file:'
        [status] = web(page_location,'-browser'); % open externally
        if ~(status==0)
            web(page_location); % open in MatLab's browser
        end
    catch
        disp(['Error opening browser  ... please open the file ' file_to_manually_open ' manually']);
    end
end

author: CelsoReyes


text:
function F=HessTaylor(deltaT,dadx,d2adx,method)
%%HESSTAYLOR Simulate a nonlinear continuous-time random process specified
%            by the Langevin equation forward in time by a step-size of
%            deltaT using a Taylor scheme.
%
%INPUTS: deltaT The size of the single step over which to generate the
%               state transition matrix.
%          dadx A xDimXxDim matrix of the derivative of the drift function
%               with respect to the state at state xCur and time curT.
%         d2adx A xDimXxDimXxDim matrix of the second derivative of the
%               drift function with respect to the state at state xCur and
%               time curT. The value at point (m,k,l) represents
%               d2a(m)/dx(k)dx(l). If not provided, this is assumed to be
%               zero.
%        method Set to 0 for the shorter Euler-Maruyama expansion,
%               otherwise will the Taylor expansion will be used (default).
%
%OUTPUTS: H The Hessian matrix of a nonlinear continuous-time random
%           process specified by the Langevin equation forward in time by a
%           step-size of deltaT using a strong Taylor scheme.
%
%The second derivative state prediction matrix is derived from the
%stochastic order 1.5 Taylor scheme described in 10.4 of [1].
%
%REFERENCES:
%[1] P. E. Kloeden and E. Platen, Numerical Solution of Stochastic
%    Differential Equations. Berlin: Springer, 1999.
%
%April 2015 David Karnick, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xDim=size(dadx,1);
if(nargin<4||isempty(method))
    method=1;
end

if(method==0)
    F=deltaT*d2adx;
else
    %Assumes 3rd derivative is zero
    term1=zeros(xDim,xDim,xDim);
    term2=zeros(xDim,xDim,xDim);
    term3=zeros(xDim,xDim,xDim);
    for m=1:xDim
        term1(:,:,m)=d2adx(:,:,m)*dadx;
        term2(:,m,:)=d2adx(:,:,m)*dadx;
        term3(:,:,m)=dadx*d2adx(:,:,m);
    end

    F=deltaT*d2adx+(deltaT^2/2)*(term1+term2+term3);
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function i4_abs_test ( )

%*****************************************************************************80
%
%% I4_ABS_TEST tests I4_ABS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2015
%
%  Author:
%
%    John Burkardt
%
  i4_lo = -100;
  i4_hi = +100;
  test_num = 10;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4_ABS_TEST\n' );
  fprintf ( 1, '  I4_ABS returns the absolute value of an I4.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      A         B=I4_ABS(A)\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
    [ a, seed ] = i4_uniform_ab ( i4_lo, i4_hi, seed );
    b = i4_abs ( a );
    fprintf ( 1, '  %10d  %10d\n', a, b );
  end

  return
end


text:
clear all; close all;
I=imread('cameraman.tif');
fun1=@dct2;
J1=blkproc(I, [8 8], fun1);
fun2=@(x) std2(x)*ones(size(x));
J2=blkproc(I, [8 8], fun2);
figure;
subplot(121);
imagesc(J1);
subplot(122);
imagesc(J2);
colormap gray;

author: HuangCongQing