text:
function [cov] = mne_pick_channels_cov(orig,include,exclude)
%
% [cov] = mne_pick_channels_cov(orig,include,exclude)
%
% Pick desired channels from a covariance matrix
%
% orig      - The original covariance matrix
% include   - Channels to include (if empty, include all available)
% exclude   - Channels to exclude (if empty, do not exclude any)
%
%

%
%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.9  2006/11/21 12:53:49  msh
%   Fixed error in picking
%
%   Revision 1.8  2006/06/29 22:12:28  msh
%   Fixed errors in channel picking
%
%   Revision 1.7  2006/05/03 18:53:05  msh
%   Approaching Matlab 6.5 backward compatibility
%
%   Revision 1.6  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.5  2006/04/18 20:44:46  msh
%   Added reading of forward solution.
%   Use length instead of size when appropriate
%
%   Revision 1.4  2006/04/15 12:21:00  msh
%   Several small improvements
%
%   Revision 1.3  2006/04/14 15:49:49  msh
%   Improved the channel selection code and added ch_names to measurement info.
%
%   Revision 1.2  2006/04/14 00:45:42  msh
%   Added channel picking and fiff_invert_transform
%
%   Revision 1.1  2006/04/13 17:05:45  msh
%   Added reading of bad channels to fiff_read_meas_info.m
%   Added mne_pick_channels_cov.m
%
%

me='MNE:mne_pick_channels_cov';

if nargin == 1
    cov = orig;
    if isempty(cov.eig) || isempty(cov.eigvec)
        decompose_eigen;
    end
    return;
elseif nargin == 2
    exclude = [];
elseif nargin ~= 3
    error(me,'Incorrect number of arguments');
end

if isempty(include) && isempty(exclude)
    cov = orig;
    if isempty(cov.eig) || isempty(cov.eigvec)
        decompose_eigen;
    end
    return;
end

if isempty(orig.names)
    error(me,'Cannot pick from a covariance matrix without channel names');
end

cov  = orig;
%
%   First do the channels to be included
%
sel = fiff_pick_channels(cov.names,include,exclude);
if isempty(sel)
   error(me,'Nothing remains after picking');
end
%
%   Select the desired stuff
%
if cov.diag
   cov.data = cov.data(sel);
else
   cov.data = cov.data(:,sel);
   cov.data = cov.data(sel,:);
end
for p = 1:size(sel,2)
   names{p} = cov.names{sel(p)};
end
cov.names = names;
cov.dim   = length(cov.names);
%
%   Eigenvalues and vectors are no longer valid
%
decompose_eigen;
%
return;

    function decompose_eigen
        if cov.diag
            cov.eig    = cov.data;
            cov.eigvec = eye(cov.dim);
        else
            [ cov.eigvec, cov.eig ] = eig(cov.data);
            %
            %   We use the convention that rows of eigvec are the
            %   eigenvectors
            %
            cov.eigvec = cov.eigvec';
            cov.eig    = diag(cov.eig);
        end
    end

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = process_ft_sourceanalysis_dics(varargin )
% PROCESS_FT_SOURCEANALYSIS Call FieldTrip function ft_sourceanalysis (DICS)

% @=============================================================================
% 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.
% =============================================================================@
%
% Authors: Vahab YoussofZadeh, Francois Tadel, 2021

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'FieldTrip: ft_sourceanalysis (DICS)';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Sources';
    sProcess.Index       = 357;
    sProcess.Description = 'https://github.com/vyoussofzadeh/DICS-beamformer-for-Brainstorm';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data'};
    sProcess.OutputTypes = {'data'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;

    % Option: Sensors selection
    sProcess.options.sensortype.Comment = 'Sensor type:';
    sProcess.options.sensortype.Type    = 'combobox_label';
    sProcess.options.sensortype.Value   = {'MEG', {'MEG', 'MEG GRAD', 'MEG MAG', 'EEG', 'SEEG', 'ECOG'; ...
                                                   'MEG', 'MEG GRAD', 'MEG MAG', 'EEG', 'SEEG', 'ECOG'}};
    % Label: Time
    sProcess.options.label1.Comment = '
Time of interest:';
    sProcess.options.label1.Type    = 'label';
    % Active time window
    sProcess.options.poststim.Comment = 'Active (post-stim):';
    sProcess.options.poststim.Type    = 'poststim';
    sProcess.options.poststim.Value   = [];
    % Baseline time window
    sProcess.options.baseline.Comment = 'Baseline (pre-stim):';
    sProcess.options.baseline.Type    = 'baseline';
    sProcess.options.baseline.Value   = [];

    % Label: Frequency
    sProcess.options.label2.Comment = '
Frequency of interset:';
    sProcess.options.label2.Type    = 'label';
    % Enter the FOI in the data in Hz, eg, 22:
    sProcess.options.foi.Comment = 'FOI:';
    sProcess.options.foi.Type    = 'value';
    sProcess.options.foi.Value   = {18, 'Hz', 0};
    % Enter the Tapering frequency in the data in Hz, eg, 4:
    sProcess.options.tpr.Comment = 'Tapering freq:';
    sProcess.options.tpr.Type    = 'value';
    sProcess.options.tpr.Value   = {4, 'Hz', 0};

    % Label: Contrast
    sProcess.options.label4.Comment = '
Contrasting analysis:'; % Contrast between pre and post, across trials
    sProcess.options.label4.Type    = 'label';
    % Contrast
    sProcess.options.method.Comment = {'Subtraction (post-pre)', 'Permutation-stats (post-vs-pre)', 'Contrasting:'; ...
                                       'subtraction', 'permutation', ''};
    sProcess.options.method.Type    = 'radio_linelabel';
    sProcess.options.method.Value   = 'subtraction';
    % ERDS
    sProcess.options.erds.Comment = {'ERD', 'ERS', 'both', 'ERD/S effects:'; ...
                                     'erd', 'ers', 'both', ''};
    sProcess.options.erds.Type    = 'radio_linelabel';
    sProcess.options.erds.Value   = 'erd';
    % Effects
    sProcess.options.effect.Comment = {'abs', 'raw', 'Absolute/raw value of the contrast:'; ...
                                       'abs', 'raw', ''};
    sProcess.options.effect.Type    = 'radio_linelabel';
    sProcess.options.effect.Value   = 'abs';

    % Label: TFR
    sProcess.options.label3.Comment = '
Time-freq response (used for data inspection):';
    sProcess.options.label3.Type    = 'label';
    % Max frequency
    sProcess.options.maxfreq.Comment = 'Max frequency (for TFR calculation):';
    sProcess.options.maxfreq.Type    = 'value';
    sProcess.options.maxfreq.Value   = {40, 'Hz', 0};
    % Plot time-freq
    sProcess.options.showtfr.Comment = 'Plot time-frequency figure for visual';
    sProcess.options.showtfr.Type    = 'checkbox';
    sProcess.options.showtfr.Value   = 1;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    OutputFiles = {};
    % Initialize FieldTrip
    [isInstalled, errMsg] = bst_plugin('Install', 'fieldtrip');
    if ~isInstalled
        bst_report('Error', sProcess, [], errMsg);
        return;
    end

    % ===== GET OPTIONS =====
    % Inverse options
    Method   = sProcess.options.method.Value;
    Modality = sProcess.options.sensortype.Value{1};
    ShowTfr  = sProcess.options.showtfr.Value;
    MaxFreq  = sProcess.options.maxfreq.Value{1};
    Baseline = sProcess.options.baseline.Value{1};
    PostStim = sProcess.options.poststim.Value{1};
    FOI      = sProcess.options.foi.Value{1};
    TprFreq  = sProcess.options.tpr.Value{1};
    % Progress bar
    bst_progress('start', 'ft_sourceanalysis', 'Loading input files...', 0, 2*length(sInputs));

    % ===== LOAD: CHANNEL FILE =====
    bst_progress('text', 'Loading input files...');

    % Load channel file
    ChannelMat = in_bst_channel(sInputs(1).ChannelFile);
    % Get selected sensors
    iChannels = channel_find(ChannelMat.Channel, Modality);
    if isempty(iChannels)
        bst_report('Error', sProcess, sInputs, ['Channels "' Modality '" not found in channel file.']);
        return;
    end

    % ===== LOAD: BAD CHANNELS =====
    % Load bad channels from all the input files
    isChannelGood = [];
    for iInput = 1:length(sInputs)
        DataFile = sInputs(1).FileName;
        DataMat = load(file_fullpath(DataFile), 'ChannelFlag');
        if isempty(isChannelGood)
            isChannelGood = (DataMat.ChannelFlag == 1);
        elseif (length(DataMat.ChannelFlag) ~= length(isChannelGood)) || (length(DataMat.ChannelFlag) ~= length(ChannelMat.Channel))
            bst_report('Error', sProcess, sInputs, 'All the input files must have the same number of channels.');
            return;
        else
            isChannelGood = isChannelGood & (DataMat.ChannelFlag == 1);
        end
    end
    % Remove bad channels
    iChannelsData = intersect(iChannels, find(isChannelGood'));
    % Error: All channels tagged as bad
    if isempty(iChannelsData)
        bst_report('Error', sProcess, sInputs, 'All the selected channels are tagged as bad.');
        return;
    elseif any(~isChannelGood)
        bst_report('Info', sProcess, sInputs, ['Found ' num2str(length(find(~isChannelGood))) ' bad channels: ', sprintf('%s ', ChannelMat.Channel(find(~isChannelGood)).Name)]);
    end

    % ===== LOAD: HEADMODEL =====
    % Get the study
    sStudyChan = bst_get('ChannelFile', sInputs(1).ChannelFile);
    % Error if there is no head model available
    if isempty(sStudyChan.iHeadModel)
        bst_report('Error', sProcess, [], ['No head model available in folder: ' bst_fileparts(sStudyChan.FileName)]);
        return;
    end
    % Load head model
    HeadModelFile = sStudyChan.HeadModel(sStudyChan.iHeadModel).FileName;
    HeadModelMat = in_bst_headmodel(HeadModelFile);
    % Convert head model to FieldTrip format
    [ftHeadmodel, ftLeadfield, iChannelsData] = out_fieldtrip_headmodel(HeadModelMat, ChannelMat, iChannelsData, 1);

    % ===== LOAD: DATA =====
    % Template FieldTrip structure for all trials
    ftData = out_fieldtrip_data(sInputs(1).FileName, ChannelMat, iChannelsData, 1);
    ftData.trial = cell(1,length(sInputs));
    ftData.time = cell(1,length(sInputs));
    % Load all the trials
    AllChannelFiles = unique({sInputs.ChannelFile});
    iChanInputs = find(ismember({sInputs.ChannelFile}, AllChannelFiles{1}));
    for iInput = 1:length(sInputs)
        DataFile = sInputs(iChanInputs(iInput)).FileName;
        DataMat = in_bst_data(DataFile);
        ftData.trial{iInput} = DataMat.F(iChannelsData,:);
        ftData.time{iInput} = DataMat.Time;
    end

    % ===== FIELDTRIP: ft_freqanalysis =====
    bst_progress('text', 'Calling FieldTrip function: ft_freqanalysis...');
    % Compute tfr-decomposition
    cfg = [];
    cfg.output     = 'pow';
    cfg.channel    = 'all';
    cfg.method     = 'mtmconvol';
    cfg.taper      = 'hanning';
    cfg.foi        = 1:2:MaxFreq;
    cfg.keeptrials = 'yes';
    cfg.t_ftimwin  = 3 ./ cfg.foi;
    cfg.tapsmofrq  = 0.8 * cfg.foi;
    cfg.toi        = Baseline(1):0.05:PostStim(2);
    tfr            = ft_freqanalysis(cfg, ftData);
    % Plot TFR
    cfg = [];
    cfg.fmax = MaxFreq;
    cfg.toi = [Baseline(1), PostStim(2)];
    cfg.baselinetime = Baseline;
    cfg.plotflag = ShowTfr;
    [time_of_interest,freq_of_interest] = do_tfr_plot(cfg, tfr);
    disp(['Global max: time:',num2str(time_of_interest),'sec']);
    disp(['Global max: freq:',num2str(freq_of_interest),'Hz']);

    % ===== FIELDTRIP: EPOCHING =====
    % Baseline
    cfg = [];
    cfg.toilim = Baseline;
    ep_data.bsl = ft_redefinetrial(cfg, ftData);
    % Post-stim
    cfg.toilim = PostStim;
    ep_data.pst = ft_redefinetrial(cfg, ftData);
    % Baseline + Post-stim
    cfg = [];
    ep_data.app = ft_appenddata(cfg, ep_data.bsl, ep_data.pst);

    % ===== TEMPORARY FOLDER =====
    % Output folder (for figures and FieldTrip structures)
    TmpDir = bst_get('BrainstormTmpDir', 0, 'dics');

    % ===== FIELDTRIP: SPECTRAL ANALYSIS =====
    cfg_main = [];
    cfg_main.fmax = MaxFreq;
    switch sProcess.options.sensortype.Value{1}
        case {'EEG', 'SEEG', 'ECOG'}
            cfg_main.sens = ftData.elec;
        case {'MEG', 'MEG GRAD', 'MEG MAG'}
            cfg_main.sens = ftData.grad;
    end
    cfg_main.outputdir = TmpDir;
    cfg_main.freq_of_interest  = freq_of_interest; % Hz

    cfg = [];
    cfg.foilim    = [2 cfg_main.fmax];
    cfg.tapsmofrq = 1;
    cfg.taper     = 'hanning';
    f_data.bsl = do_fft(cfg, ep_data.bsl); f_data.bsl.elec = cfg_main.sens;
    f_data.pst = do_fft(cfg, ep_data.pst); f_data.pst.elec = cfg_main.sens;

    % ===== FIELDTRIP: PSD SENSOR SPACE =====
    outputdir_dics = cfg_main.outputdir;
    if exist(outputdir_dics, 'file') == 0, mkdir(outputdir_dics), end

    f_sugg = round(cfg_main.freq_of_interest);
    disp(['Suggested by TFR: ', num2str(f_sugg),'(+-3Hz)']);
    disp(['Select foi,eg ,', num2str(f_sugg),':']);

    cfg = [];
    cfg.foilim = [FOI, FOI];
    if (FOI >= 4)
        cfg.taper = 'dpss';
        cfg.tapsmofrq = TprFreq;
    else
        cfg.taper = 'hanning';
        cfg.tapsmofrq = 1;
    end

    f_data.app = do_fft(cfg, ep_data.app); f_data.app.elec = cfg_main.sens;
    f_data.bsl = do_fft(cfg, ep_data.bsl); f_data.bsl.elec = cfg_main.sens;
    f_data.pst = do_fft(cfg, ep_data.pst); f_data.pst.elec = cfg_main.sens;

    % ===== SOURCE ANALYSIS =====
    switch Method
        case 'subtraction'
            cfg = [];
            cfg.method = 'dics';
            cfg.dics.lambda = '100%';
            cfg.sourcemodel  = ftLeadfield;
            cfg.frequency    = f_data.app.freq;
            cfg.headmodel = ftHeadmodel;
            cfg.dics.keepfilter = 'yes';
            cfg.dics.fixedori    = 'yes'; % project on axis of most variance using SVD
            sourceavg = ft_sourceanalysis(cfg, f_data.app);

            cfg = [];
            cfg.method = 'dics';
            cfg.dics.lambda = '0%';
            cfg.sourcemodel        = ftLeadfield;
            cfg.sourcemodel.filter = sourceavg.avg.filter;
            cfg.dics.fixedori    = 'yes'; % project on axis of most variance using SVD
            cfg.headmodel = ftHeadmodel;
            s_data.bsl      = ft_sourceanalysis(cfg, f_data.bsl);
            s_data.pst      = ft_sourceanalysis(cfg, f_data.pst);

            switch sProcess.options.erds.Value
                case 'erd'
                    cfg = [];
                    cfg.parameter = 'pow';
                    cfg.operation = 'log10(x1/x2)'; % sourceA divided by sourceB
                    source_diff_dics = ft_math(cfg, s_data.pst, s_data.bsl);
                    source_diff_dics.pow(isnan(source_diff_dics.pow))=0;
                    source_diff_dics.pow(source_diff_dics.pow>0)=0;
                case 'ers'
                    cfg = [];
                    cfg.parameter = 'pow';
                    cfg.operation = 'log10(x1/x2)'; % sourceA divided by sourceB
                    source_diff_dics = ft_math(cfg, s_data.pst, s_data.bsl);
                    source_diff_dics.pow(isnan(source_diff_dics.pow))=0;
                    source_diff_dics.pow(source_diff_dics.pow<0)=0;
                case 'both'
                    cfg = [];
                    cfg.parameter = 'pow';
                    cfg.operation = 'log10(x1/x2)'; % sourceA divided by sourceB
                    source_diff_dics = ft_math(cfg,s_data.pst,s_data.bsl);
                    source_diff_dics.pow(isnan(source_diff_dics.pow))=0;
            end

        case 'permutation'
            cfg = [];
            cfg.method = 'dics';
            cfg.dics.lambda = '100%';
            cfg.frequency    = f_data.app.freq;
            cfg.headmodel = ftHeadmodel;
            cfg.sourcemodel  = ftLeadfield;
            cfg.dics.keepfilter = 'yes';
            cfg.dics.fixedori    = 'yes'; % project on axis of most variance using SVD
            sourceavg = ft_sourceanalysis(cfg, f_data.app);

            cfg = [];
            cfg.method = 'dics';
            cfg.sourcemodel        = ftLeadfield;
            cfg.sourcemodel.filter = sourceavg.avg.filter;
            cfg.dics.fixedori    = 'yes'; % project on axis of most variance using SVD
            cfg.rawtrial = 'yes';
            cfg.headmodel = ftHeadmodel;
            s_data.bsl      = ft_sourceanalysis(cfg, f_data.bsl);
            s_data.pst      = ft_sourceanalysis(cfg, f_data.pst);

            stat = do_source_stat_montcarlo(s_data);

            tmp = stat.stat;
            tmp2 = zeros(size(stat.pos,1),1);
            tmp2(stat.inside) = tmp;

            stats1  = stat;
            stats1.stat =  tmp2;
            stats1.mask = stat.inside;
            stats2 = stats1;
            stats2.stat(stats2.stat>0)=0;
            stats2.stat(isnan(stats2.stat))=0;
    end

    % ===== SAVE RESULTS =====
    % === CREATE OUTPUT STRUCTURE ===
    bst_progress('text', 'Saving source file...');
    bst_progress('inc', 1);
    % Output study
    if (length(sInputs) == 1)
        iStudyOut = sInputs(1).iStudy;
        RefDataFile = sInputs(iChanInputs(iInput)).FileName;
    else
        [tmp, iStudyOut] = bst_process('GetOutputStudy', sProcess, sInputs);
        RefDataFile = [];
    end
    % Create structure
    ResultsMat = db_template('resultsmat');
    ResultsMat.ImagingKernel = [];
    switch Method
        case 'subtraction'
            switch sProcess.options.effect.Value
                case 'abs'
                    source_diff_dics.pow = abs((source_diff_dics.pow));
                case 'raw'
                    source_diff_dics.pow = source_diff_dics.pow;
            end
            ResultsMat.ImageGridAmp  = source_diff_dics.pow;
            ResultsMat.cfg           = source_diff_dics.cfg;

        case 'permutation'
            switch sProcess.options.effect.Value
                case 'abs'
                    stats2.stat = abs((stats2.stat));
                case 'raw'
                    stats2.stat = stats2.stat;
            end
            ResultsMat.ImageGridAmp  = stats2.stat;
            ResultsMat.cfg           = stat.cfg;
    end
    ResultsMat.nComponents   = 1;
    ResultsMat.Function      = Method;
    ResultsMat.Time          = 1;
    ResultsMat.DataFile      = RefDataFile;
    ResultsMat.HeadModelFile = HeadModelFile;
    ResultsMat.HeadModelType = HeadModelMat.HeadModelType;
    ResultsMat.ChannelFlag   = DataMat.ChannelFlag;
    ResultsMat.GoodChannel   = iChannelsData;
    ResultsMat.SurfaceFile   = HeadModelMat.SurfaceFile;
    ResultsMat.nAvg          = DataMat.nAvg;
    ResultsMat.Leff          = DataMat.Leff;
    ResultsMat.Comment       = ['DICS: ' Method, ' ',num2str(FOI),'Hz ', sprintf('%1.3fs-%1.3fs', PostStim)];
    switch lower(ResultsMat.HeadModelType)
        case 'volume'
            ResultsMat.GridLoc    = HeadModelMat.GridLoc;
        case 'surface'
            ResultsMat.GridLoc    = [];
        case 'mixed'
            ResultsMat.GridLoc    = HeadModelMat.GridLoc;
            ResultsMat.GridOrient = HeadModelMat.GridOrient;
    end
    ResultsMat = bst_history('add', ResultsMat, 'compute', ['ft_sourceanalysis: ' Method ' ' Modality ' ']);

    % === SAVE OUTPUT FILE ===
    % Output filename
    OutputDir = bst_fileparts(file_fullpath(DataFile));
    ResultFile = bst_process('GetNewFilename', OutputDir, ['results_', Method, '_', Modality]);
    % Save new file structure
    bst_save(ResultFile, ResultsMat, 'v6');

    % ===== REGISTER NEW FILE =====
    % Create new results structure
    newResult = db_template('results');
    newResult.Comment       = ResultsMat.Comment;
    newResult.FileName      = file_short(ResultFile);
    newResult.DataFile      = ResultsMat.DataFile;
    newResult.isLink        = 0;
    newResult.HeadModelType = ResultsMat.HeadModelType;
    % Get output study
    sStudyOut = bst_get('Study', iStudyOut);
    % Add new entry to the database
    iResult = length(sStudyOut.Result) + 1;
    sStudyOut.Result(iResult) = newResult;
    % Update Brainstorm database
    bst_set('Study', iStudyOut, sStudyOut);
    % Store output filename
    OutputFiles{end+1} = newResult.FileName;
    % Expand data node
    panel_protocols('SelectNode', [], newResult.FileName);

    % Delete the temporary files
    file_delete(TmpDir, 1, 1);
    % Save database
    db_save();
    % Hide progress bar
    bst_progress('stop');
end

%% ===== TIME-FREQUENCY =====
function [time_of_interest,freq_of_interest] = do_tfr_plot(cfg_main, tfr)
    % First compute the average over trials:
    cfg = [];
    freq_avg = ft_freqdescriptives(cfg, tfr);

    % And baseline-correct the average:
    cfg = [];
    cfg.baseline = cfg_main.baselinetime;
    cfg.baselinetype = 'db'; % Use decibel contrast here
    freq_avg_bsl = ft_freqbaseline(cfg, freq_avg);

    freq_avg_bsl.powspctrm(isnan(freq_avg_bsl.powspctrm))=0;
    meanpow = squeeze(mean(freq_avg_bsl.powspctrm, 1));

    tim_interp = linspace(cfg_main.toi(1), cfg_main.toi(2), 512);
    freq_interp = linspace(1, cfg_main.fmax, 512);

    % We need to make a full time/frequency grid of both the original and
    % interpolated coordinates. Matlab's meshgrid() does this for us:
    [tim_grid_orig, freq_grid_orig] = meshgrid(tfr.time, tfr.freq);
    [tim_grid_interp, freq_grid_interp] = meshgrid(tim_interp, freq_interp);

    % And interpolate:
    pow_interp = interp2(tim_grid_orig, freq_grid_orig, meanpow, tim_grid_interp, freq_grid_interp, 'spline');

    % while n==1
    pow_interp1  = pow_interp(50:end,50:end);
    tim_interp1  = tim_interp(50:end);
    freq_interp1 = freq_interp(50:end);

    [~,idx] = min(pow_interp1(:));
    [row,col] = ind2sub(size(pow_interp1),idx);

    time_of_interest = tim_interp1(col);
    freq_of_interest = freq_interp1(row);

    timind = nearest(tim_interp, time_of_interest);
    freqind = nearest(freq_interp, freq_of_interest);
    pow_at_toi = pow_interp(:,timind);
    pow_at_foi = pow_interp(freqind,:);

    % Plot figure
    if cfg_main.plotflag
        figure();
        ax_main  = axes('Position', [0.1 0.2 0.55 0.55]);
        ax_right = axes('Position', [0.7 0.2 0.1 0.55]);
        ax_top   = axes('Position', [0.1 0.8 0.55 0.1]);

        axes(ax_main);
        imagesc(tim_interp, freq_interp, pow_interp);
        % note we're storing a handle to the image im_main, needed later on
        xlim([cfg_main.toi(1), cfg_main.toi(2)]);
        axis xy;
        xlabel('Time (s)');
        ylabel('Frequency (Hz)');
        clim = max(abs(meanpow(:)));
        caxis([-clim clim]);
        % colormap(brewermap(256, '*RdYlBu'));
        hold on;
        plot(zeros(size(freq_interp)), freq_interp, 'k:');

        axes(ax_top);
        area(tim_interp, pow_at_foi, ...
            'EdgeColor', 'none', 'FaceColor', [0.5 0.5 0.5]);
        xlim([cfg_main.toi(1), cfg_main.toi(2)]);
        ylim([-clim clim]);
        box off;
        ax_top.XTickLabel = [];
        ylabel('Power (dB)');
        hold on;
        plot([0 0], [-clim clim], 'k:');

        axes(ax_right);
        area(freq_interp, pow_at_toi,...
            'EdgeColor', 'none', 'FaceColor', [0.5 0.5 0.5]);
        view([270 90]); % this rotates the plot
        ax_right.YDir = 'reverse';
        ylim([-clim clim]);
        box off;
        ax_right.XTickLabel = [];
        ylabel('Power (dB)');

        h = colorbar(ax_main, 'manual', 'Position', [0.85 0.2 0.05 0.55]);
        ylabel(h, 'Power vs baseline (dB)');

        % Main plot:
        axes(ax_main);
        plot(ones(size(freq_interp))*time_of_interest, freq_interp,...
            'Color', [0 0 0 0.1], 'LineWidth', 3);
        plot(tim_interp, ones(size(tim_interp))*freq_of_interest,...
            'Color', [0 0 0 0.1], 'LineWidth', 3);

        % Marginals:
        axes(ax_top);
        plot([time_of_interest time_of_interest], [0 clim],...
            'Color', [0 0 0 0.1], 'LineWidth', 3);
        axes(ax_right);
        hold on;
        plot([freq_of_interest freq_of_interest], [0 clim],...
            'Color', [0 0 0 0.1], 'LineWidth', 3);
    end
end

function [freq,ff, psd,tapsmofrq] = do_fft(cfg_mian, data)
    cfg              = [];
    cfg.method       = 'mtmfft';
    cfg.output       = 'fourier';
    cfg.keeptrials   = 'yes';
    cfg.foilim       = cfg_mian.foilim;
    cfg.tapsmofrq    = cfg_mian.tapsmofrq;
    cfg.taper        = cfg_mian.taper;
    cfg.pad          = 4;
    freq             = ft_freqanalysis(cfg, data);
    psd = squeeze(mean(mean(abs(freq.fourierspctrm),2),1));
    ff = linspace(1, cfg.foilim(2), length(psd));

    tapsmofrq = cfg.tapsmofrq;
end

function stat = do_source_stat_montcarlo(s_data)
    cfg = [];
    cfg.parameter        = 'pow';
    cfg.method           = 'montecarlo';
    cfg.statistic        = 'depsamplesT';
    cfg.correctm         = 'fdr';
    cfg.clusteralpha     = 0.001;
    cfg.tail             = 0;
    cfg.clustertail      = 0;
    cfg.alpha            = 0.05;
    cfg.numrandomization = 5000;

    ntrials                       = numel(s_data.bsl.trial);
    design                        = zeros(2,2*ntrials);
    design(1,1:ntrials)           = 1;
    design(1,ntrials+1:2*ntrials) = 2;
    design(2,1:ntrials)           = 1:ntrials;
    design(2,ntrials+1:2*ntrials) = 1:ntrials;

    cfg.design   = design;
    cfg.ivar     = 1;
    cfg.uvar     = 2;
    stat         = ft_sourcestatistics(cfg,s_data.pst,s_data.bsl);
end

author: brainstorm-tools


text:
% Compute y = log( sum(exp(x),2) ), the softmax in a numerically safe way by
%  subtracting the row maximum to avoid cancelation after taking the exp
%  the sum is done along the rows.
%
% Copyright (c) by Hannes Nickisch, 2013-10-16.

function [y,x] = logsumexp2(logx)
  N = size(logx,2); max_logx = max(logx,[],2);
  % we have all values in the log domain, and want to calculate a sum
  x = exp(logx-max_logx*ones(1,N));
  y = log(sum(x,2)) + max_logx;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% Euclidean projection on a rectangle
% Section 8.1.1, Boyd & Vandenberghe "Convex Optimization"
% Joelle Skaf - 10/07/05
%
% The projection of x0 on a rectangle C = {x | l <= x <= u} is given by
%           minimize || x - x0 ||^2
%               s.t.    l <= x <= u
% It is also given by P_C(x0)_k = l_k       if  x0_k <= l_k
%                                 x0_k      if  l_k <= x0_k <= u_k
%                                 u_k       if  x0_k >= u_k

% Input data: generate vectors l and u such that l < 0 < u
n  = 10;
l  = -rand(n,1);
u  = rand(n,1);
x0 = randn(n,1);

% Analytical solution
fprintf(1,'Computing the analytical solution ...');
pc_x0 = x0;
pc_x0(find(x0<=l)) = l(find(x0<=l));
pc_x0(find(x0>=u)) = u(find(x0>=u));
fprintf(1,'Done! \n');

% Solution via QP
fprintf(1,'Computing the optimal solution by solving a QP ...');

cvx_begin quiet
    variable x(n)
    minimize ( norm(x-x0) )
    x <= u;
    x >= l;
cvx_end

fprintf(1,'Done! \n');

% Verification
disp('-----------------------------------------------------------------');
disp('Verifying that the analytical solution and the solution obtained via QP are equal: ');
[pc_x0 x]

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function createVelocititesControlButtons(panel_x0,panel_y0,fig_h,panel_Width,panel_height)
%% An interface for velocities magitueds control, this interface is a part of the Virtual smartPad interface
% Copyright: Mohammad SAFEEA, 2018-July-27th

global virtualTeachPendantLinearVelocityMagnitude;
global virtualTeachPendantAngularVelocityMagnitude;
global virtualTeachPendantJointVelocityMagnitude;

virtualTeachPendantLinearVelocityMagnitude=5;
virtualTeachPendantAngularVelocityMagnitude=2;
virtualTeachPendantJointVelocityMagnitude=2;

initialVals={virtualTeachPendantLinearVelocityMagnitude;
    virtualTeachPendantAngularVelocityMagnitude;
    virtualTeachPendantJointVelocityMagnitude};

% Draw a panel, the container of interface
pnl_h=uipanel(fig_h,'Title','Motion velocity','FontSize',12,...
            'Units','pixels');
pnl_h.Position=[panel_x0 panel_y0 panel_Width panel_height];

interMargin=2;
basic_width=(panel_Width-5*interMargin)/10;

upperMargin=60;
interCentersVerticalDistance=(panel_height-upperMargin)/2.7;

text_height=30;
labelsHeight=25;
control_height=interCentersVerticalDistance*.75;

centegradeSymbol=char(176);
st1Temp=['EEF ang... [',centegradeSymbol,'/sec]'];
st2Temp=['Joint vel [',centegradeSymbol,'/sec]'];
labels={'EEF lin... [mm/sec]';st1Temp;st2Temp};

    for i=1:3
        yCenter=panel_height-upperMargin-interCentersVerticalDistance*(i-1);
        % joint label
        lbl=uicontrol(pnl_h,'Style','text',...
            'Units','pixels',...
            'String',labels{i});
        cWidth=5*basic_width;
        cHeight=labelsHeight;
        x0=interMargin;
        y0=yCenter-control_height/2;
        lbl.Position=[x0 y0 cWidth cHeight];
        set(lbl,'HorizontalAlignment','left')
        % minus button
        minusButHandle=uicontrol(pnl_h,'Style', 'pushbutton', 'String', '-',...
                'Units','pixels');
        x0=x0+cWidth+interMargin ;
        cWidth=basic_width;
        cHeight=control_height;
        y0=yCenter-cHeight/2;
        minusButHandle.Position=[x0 y0 cWidth cHeight];
        % text of velocity magnitude
        str=sprintf('%0.1f',initialVals{i});
        textHandle=uicontrol(pnl_h,'Style', 'edit', 'String', str,...
                'Units','pixels');
        x0=x0+cWidth+interMargin ;
        cWidth=2*basic_width;
        cHeight=text_height;
        y0=yCenter-cHeight/2;
        textHandle.Position=[x0 y0 cWidth cHeight];
        % plus button
        plusButHandle=uicontrol(pnl_h,'Style', 'pushbutton', 'String', '+',...
                'Units','pixels');
        x0=x0+cWidth+interMargin ;
        cWidth=basic_width;
        cHeight=control_height;
        y0=yCenter-cHeight/2;
        plusButHandle.Position=[x0 y0 cWidth cHeight];
        % set callbacks of buttons
        if i==1
            fun1=@(object,event)velLinearCallBack(object,event,textHandle,+1);
            fun2=@(object,event)velLinearCallBack(object,event,textHandle,-1);
        elseif i==2
            fun1=@(object,event)velAngularCallBack(object,event,textHandle,+0.2);
            fun2=@(object,event)velAngularCallBack(object,event,textHandle,-0.2);
        elseif i==3
            fun1=@(object,event)jointVelCallBack(object,event,textHandle,+0.2);
            fun2=@(object,event)jointVelCallBack(object,event,textHandle,-0.2);
        end
        set(minusButHandle,'Callback', fun2);
        set(plusButHandle,'Callback', fun1);

    end
end

function velLinearCallBack(object,event,textHandle,increment)
global virtualTeachPendantLinearVelocityMagnitude;
virtualTeachPendantLinearVelocityMagnitude=virtualTeachPendantLinearVelocityMagnitude+increment;
maxLimit=10;
minLimit=1;
if virtualTeachPendantLinearVelocityMagnitude>maxLimit
    virtualTeachPendantLinearVelocityMagnitude=maxLimit;
end
if virtualTeachPendantLinearVelocityMagnitude<minLimit
    virtualTeachPendantLinearVelocityMagnitude=minLimit;
end
% Update the textbox
str=sprintf('%0.1f',virtualTeachPendantLinearVelocityMagnitude);
set(textHandle,'String',str);
end

function velAngularCallBack(object,event,textHandle,increment)
global virtualTeachPendantAngularVelocityMagnitude;
virtualTeachPendantAngularVelocityMagnitude=virtualTeachPendantAngularVelocityMagnitude+increment;
maxLimit=10;
minLimit=1;
if virtualTeachPendantAngularVelocityMagnitude>maxLimit
    virtualTeachPendantAngularVelocityMagnitude=maxLimit;
end
if virtualTeachPendantAngularVelocityMagnitude<minLimit
    virtualTeachPendantAngularVelocityMagnitude=minLimit;
end
% Update the textbox
str=sprintf('%0.1f',virtualTeachPendantAngularVelocityMagnitude);
set(textHandle,'String',str);
end

function jointVelCallBack(object,event,textHandle,increment)
global virtualTeachPendantJointVelocityMagnitude;
virtualTeachPendantJointVelocityMagnitude=virtualTeachPendantJointVelocityMagnitude+increment;
maxLimit=10;
minLimit=1;
if virtualTeachPendantJointVelocityMagnitude>maxLimit
    virtualTeachPendantJointVelocityMagnitude=maxLimit;
end
if virtualTeachPendantJointVelocityMagnitude<minLimit
    virtualTeachPendantJointVelocityMagnitude=minLimit;
end
% Update the textbox
str=sprintf('%0.1f',virtualTeachPendantJointVelocityMagnitude);
set(textHandle,'String',str);
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function value = ch_is_printable ( ch )

%*****************************************************************************80
%
%% CH_IS_PRINTABLE determines if a character is printable.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character CH, a character to check.
%
%    Output, logical VALUE is TRUE if C is a printable character.
%
  if ( 32 <= ch & ch <= 127 )
    value = 1;
  else
    value = 0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
tracker_label = 'Demo_py';

% Note: be carefull for double backslashes on Windows
tracker_command = generate_python_command('python_static_rgbd', ...
    {'<path-to-toolkit-dir>\\tracker\\examples\\python', ...  % tracker source and vot.py are here
    '<path-to-trax-dir>\\support\\python'});

tracker_interpreter = 'python';

tracker_linkpath = {'<path-to-trax-build-dir>\bin', ...
    '<path-to-trax-build-dir>\lib'};

author: votchallenge

repo: toolkit-legacy


text:
function rv = mostver(varargin)
%MOSTVER  Prints or returns MOST version info for current installation.
%   V = MOSTVER returns the current MOST version number.
%   V = MOSTVER('all') returns a struct with the fields Name, Version,
%   Release and Date (all strings). Calling MOSTVER without assigning the
%   return value prints the version and release date of the current
%   installation of MOST.
%
%   See also MPVER.

%   MOST
%   Copyright (c) 2010-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MOST.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/most for more info.

v = struct( 'Name',     'MOST', ...
            'Version',  '1.2', ...
            'Release',  '', ...
            'Date',     '13-Dec-2022' );
if nargout > 0
    if nargin > 0
        rv = v;
    else
        rv = v.Version;
    end
else
    fprintf('%-22s Version %-9s  %11s\n', v.Name, v.Version, v.Date);
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function stirling1_values_test ( )

%*****************************************************************************80
%
%% STIRLING1_VALUES_TEST demonstrates the use of STIRLING1_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'STIRLING1_VALUES_TEST:\n' );
  fprintf ( 1, '  STIRLING1_VALUES returns values of \n' );
  fprintf ( 1, '  the Stirling numbers of the first kind.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      N      M        S1\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, m, n, s1 ] = stirling1_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %4d  %4d  %10d\n', m, n, s1 );

  end

  return
end


text:
function [ u, it_num ] = monogrid_poisson_1d ( n, a, b, ua, ub, force, exact )

%*****************************************************************************80
%
%% MONOGRID_POISSON_1D solves a 1D PDE, using the Gauss-Seidel method.
%
%  Discussion:
%
%    This routine solves a 1D boundary value problem of the form
%
%      - U''(X) = F(X) for A < X < B,
%
%    with boundary conditions U(A) = UA, U(B) = UB.
%
%    The Gauss-Seidel method is used.
%
%    This routine is provided primarily for comparison with the
%    multigrid solver.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 July 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    William Hager,
%    Applied Numerical Linear Algebra,
%    Prentice-Hall, 1988,
%    ISBN13: 978-0130412942,
%    LC: QA184.H33.
%
%  Parameters:
%
%    Input, integer N, the number of intervals.
%
%    Input, real A, B, the left and right endpoints of the region.
%
%    Input, real UA, UB, the left and right boundary values.
%
%    Input, function value = FORCE ( x ), the name of the function
%    which evaluates the right hand side.
%
%    Input, function value = EXACT ( x ), the name of the function
%    which evaluates the exact solution.
%
%    Output, integer IT_NUM, the number of iterations.
%
%    Output, real U(N+1,1), the computed solution.
%

%
%  Initialization.
%
  tol = 0.0001;
%
%  Set the nodes.
%
  x = ( linspace ( a, b, n + 1 ) )';
%
%  Set the right hand side.
%
  r = zeros ( n + 1, 1 );

  r(1) = ua;
  r(2:n) = force ( x(2:n) ) / n / n;
  r(n+1) = ub;

  u = zeros ( n + 1, 1 );

  it_num = 0;
%
%  Gauss-Seidel iteration.
%
  while ( 1 )

    it_num = it_num + 1;

    [ u, d1 ] = gauss_seidel ( n + 1, r, u );

    if ( d1 <= tol )
      break
    end

  end

  return
end


text:
function lmax = lebesgue_constant ( n, x, xfun )

%*****************************************************************************80
%
%% LEBESGUE_CONSTANT estimates the Lebesgue constant for a set of points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 January 2014
%
%  Author:
%
%    John Burkardt.
%
%  Parameters:
%
%    Jean-Paul Berrut, Lloyd Trefethen,
%    Barycentric Lagrange Interpolation,
%    SIAM Review,
%    Volume 46, Number 3, September 2004, pages 501-517.
%
%  Parameters:
%
%    Input, integer N, the number of interpolation points.
%
%    Input, real X(N), the interpolation points.
%
%    Input, real XFUN(*), the evaluation points.
%
%    Output, real LMAX, an estimate of the Lebesgue constant for the points.
%
  lfun = lebesgue_function ( n, x, xfun );

  lmax = max ( lfun );

  return
end


text:
function value = wedge01_integral ( e )

%*****************************************************************************80
%
%% WEDGE01_INTEGRAL returns the integral of a monomial in the unit wedge in 3D.
%
%  Discussion:
%
%    This routine returns the integral of
%
%      product ( 1 <= I <= 3 ) X(I)^E(I)
%
%    over the unit wedge.
%
%    The integration region is:
%
%      0 <= X
%      0 <= Y
%      X + Y <= 1
%      -1 <= Z <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 August 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971,
%    ISBN: 0130438936,
%    LC: QA311.S85.
%
%  Parameters:
%
%    Input, integer E(3), the exponents.
%
%    Output, real VALUE, the integral of the monomial.
%
  value = 1.0;

  k = e(1);

  for i = 1 : e(2)
    k = k + 1;
    value = value * i / k;
  end

  k = k + 1;
  value = value / k;

  k = k + 1;
  value = value / k;
%
%  Now account for integration in Z.
%
  if ( e(3) == - 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'WEDGE01_INTEGRAL - Fatal error!\n' );
    fprintf ( 1, '  E(3) = -1 is not a legal input.\n' );
    error ( 'WEDGE01_INTEGRAL - Fatal error!' );
  elseif ( mod ( e(3), 2 ) == 1 )
    value = 0.0;
  else
    value = value * 2.0 / ( e(3) + 1 );
  end

  return
end


text:
function q = itimes(q1,q2,takeRight)
% computes inv(o1) .* o2
%
% Syntax
%   q = q1 .* q2
%
%   q = times(q1,q2)
%   q = itimes(q1,q2,1) % inv(o1) .* o2
%   q = itimes(q1,q2,0) % o1 .* inv(o2)
%
% Input
%  q1 - @quaternion
%  q2 - @quaternion
%  takeRight - logical, use as output left or right input
%
% Output
%  q  - @quaternion

% which input will become the output?
if takeRight
  q = q2;
  a1 = q1.a; b1 = -q1.b; c1 = -q1.c; d1 = -q1.d;
  a2 = q2.a; b2 = q2.b; c2 = q2.c; d2 = q2.d;
else
  q = q1;
  a1 = q1.a; b1 = q1.b; c1 = q1.c; d1 = q1.d;
  a2 = q2.a; b2 = -q2.b; c2 = -q2.c; d2 = -q2.d;
end

%standart algorithm
q.a = a1 .* a2 - b1 .* b2 - c1 .* c2 - d1 .* d2;
q.b = b1 .* a2 + a1 .* b2 - d1 .* c2 + c1 .* d2;
q.c = c1 .* a2 + d1 .* b2 + a1 .* c2 - b1 .* d2;
q.d = d1 .* a2 - c1 .* b2 + b1 .* c2 + a1 .* d2;

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
% MODULUS_LAYER Calculate the modulus of coefficients in a layer
%
% Usage
%    U = MODULUS_LAYER(W)
%
% Input
%    W (struct): A scattering layer, as output by WAVELET_LAYER_*, for exam-
%       ple.
%
% Output
%    U (struct): The same layer with all coefficients set to their absolute
%       values.
%
% See Also
%    PAD_SIGNAL, UNPAD_SIGNAL

function U = modulus_layer(W)
	U.signal = cellfun(@abs, W.signal, 'UniformOutput', 0);
	U.meta = W.meta;
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
% predicting numerical values using Linear Regression
clear all;
format long
disp('===== Linear Regression ====');
disp('Reading featur vector');

Indices = crossvalind('Kfold', 1107, 10);
figure;

subplot(3,2,1);

for feat = 1:3
      MSEarray =[];
    elapsedarray =[];
    for crossvalidateIter = 1:10
        (fprintf('%d',crossvalidateIter));
    featurs = csvread('data\productsdesc.csv');
    num_data = size(featurs,1); %5000;
    %disp(sprintf('Number of datapoints %d',num_data))

    possiblefeaturizations =  {'bernouli', 'tfidf','multinomial'};
    %featurization = 'bernouli'%'tfidf'%'tfidf'%'multinomial'%'tfidf' %'multinomial'; % 'bernouli', 'tfidf'
    featurization  = possiblefeaturizations{feat};

    featurs = featurs(:,2:size(featurs,2));
    if strcmp(featurization,'multinomial')
        %just pass
    elseif strcmp(featurization,'bernouli')
        featurs = bernoulli(featurs);
    elseif strcmp(featurization,'tfidf')
        featurs = tfidf(featurs);
    end

     size_training = floor(.9*num_data);

        trainingset = featurs(Indices~=crossvalidateIter,:);
        testset = featurs(Indices==crossvalidateIter,:);

    %disp('Splitting up data into training/test sets');
   [num,txt,raw] = xlsread('data\productsratings.csv');

        responsevals = num;

        responsevals_training = responsevals(Indices~=crossvalidateIter,:);
        responsevals_test = responsevals(Indices==crossvalidateIter,:);

    %disp('Linear Regression');
    %

    tic;

    predictions = [];
    actual = [];
    for i =1:3
        a = responsevals_training(:,i);
        b = responsevals_test(:,i)';
        regresscoeff = regress(a, trainingset);
        C2 = (regresscoeff'*(testset'));
        predictions = [predictions, C2'];
        actual = [actual, b'];
    end

    MSE = mean(sum(((predictions-actual).^2)'))
     elapsed = toc;
        MSEarray = [MSEarray MSE];
        elapsedarray = [elapsedarray elapsed];

    end
    featurization  = possiblefeaturizations{feat}
    fprintf('Avergae MSE for Linear Regression is %0.10f\n', mean(MSEarray));
    subplot(3,2,feat*2-1)
    plot(MSEarray,'-o');
    title(strcat('MSE for Linear Regression (',featurization,') ', sprintf(' avergae MSE = %0.10f\n', mean(MSEarray))));
    xlabel('10 fold cross-validation (iteration no)')
    ylabel('MSE')
    subplot(3,2,feat*2)
    plot(elapsedarray,'-o');
    title(strcat('Elapsed time for Linear Regression (',featurization, ') ' , sprintf(' avergae elapsed time = %0.10f\n', mean(elapsedarray))));
    xlabel('10 fold cross-validation (iteration no)')
    ylabel('Elapsed Time')
    drawnow;
end

author: faridani

repo: MatlabNLP

save_path: github-repos/MATLAB/faridani-MatlabNLP


text:
function[z1,z2,z3] = primtoriem(u1,u2,u3)
	% Transformation from primitive to Riemann invariants
global gamma
u = u2./u1; 			       c2 = gamma*(gamma-1)*(u3./u1 - 0.5*u.^2);
z2 = log((1/gamma)*c2.*u1.^(1-gamma)); c2 = sqrt(c2);
z1 = u - 2*c2/(gamma-1);	       z3 = u + 2*c2/(gamma-1);

repo: Aero-matlab


text:
fprintf('METODO NEWTON-RAPHSON\n'); 			%TITULO
syms f(x) %crea una variable simbolica
f(x)=input('Ingrese la funcion: ');            		% se ingresara la funcion
r0=input('ingrese la primera aproximacion: ');  	%se ingresa el primer punto o raiz
tol=input('Ingrese la tolerancia: ');           	%valor del error minimo aceptable
n=input('Ingrese numero maximo de iteraciones: '); 	%hasta cuantas iteraciones
df=diff(f,x);       %derivada de f(x)
error=100;          %para ingresar al while
i=0;                 %contador
ri=r0;
fprintf('iteraciones\t\t\t\tri\t\t\t\t\tri+1\t\t\t\terror\n');    	%el encabezado a imprimir
fprintf('%i\t\t\t\t%4.11f\t\t\t\t---------\t\t\t\t-----------\n',i,r0); %imprime nuestros datos inciales
while(error>=tol && i<=n)
   i=i+1;
   rim1=ri-(f(ri)/df(ri)); %hallamos la nueva raiz
   error=abs(((ri-rim1)/rim1)*100);     				%hallamos el error al tratar de encontrar la raiz
   fprintf('%i\t\t\t\t%4.11f\t\t\t\t%4.11f\t\t\t\t%4.11f\n',i,r0,rim1,error); %se imprimen los valores hallados por el metodo newton
   ri=rim1; 								%en caso no termine el bucle se vuelve a repetir con una nueva raiz
end
fprintf('La raiz es %4.11f y el error es %4.11f\n',rim1,error);     	%se imprimen los datos finales

repo: ARTE


function [A_lon,B_lon,A_lat,B_lat] = compute_ss_model(filename,x_trim,u_trim)
% x_trim is the trimmed state,
% u_trim is the trimmed input

% add stuff here
[A,B,C,D] = linmod(filename, x_trim, u_trim);

E1 = [...
      0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0;...
      0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0;...
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1;...
      0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0;...
      0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0;...
      ];
E2 = [...
      0, 1, 0, 0;...
      0, 0, 1, 0;...
      ];
A_lat = E1 * A * E1';
B_lat = E1 * B * E2';

E3 = [...
      0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0;...
      0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0;...
      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0;...
      0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0;...
      0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0;...
      ];
E4 = [...
      1, 0, 0, 0;...
      0, 0, 0, 1;...
      ];

A_lon = E3 * A * E3';
B_lon = E3 * B * E4';


text:
% This function output the support polygon depending on the current contact
% status

function CoM_limits = SupportPolygonWRTContactStatusFCN(contactStatus,Config)
CoM_limits = Config.CoM_limit_double;

if(and(contactStatus(1),contactStatus(2))) % Double Support
    CoM_limits = Config.CoM_limit_double;
elseif(and(~contactStatus(1),~contactStatus(2))) % Flying phase
    CoM_limits = Config.CoM_limit_double;
elseif(contactStatus(1)) % Left Foot contact
    CoM_limits =Config.CoM_limit_left;
elseif(contactStatus(2)) % Right Foot contact
    CoM_limits = Config.CoM_limit_right;
end
end

author: robotology

repo: whole-body-controllers


text:
function [ accuracies, F1s, corrs, ccc, rms, classes ] = evaluate_au_prediction_results( labels, ground_truth )
%EVALUATE_CLASSIFICATION_RESULTS Summary of this function goes here
%   Detailed explanation goes here

    classes = sort(unique(ground_truth));
    accuracies = zeros(numel(classes),1);
    F1s = zeros(numel(classes),1);

    corrs = corr(labels, ground_truth);

    rms = sqrt(mean((labels-ground_truth).^2));

    std_g = std(ground_truth);
    std_p = std(labels);

    ccc = 2 * corrs * std_g * std_p / (std_g^2 + std_p^2 + (mean(labels) - mean(ground_truth))^2);

    % the label is taken to belong to a class it is closest to
    label_dists = zeros(numel(labels), numel(classes));

    for i=1:numel(classes)
        label_dists(:,i) = abs(labels - classes(i));
    end

    [~, labels] = min(label_dists');
    labels = labels';

    for i=1:numel(classes)
       labels(labels==i) = classes(i);
    end

    for i=1:numel(classes)

        pos_samples = ground_truth == classes(i);
        neg_samples = ground_truth ~= classes(i);

        pos_labels = labels == classes(i);
        neg_labels = labels ~= classes(i);

        TPs = sum(pos_samples & pos_labels);
        TNs = sum(neg_samples & neg_labels);

        FPs = sum(pos_labels & neg_samples);
        FNs = sum(neg_labels & pos_samples);

        accuracies(i) = (TPs + TNs) / numel(pos_samples);

        F1s(i) = 2 * TPs / (2*TPs + FNs + FPs);

    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function Y=GeneLableY(TrainLabel,Num)

N=length(TrainLabel);
Y=zeros(Num,N);

for i=1:N
    Y(TrainLabel(:,i),i)=1;
end
end

author: BehnoodRasti


text:
if (~isdeployed)
    addpath('lib/pose');
    addpath('lib/pose/multicut');
    addpath('lib/utils');
    addpath('lib/vis');
    addpath('lib/eval');
    addpath('lib/multicut');
    addpath('lib/multicut/hdf5');
    caffe_dir = 'external/caffe/matlab/';
    if exist(caffe_dir)
        addpath(caffe_dir);
    else
        warning('Please install Caffe in ./external/caffe');
    end
    addpath('external/liblinear-1.94/matlab/')
    fprintf('Pose startup done\n');
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function f = or(f, g)
%|   CHEBTECH logical OR.
%   F | G performs a logical OR of the CHEBTECH objects F and G and returns a
%   CHEBTECH containing elements set to either logical 1 (TRUE) or logical 0
%   (FALSE).  An element of the output CHEBTECH is set to 1 if either input
%   CHEBTECH contains a non-zero element at that point, otherwise it is set to
%   0.  F and G must either be identically zero or have roots in their domains.
%   If this is not the case, garbage is returned with no warning.

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

arbitraryPoint = 0.1273881594;
f.coeffs = feval(f, arbitraryPoint) | feval(g, arbitraryPoint);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function Training_Data = ReadFace(Training_Path)
% ---------- Construct 2D matrix from all of the 1D image vectors in the training data file ------------
flist = dir(strcat(Training_Path,'\*.jpg'));
Training_Data = [];
for imidx = 1:length(flist)
    fprintf('Constructing Training Image Data Space [%d] \n', imidx);
    path = strcat(Training_Path,strcat('\',int2str(imidx),'.jpg'));
    img = imread(path);
    [irow icol] = size(img);
    temp = reshape(img',irow*icol,1);   % Reshaping 2D images to 1D image vectors
    Training_Data = [Training_Data temp];
end
fprintf('\n');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


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

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

% Last Modified by GUIDE v2.5 26-Aug-2012 15:11:43

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

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

% --- Executes just before SCARA_GUI is made visible.
function SCARA_GUI_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to SCARA_GUI (see VARARGIN)

% Choose default command line output for SCARA_GUI
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes SCARA_GUI wait for user response (see UIRESUME)
% uiwait(handles.figure1);
SCARA28 = vrworld('SCARA28.wrl');
open(SCARA28)
view(SCARA28);

dist=-get(handles.slider5,'value');
y1=-.5*dist;
y2=.5*dist;
SCARA28.EndV1.translation = [-.35, y1, 0];
SCARA28.EndV2.translation = [-.35, y2, 0];

d4=-get(handles.slider4,'value');
SCARA28.d4b.translation=[0,d4,0];

T4=get(handles.slider3,'value');
SCARA28.d4b.rotation = [0, 1, 0,T4*pi/180];

T2=-get(handles.slider2,'value');
SCARA28.a2.rotation = [1, 0, 0, T2*pi/180];

T1=get(handles.slider1,'value');
SCARA28.d1.rotation = [0, 1, 0, T1*pi/180];
clc
% --- Outputs from this function are returned to the command line.
function varargout = SCARA_GUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

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

% --- Executes on slider movement.
function slider1_Callback(hObject, eventdata, handles)
% hObject    handle to slider1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of
%        slider

SCARA28 = vrworld('SCARA28.wrl');
open(SCARA28)
T1=get(handles.slider1,'value');
if T1==360
   set(handles.slider1,'value',0);
end
SCARA28.d1.rotation = [0, 1, 0, T1*pi/180];

% --- Executes during object creation, after setting all properties.
function slider1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function slider2_Callback(hObject, eventdata, handles)
% hObject    handle to slider2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
SCARA28 = vrworld('SCARA28.wrl');
open(SCARA28)
T2=-get(handles.slider2,'value');
SCARA28.a2.rotation = [1, 0, 0, T2*pi/180];

% --- Executes during object creation, after setting all properties.
function slider2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function slider3_Callback(hObject, eventdata, handles)
% hObject    handle to slider3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
SCARA28 = vrworld('SCARA28.wrl');
open(SCARA28)
T4=get(handles.slider3,'value');
if T4==360
   set(handles.slider3,'value',0);
end
SCARA28.d4b.rotation = [0, 1, 0,T4*pi/180];

% --- Executes during object creation, after setting all properties.
function slider3_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function slider4_Callback(hObject, eventdata, handles)
% hObject    handle to slider4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
SCARA28 = vrworld('SCARA28.wrl');
open(SCARA28)
d4=-get(handles.slider4,'value');
SCARA28.d4b.translation=[0,d4,0];

% --- Executes during object creation, after setting all properties.
function slider4_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --- Executes on slider movement.
function slider5_Callback(hObject, eventdata, handles)
% hObject    handle to slider5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
SCARA28 = vrworld('SCARA28.wrl');
open(SCARA28)
dist=-get(handles.slider5,'value');
y1=-.5*dist;
y2=.5*dist;
SCARA28.EndV1.translation = [-.35, y1, 0];
SCARA28.EndV2.translation = [-.35, y2, 0];

% --- Executes during object creation, after setting all properties.
function slider5_CreateFcn(hObject, eventdata, handles)
% hObject    handle to slider5 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background.
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

% --------------------------------------------------------------------
function File_tag_Callback(hObject, eventdata, handles)
% hObject    handle to File_tag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function Edit_tag_Callback(hObject, eventdata, handles)
% hObject    handle to Edit_tag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------
function Reset_tag_Callback(hObject, eventdata, handles)
% hObject    handle to Reset_tag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
SCARA28 = vrworld('SCARA28.wrl');
open(SCARA28)
dist=-1.3;
y1=-.5*dist;
y2=.5*dist;
SCARA28.EndV1.translation = [-.35, y1, 0];
SCARA28.EndV2.translation = [-.35, y2, 0];
d4=0;
SCARA28.d4b.translation=[0,d4,0];
T4=0;
SCARA28.d4b.rotation = [0, 1, 0,T4*pi/180];
T2=0;
SCARA28.a2.rotation = [1, 0, 0, T2*pi/180];
T1=0;
SCARA28.d1.rotation = [0, 1, 0, T1*pi/180];
clc

set(handles.slider5,'value',-dist);
set(handles.slider4,'value',0);
set(handles.slider3,'value',0);
set(handles.slider2,'value',0);
set(handles.slider1,'value',0);

% --------------------------------------------------------------------
function About_us_tag_Callback(hObject, eventdata, handles)
% hObject    handle to About_us_tag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

about_us

% --- Executes during object creation, after setting all properties.
function figure1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% --------------------------------------------------------------------
function Exit_Tag_Callback(hObject, eventdata, handles)
% hObject    handle to Exit_Tag (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

author: Sable

repo: mcbench-benchmarks


text:
function draw_points (pnts)
%DRAW_POINTS
% draw points one by one in a new figure
% usages:
%   draw_points(pnts);
% pnts must be a Mx2 double matrix

M = size(pnts, 2);
if (M~=2 || ~isa(pnts,'double'))
    error('pnts must be a Mx2 double matrix');
end
hfig = figure;
hax = axes('parent', hfig);
plot(pnts(:,1),pnts(:,2),...
    'parent', hax,...
    'linestyle', 'none',...
    'marker', 'o');
axis('square');
set(hfig, 'name', 'example: view user-defined data type');


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
    = m0_start(fid, data, varargins, opts)
% This file is an entry for the M0 strategy.
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
%           = m0_start(fid, data, varargins, opts)
% cum_ret: cumulative wealth achived at the end of a period.
% cumprod_ret: cumulative wealth achieved till the end each period.
% daily_ret: daily return achieved by a strategy.
% daily_portfolio: daily portfolios
%
% data: market sequence vectors
% fid: handle for write log file
% varargins: variable parameters
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
%          = m0_start(fid, data, {0.5, 0}, opts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract the parameters
beta = varargins{1};    %
tc = varargins{2};      % transaction cost fee rate

% Run the M0 algorithm
[cum_ret, cumprod_ret, daily_ret, daily_portfolio] ...
    = m0_run(fid, data, beta, tc, opts);

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function r = drihaczekdist(f)
%DRIHACZEKDIST discrete Rihaczek distribution
%   Usage r = drihaczekdist(f);
%
%
%   `drihaczekdist(f)` computes a discrete Rihaczek distribution of vector
%   *f*. The discrete Rihaczek distribution is computed by
%
%   .. math:: r\left( k+1,\; l+1 \right)\; =\; f\left( l+1 \right)\; \overline{c\left( k+1 \right)}e^{-2\pi ikl/L}
%
%   where $k, l=0,\ldots,L-1$ and $c$ is the Fourier transform of $f$.
%
%   **WARNING**: The quadratic time-frequency distributions are highly
%   redundant. For an input vector of length L, the quadratic time-frequency
%   distribution will be a $L \times L$ matrix. If *f* is multichannel
%   ($L\times W$ matrix), the resulting distributions are stacked along
%   the third dimension such that the result is $L\times L \times W$ cube.

%   AUTHOR: Jordy van Velthoven
%   TESTING: TEST_DRIHACZEKDIST
%   REFERENCE: REF_DRIHACZEKDIST

complainif_notenoughargs(nargin, 1, 'DRIHACZEKDIST');

[f,Ls]=comp_sigreshape_pre(f,upper(mfilename));

c = dgt(f, f, 1, Ls);

r = dsft(c);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function S = SeparateScales(X,L)
% SeparatesScales: Separates an image into an (orthonormal)
%                                    series of disjoint scales.
%  Usage:
%    S =  SeparateScales(X,L);
%  Inputs:
%    X     n by n image
%    L     scale of the coarsest coefficients
%  Outputs: S data structure which contains the FT of X tapered at
%           dyadic scales.
%
% Description
%   For each s = 1, ..., J - L,
%   S{s} is  a structure which contains the DFT of X tapered with a
%   Meyer-style window which isolates frequencies near a dyadic
%   frequency subband
%
%   S{s}         m by m matrix with m = 2*2^(j+1)); matrix
%                entries are Fourier coefficients of Img at
%                -2^(j+1) <= k1,k2 < 2^(j+1) after windowing with a
%                bandpass Meyer window which isolates frequencies
%                in the range 2^(j-1) <= |k| <= 2^j.
%
%                j = L - 2 + s
% See Also
%   Adj_ Separatescales, Curvelet02Xform
%
% By Emmanuel Candes, 2003-2004

  n = size(X,1);
  n2 = n/2;
  J = log2(n);
  deg = 3;

  scale = (L-1):(J-1);
  nscales = length(scale);
  S = cell(1,nscales);  % Create data structure

  F = fft2_mid0(X)/sqrt(prod(size(X)));

  %      Partition at Coarse Level
  [ix,w] = CoarseMeyerWindow(L-1,deg);
  w = [0 reverse(w(2:length(w))) w];

  w2 = w'*w; % Build 2D-window

  l = 2^L;
  lx = (n2 - l + 1):(n2 + l);
  S{1} = F(lx,lx).*w2;

  %      Loop to Get Partition for  j = L, ..., J - 1;
  for j = L:(J-1),

    l = min(2^(j+1),n2);
    wlo = zeros(1,l); whi = wlo;

    [ixf, wf] = CoarseMeyerWindow(j-1,deg);
    wlo(ixf+1)  = wf;
    wlo = [0 reverse(wlo(2:length(wlo))) wlo];

    if j < J -1,
      [ixp, wp] = CoarseMeyerWindow(j,deg);
      whi(ixp+1) = wp;
      whi = [0 reverse(whi(2:length(whi))) whi];
    end

    if j == J - 1,
      whi = ones(1,2*l);
    end

    w2 = sqrt((whi'*whi).^2 - (wlo'*wlo).^2); % Build 2D-window

    lx = (n2 - l + 1):(n2 + l);
    S{j - L + 2} = F(lx,lx).*w2;
  end

author: thomaskuestner


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function [eu, F] = sp_eval_stress_msh (u, space, msh, lambda_lame, mu_lame)

warning ('geopdes:obsolete', 'The function SP_EVAL_STRESS_MSH is obsolete. Using SP_EVAL_MSH instead. Read the help for the usage')

[eu, F] = sp_eval_msh (u, space, msh, 'stress', lambda_lame, mu_lame);

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function MRSIStruct = setFov(MRSIStruct, dimension, value)
    switch(lower(dimension))
        case('x')
            MRSIStruct.fov.x = value;
        case('y')
            MRSIStruct.fov.y = value;
        case('z')
            MRSIStruct.fov.z = value;
    end
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function obj = updateRectangle(obj, rectIndex)

%----RECTANGLE FIELDS----%

% x - [DONE]
% y - [DONE]
% mode - [DONE]
% name - [DONE]
% text - [NOT SUPPORTED IN MATLAB]
% error_y - [HANDLED BY ERRORBAR]
% error_x - [HANDLED BY ERRORBAR]
% line.color - [DONE]
% line.width - [DONE]
% line.dash - [DONE]
% line.opacity - [NOT SUPPORTED IN MATLAB]
% line.smoothing - [NOT SUPPORTED IN MATLAB]
% line.shape - [NOT SUPPORTED IN MATLAB]
% connectgaps - [NOT SUPPORTED IN MATLAB]
% fill - [HANDLED BY RECTANGLE]
% fillcolor - [HANDLED BY RECTANGLE]
% opacity - [NOT SUPPORTED IN MATLAB]
% textfont - [NOT SUPPORTED IN MATLAB]
% textposition - [NOT SUPPORTED IN MATLAB]
% xaxis [DONE]
% yaxis [DONE]
% showlegend [DONE]
% stream - [HANDLED BY PLOTLYSTREAM]
% visible [DONE]
% type [DONE]

%-AXIS INDEX-%
axIndex = obj.getAxisIndex(obj.State.Plot(rectIndex).AssociatedAxis);

%-RECTANGLE DATA STRUCTURE- %
rect_data = get(obj.State.Plot(rectIndex).Handle);

%-CHECK FOR MULTIPLE AXES-%
[xsource, ysource] = findSourceAxis(obj,axIndex);

%-AXIS DATA-%
eval(['xaxis = obj.layout.xaxis' num2str(xsource) ';']);
eval(['yaxis = obj.layout.yaxis' num2str(ysource) ';']);

%-------------------------------------------------------------------------%

%-rectangle xaxis-%
obj.data{rectIndex}.xaxis = ['x' num2str(xsource)];

%-------------------------------------------------------------------------%

%-rectangle yaxis-%
obj.data{rectIndex}.yaxis = ['y' num2str(ysource)];

%-------------------------------------------------------------------------%

%-rectangle type-%
obj.data{rectIndex}.type = 'scatter';

%-------------------------------------------------------------------------%

%-rectangle x-%
obj.data{rectIndex}.x = [rect_data.Position(1) rect_data.Position(1) ...
    rect_data.Position(1) + rect_data.Position(3) ...
    rect_data.Position(1) + rect_data.Position(3) ...
    rect_data.Position(1)];

%-------------------------------------------------------------------------%

%-rectangle y-%
obj.data{rectIndex}.y = [rect_data.Position(2) rect_data.Position(2) + rect_data.Position(4) ...
    rect_data.Position(2) + rect_data.Position(4) ...
    rect_data.Position(2) ...
    rect_data.Position(2)];

%-------------------------------------------------------------------------%

%-rectangle name-%
obj.data{rectIndex}.name = rect_data.DisplayName;

%-------------------------------------------------------------------------%

%-rectangle mode-%
obj.data{rectIndex}.mode = 'lines';

%-------------------------------------------------------------------------%

%-rectangle fill-%
obj.data{rectIndex}.visible = strcmp(rect_data.Visible,'on');

%-------------------------------------------------------------------------%

%-rectangle fill-%
obj.data{rectIndex}.fill = 'tonexty';

%-------------------------------------------------------------------------%

%-rectangle line-%
obj.data{rectIndex}.line = extractPatchLine(rect_data);

%-------------------------------------------------------------------------%

%-rectangle fillcolor-%
fill = extractPatchFace(rect_data);
obj.data{rectIndex}.fillcolor = fill.color;

%-------------------------------------------------------------------------%

%-rectangle showlegend-%
leg = get(rect_data.Annotation);
legInfo = get(leg.LegendInformation);

switch legInfo.IconDisplayStyle
    case 'on'
        showleg = true;
    case 'off'
        showleg = false;
end

obj.data{rectIndex}.showlegend = showleg;

%-------------------------------------------------------------------------%

end

author: plotly


text:
function x = simuARMA(a,b,nobs)
% x = simuARMA(a,b,nobs)
% simulation of nobs observations of an ARMA process with parameters
% a,b en standard deviation 1.
%
% No transients in finite signal
% Paper:
%
%   Broersen, P.M.T. and S. de Waele
% 	Generating Data with Prescribed Power Spectral Density.
%	IEEE Trans. on Instrumentation and Measurement, vol. 52, no. 4, p 1061-1067, 2003.
%
%
% uses randn( , )
%

aro = length(a)-1;
mao = length(b)-1;
[cor,g]=arma2cor(a,b);
[cor2,gar]=arma2cor(a,1);

if ~aro,
   v = randn(nobs+mao,1)/sqrt(g);
else
   z = randn(aro,1);
   vs = zeros(aro,1);
   [at rc] = ar2arset(a);
   f = 1;
   vs(1) = f*z(1);
   al = [];
   for i = 2:aro,
   	f = sqrt(1-rc(i)^2)*f;
      al = [al 0] + rc(i)*[fliplr(al) 1];
      vs(i) = f*z(i)-al*vs(i-1:-1:1);
   end %for i = 2:aro,
   vs = vs*sqrt(gar/g);
   v = armafilter(randn(nobs+mao,1)/sqrt(g),a,1,vs,z);
end %if ~aro,
x = armafilter(v,1,b);
x = x(mao+1:mao+nobs);


text:
function  jitterLags = realized_kernel_jitter_lag_length(noiseEstimate,iqEstimate,kernel,N)
% Computes the optimal amount of end point jitter given a kernel, noise estimate, integrated
% quarticity estimate and number of observations.
%
% USAGE:
%   [JITTERLAGS] = realized_kernel_weights(NOISEESTIMATE,IQESTIMATE,KERNEL,N)
%
% INPUTS:
%   NOISEESTIMATE - Estimated variance of the noise present in a high freuqency asset price (see
%                     realized_kernel_select_lag_length)
%   IQESTIMATE    - Estimated integrated quarticity.  For the purposes of estimating the number of
%                     points to jitter it is often reasonable to use the square of a low frequence
%                     IV measure.
%   KERNEL        - String containing one of the supported kernel types:
%                     Non-flat-top, weakly positive, n^(1/5) rate:
%                     - 'nonflatparzen' [RECOMMENDED] Parzen kernel applied to non-flat-top case
%                     - 'qs' Quadratic Spectral
%                     - 'fejer' Fejer kernel
%                     - 'thinf' Tukey-Hanning kernel with infinite lag order
%                     - 'bnhls' Kernel proposed by Barndorf-Neilsen, Hansen, Lunde and Shephard
%                     Flat-top, n^(1/4) rate:
%                     - 'parzen' Parzen's kernel
%                     - 'th1' Tukey Hanning kernel with power 1
%                     - 'th2' Tukey Hanning kernel with power 2
%                     - 'th5' Tukey Hanning kernel with power 5
%                     - 'th16' Tukey Hanning kernel with power 16
%                     - 'cubic','multiscale' 3rd order (cubic) kernel, which
%                         is asymptotically equivalent to the multiscale estimator
%                     - '5thorder' 5th order Kernel
%                     - '6thorder' 6th order Kernel
%                     - '7thorder' 7th order Kernel
%                     - '8thorder' 8th order Kernel
%                     Flat-top, n^(1/6) rate:
%                     - 'bartlett','twoscale' Bartlett kernel, which is
%                         asymptotically equivalent to the two-scale estimator
%                     - '2ndorder' 2nd order (quadratic) kernel
%                     - 'epanechnikov' Epanechnikov kernel
%   N                 - Number of observation of the price
%
% OUTPUTS:
%   JITTERLAGS    - Estimated optimal number of lags to jitter the end points of a kernel.  If 1
%                     then no pre-averaging is needindicated
%
% COMMENTS:
%   This is a helper function for REALIZED_KERNEL. See Barndorf-Nielsen,
%   Hansen, Lunde and Shephard (2008a, 2008b) for details about the kernel and
%   their properties.
%
%  See also REALIZED_KERNEL_SELECT_LAG_LENGTH, REALIZED_VARIANCE, REALIZED_KERNEL,
%  REALIZED_PRICE_FILTER, REALIZED_KERNEL_WEIGHTS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin~=4
    error('Four inputs required.')
end
if length(noiseEstimate)>1 || noiseEstimate<0
    error('NOISEESTIMATE must be a non-negative scalar.');
end
if length(iqEstimate)>1 || iqEstimate<0
    error('IQESTIMATE must be a non-negative scalar.');
end
if ~ismember(kernel,{'bartlett','twoscale','2ndorder','epanechnikov','cubic','multiscale','5thorder','6thorder','7thorder','8thorder','parzen','th1','th2','th5','th16','nonflatparzen','qs','fejer','thinf','bnhls'})
    error('KERNEL must be one of the listed types.')
end

if ~isscalar(N) || N<1 || floor(N)~=N
    error('N must be an integer greater than 1.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Look up table
switch lower(kernel)
    case {'bartlett','twoscale'}
        cStar = 2.28;
        k00 = 1/3;
        kernelType = 2;
    case '2ndorder'
        cStar = 3.42;
        k00 = 1/5;
        kernelType = 2;
    case 'epanechnikov'
        cStar = 2.46;
        k00 = 8/15;
        kernelType = 2;
    case {'cubic','multiscale'}
        cStar = 3.68;
        k00 = 0.371;
        k11 = 1.20;
        k22 = 12.0;
        kernelType = 1;
    case '5thorder'
        k00 = 0.391;
        k11 = 1.42;
        k22 = 17.1;
        cStar = 3.70;
        kernelType = 1;
    case '6thorder'
        k00 = 0.471;
        k11 = 1.55;
        k22 = 22.8;
        cStar = 3.97;
        kernelType = 1;
    case '7thorder'
        k00 = 0.533;
        k11 = 1.71;
        k22 = 31.8;
        cStar = 4.11;
        kernelType = 1;
    case '8thorder'
        k00 = 0.582;
        k11 = 1.87;
        k22 = 43.8;
        cStar = 4.31;
        kernelType = 1;
    case 'parzen'
        k00 = 0.269;
        k11 = 1.50;
        k22 = 24.0;
        cStar = 4.77;
        kernelType = 1;
    case 'th1'
        k00 = 0.375;
        k11 = 1.23;
        k22 = 12.1;
        cStar = 3.70;
        kernelType = 1;
    case 'th2'
        k00 = 0.219;
        k11 = 1.71;
        k22 = 41.7;
        cStar = 5.74;
        kernelType = 1;
    case 'th5'
        k00 = 0.097;
        k11 = 3.50;
        k22 = 489.0;
        cStar = 8.07;
        kernelType = 1;
    case 'th16'
        k00 = 0.032;
        k11 = 10.26;
        k22 = 14374.0;
        cStar = 39.16;
        kernelType = 1;
    case 'nonflatparzen'
        cStar = ((12)^2/0.269)^(1/5);
        k00 = .269;
        kernelType = 3;
    case 'qs'
        cStar = ((1/5)^2/(3*pi/5))^(1/5);
        k00 = 3*pi/5;
        kernelType = 3;
    case 'fejer'
        cStar = ((2/3)^2/(pi/3))^(1/5);
        k00 = pi/3;
        kernelType = 3;
    case 'thinf'
        cStar = ((pi^2/2)^2/(.52))^(1/5);
        k00 = 0.52;
        kernelType = 3;
    case 'bnhls'
        cStar = (1^2/(5/4))^(1/5);
        k00 = 5/4;
        kernelType = 3;
    otherwise
        error('KERNEL must be one of the listed types.')
end

% Compute the quantities necessary
if kernelType ==1
    % n^1/4 flat top
    % Approximate at the constant variance solution
    d=k00*k22/k11^2;
    fd = sqrt(1+sqrt(1+3*d));
    g = sqrt(k00*k11)*(1/fd + fd);
    avar = 16/3 * g * noiseEstimate * iqEstimate^(3/4);
    power = 1/4;
elseif kernelType == 2
    % n^1/6 flat top
    avar = 6 * cStar * k00 * noiseEstimate^(4/3) * iqEstimate^(2/3);
    power = 1/6;
elseif kernelType == 3
    % n^1/5 non flat top
    avar = 5 * cStar * k00 * noiseEstimate^(4/5) * iqEstimate^(4/5);
    power = 1/5;
end

% Problem is to min 8*noise^2*jitterLag^(-2) + avar * (N-jitterLag)^(-2*power)
% jittering over m data points, N observations in total
MSE = inf*ones(N-1,1);
for m=1:(N-1)
    MSE(m) = 8 * noiseEstimate^2 * m^(-2) + avar * (N-m)^(-2*power);
end

% Find the minimum, which may be 1
[temp,jitterLags] = min(MSE);

author: bashtage

repo: mfe-toolbox


text:
function [Signal, Sigma] = MPdenoising(data, mask, kernel, sampling, centering)
    %
    % "MPPCA": 4d image denoising and noise map estimation by exploiting  data redundancy in the PCA domain using universal properties of the eigenspectrum of
    % random covariance matrices, i.e. Marchenko Pastur distribution
    %
    %  [Signal, Sigma] = MPdenoising(data, mask, kernel, sampling)
    %       output:
    %           - Signal: [x, y, z, M] denoised data matrix
    %           - Sigma: [x, y, z] noise map
    %       input:
    %           - data: [x, y, z, M] data matrix
    %           - mask:   (optional)  region-of-interest [boolean]
    %           - kernel: (optional)  window size, typically in order of [5 x 5 x 5]
    %           - sampling:
    %                    1. full: sliding window (default for noise map estimation, i.e. [Signal, Sigma] = MPdenoising(...) )
    %                    2. fast: block processing (default for denoising, i.e. [Signal] = MPdenoising(...))
    %
    %  Authors: Jelle Veraart (jelle.veraart@nyumc.org)
    % Copyright (c) 2016 New York Universit and University of Antwerp
    %
    %      Permission is hereby granted, free of charge, to any non-commercial entity
    %      ('Recipient') obtaining a copy of this software and associated
    %      documentation files (the 'Software'), to the Software solely for
    %      non-commercial research, including the rights to use, copy and modify the
    %      Software, subject to the following conditions:
    %
    %        1. The above copyright notice and this permission notice shall be
    %      included by Recipient in all copies or substantial portions of the
    %      Software.
    %
    %        2. THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND,
    %      EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIESOF
    %      MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
    %      NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BELIABLE FOR ANY CLAIM,
    %      DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
    %      OTHERWISE, ARISING FROM, OUT OF ORIN CONNECTION WITH THE SOFTWARE OR THE
    %      USE OR OTHER DEALINGS IN THE SOFTWARE.
    %
    %        3. In no event shall NYU be liable for direct, indirect, special,
    %      incidental or consequential damages in connection with the Software.
    %      Recipient will defend, indemnify and hold NYU harmless from any claims or
    %      liability resulting from the use of the Software by recipient.
    %
    %        4. Neither anything contained herein nor the delivery of the Software to
    %      recipient shall be deemed to grant the Recipient any right or licenses
    %      under any patents or patent application owned by NYU.
    %
    %        5. The Software may only be used for non-commercial research and may not
    %      be used for clinical care.
    %
    %        6. Any publication by Recipient of research involving the Software shall
    %      cite the references listed below.
    %
    % REFERENCES
    %      Veraart, J.; Fieremans, E. & Novikov, D.S. Diffusion MRI noise mapping
    %      using random matrix theory Magn. Res. Med., 2016, early view, doi:
    %      10.1002/mrm.26059

    if isa(data,'integer')
        data = single(data);
    end
    [sx, sy, sz, M] = size(data);

    if ~exist('mask', 'var') || isempty(mask)
        mask = true([sx, sy, sz]);
    end
    if ~isa(mask,'boolean')
        mask = mask>0;
    end

    if ~exist('kernel', 'var') || isempty(kernel)
        kernel = [5 5 5];
    end

    if isscalar(kernel)
        kernel = [kernel, kernel, kernel];
    end
    kernel = kernel + (mod(kernel, 2)-1);   % needs to be odd.
    k = (kernel-1)/2; kx = k(1); ky = k(2); kz = k(3);
    N = prod(kernel);

    if ~exist('sampling', 'var') || isempty(sampling)
        if nargout > 1
            sampling = 'full';
        else
            sampling = 'fast';
        end
    end

    % create mask
    if ~exist('mask', 'var') || isempty(mask)
        mask = true(sx, sy, sz);
    end

    if ~exist('centering', 'var') || isempty(centering)
        centering = false;
    end

    if strcmp(sampling, 'fast')
        if nargout>1
            warning('undersampled noise map will be returned')
        end
        % compute center points of patches
        stats = regionprops(mask, 'BoundingBox');
        n = ceil(stats.BoundingBox(4:6) ./ kernel);

        x = linspace(ceil(stats.BoundingBox(1))+k(1), floor(stats.BoundingBox(1))-k(1) + stats.BoundingBox(4), n(1)); x = round(x);
        y = linspace(ceil(stats.BoundingBox(2))+k(2), floor(stats.BoundingBox(2))-k(2) + stats.BoundingBox(5), n(2)); y = round(y);
        z = linspace(ceil(stats.BoundingBox(3))+k(3), floor(stats.BoundingBox(3))-k(3) + stats.BoundingBox(6), n(3)); z = round(z);

        [y, x, z] = meshgrid(x, y, z); x = x(:); y = y(:); z = z(:);
    end

    if strcmp(sampling, 'full')
        warning('image boundaries are not processed.')
        mask(1:k(1), :, :) = 0;
        mask(sx-k(1):sx, :, :) = 0;

        mask(:, 1:k(2), :) = 0;
        mask(:, sy-k(2):sy, :, :) = 0;
        mask(:,:,1:k(3)) = 0;
        mask(:,:,sz-k(3)) = 0;

        x = []; y = []; z = [];
        for i = k(3)+1:sz-k(3)
            [x_, y_] = find(mask(:,:,i) == 1);
            x = [x; x_]; y = [y; y_];  z = [z; i*ones(size(y_))];
        end
        x = x(:); y = y(:); z = z(:);
    end

    % Declare variables:
    if logical(exist('OCTAVE_VERSION', 'builtin')) % for Octave, zeros(..,'like',data) not implemented. replace with class(data).
        zerosoption = {class(data)};
    else
        zerosoption = {'like', data};
    end
    sigma = zeros(1, numel(x), zerosoption{:});
    npars = zeros(1, numel(x), zerosoption{:});
    signal = zeros(M, prod(kernel), numel(x), zerosoption{:});

    Sigma = zeros(sx, sy, sz, zerosoption{:});
    Npars = zeros(sx, sy, sz, zerosoption{:});
    Signal = zeros(sx, sy, sz, M, zerosoption{:});

    % compute scaling factor for in case N<M
    R = min(M, N);
    scaling = (max(M, N) - (0:R-centering-1)) / N;
    scaling = scaling(:);

    % start denoising
    for nn = 1:numel(x)

        % create data matrix
        X = data(x(nn)-kx:x(nn)+kx, y(nn)-ky:y(nn)+ky, z(nn)-kz:z(nn)+kz, :);
        X = reshape(X, N, M); X = X';

        if centering
            colmean = mean(X, 1);
            X = X - repmat(colmean, [M, 1]);
        end
        % compute PCA eigenvalues
        [u, vals, v] = svd(X, 'econ');
        vals = diag(vals).^2 / N;

        % First estimation of Sigma^2;  Eq 1 from ISMRM presentation
        csum = cumsum(vals(R-centering:-1:1)); cmean = csum(R-centering:-1:1)./(R-centering:-1:1)'; sigmasq_1 = cmean./scaling;

        % Second estimation of Sigma^2; Eq 2 from ISMRM presentation
        gamma = (M - (0:R-centering-1)) / N;
        rangeMP = 4*sqrt(gamma(:));
        rangeData = vals(1:R-centering) - vals(R-centering);
        sigmasq_2 = rangeData./rangeMP;

        % sigmasq_2 > sigma_sq1 if signal-components are represented in the
        % eigenvalues

        t = find(sigmasq_2 < sigmasq_1, 1);

        if isempty(t)
            sigma(nn) = NaN;
            signal(:, :, nn) = X;
            t = R+1;
        else
            sigma(nn) = sqrt(sigmasq_1(t));
            vals(t:R) = 0;
            s = u*diag(sqrt(N*vals))*v';
            if centering
               s = s + repmat(colmean, [M, 1]);
            end

            signal(:, :, nn) = s;
        end
        npars(nn) = t-1;
    end

    for nn = 1:numel(x)
        Sigma(x(nn), y(nn), z(nn)) = sigma(nn);
        Npars(x(nn), y(nn), z(nn)) = npars(nn);
        if strcmp(sampling, 'fast')
            Signal(x(nn)-k(1):x(nn)+k(1),y(nn)-k(2):y(nn)+k(2),z(nn)-k(3):z(nn)+k(3), :) = unpatch(signal(:,:,nn), k);
        elseif strcmp(sampling, 'full')
            Signal(x(nn), y(nn),z(nn), :) = signal(:,ceil(prod(kernel)/ 2),nn);
        end
    end
end

function data = unpatch(X, k)
    kernel=k+k+1;
    data = zeros([kernel, size(X, 1)]);
    tmp = zeros(kernel);
    for i = 1:size(X, 1);
        tmp(:) = X(i, :);
        data(:,:,:,i) = tmp;
    end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function calpak_test0005 ( )

%*****************************************************************************80
%
%% CALPAK_TEST0005 tests CWS_TO_JED_GPS and JED_TO_CWS_GPS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 June 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST0005\n' );
  fprintf ( 1, '  For the GPS calendar:\n' );
  fprintf ( 1, '  JED_TO_CWS_GPS: JED -> CWS.\n' );
  fprintf ( 1, '  CWS_TO_JED_GPS: CWS -> JED.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   JED (in)       CWS                            JED (out)\n' );
  fprintf ( 1, '\n' );

  jed_epoch = epoch_to_jed_gps ( );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0.0 )
      break
    end

    if ( jed_epoch <= jed1 )

      [ c2, w2, sec2 ] = jed_to_cws_gps ( jed1 );

      s2 = cws_to_s_gps ( c2, w2, sec2 );

      jed3 = cws_to_jed_gps ( c2, w2, sec2 );

      fprintf ( 1, '  %11.2f     %s  %11.2f\n', jed1, s2, jed3 );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


function fn = facenormals( V,F )
%FACENORMALS Compute face normals of triangular mesh
%   V - nverts by 3 matrix containing vertex positions
%   F - ntri by 3 matrix containing triangle vertex indices

% Get the triangle vertices
v1      = F(:, 1);
v2      = F(:, 2);
v3      = F(:, 3);

% Compute the edge vectors
e1s = V(v2, :) - V(v1, :);
e2s = V(v3, :) - V(v1, :);

% Compute cross products between edge vectors
fn    = cross(e1s, e2s, 2);

% Normalise face normals to unit length
norms = sqrt(fn(:,1).^2+fn(:,2).^2+fn(:,3).^2);
fn=fn./repmat(norms,[1 3]);

end


text:
%LMNC Levenberg-Marquardt trained feed-forward neural net classifier
%
%  [W,HIST] = LMNC (A,UNITS,ITER,W_INI,T)
%
% INPUT
%  A        Dataset
%  UNITS    Vector with numbers of units in each hidden layer (default: [5])
%  ITER     Number of iterations to train (default: inf)
%  W_INI    Weight initialization network mapping (default: [], meaning
%           initialization by Matlab's neural network toolbox)
%  T        Tuning set (default: [], meaning use A)
%
% OUTPUT
%  W        Trained feed-forward neural network mapping
%  HIST     Progress report (see below)
%
% DESCRIPTION
% A feed-forward neural network classifier with length(N) hidden layers with
% N(I) units in layer I is computed for the dataset A. Training is stopped
% after ITER epochs (at least 50) or if the iteration number exceeds twice
% that of the best classification result. This is measured by the labeled
% tuning set T. If no tuning set is supplied A is used. W_INI is used, if
% given, as network initialization. Use [] if the standard Matlab
% initialization is desired. Progress is reported in file FID (default 0).
%
% The entire training sequence is returned in HIST (number of epochs,
% classification error on A, classification error on T, MSE on A, MSE on T,
% mean of squared weights).
%
% Uses the Mathworks' Neural Network toolbox.
%
% SEE ALSO (PRTools Guide)
% MAPPINGS, DATASETS, BPXNC, NEURC, RNNC, RBNC, PRPROGRESS

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: lmnc.m,v 1.3 2007/06/15 09:58:30 duin Exp $

function [w,hist] = lmnc(varargin)

		[w,hist] = ffnc(mfilename,varargin{:});

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function value = c8_mul ( z1, z2 )

%*****************************************************************************80
%
%% C8_MUL multiplies two C8's.
%
%  Discussion:
%
%    A C8 is a complex value.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, complex Z1, Z2, the values to multiply.
%
%    Output, complex VALUE, the function value.
%
  value = z1 * z2;

  return
end