text:
function [interp1_fs, interp2_fs] = get_interp_fourier(sz, params)

% Compute the Fourier series of the interpolation function (b_d in the
% paper). The interpolation method is set using params.interpolation_method:
% - 'ideal' performs ideal reconstruction, which corresponds to a periodic
%   summation of a sinc function b_d in the spatial domain.
% - 'bicubic' uses a bicubic kernel b_d in the spatial domain.

switch lower(params.interpolation_method)
    case 'none'
        % Do nothing
        interp1_fs = ones(sz(1),1);
        interp2_fs = ones(1,sz(2));
    case 'ideal'
        % Ideal reconstruction (flat in the frequency domain)
        interp1_fs = ones(sz(1),1) / sz(1);
        interp2_fs = ones(1,sz(2)) / sz(2);
    case 'bicubic'
        % Take the truncated fourier series from the cubic spline
        a = params.interpolation_bicubic_a;
        interp1_fs = real(1/sz(1) * cubic_spline_fourier((-(sz(1)-1)/2:(sz(1)-1)/2)'/sz(1), a));
        interp2_fs = real(1/sz(2) * cubic_spline_fourier((-(sz(2)-1)/2:(sz(2)-1)/2)/sz(2), a));
    otherwise
        error('Unknown dft interpolation method');
end

if params.interpolation_centering
    % Center the feature grids by shifting the interpolated features
    % Multiply Fourier coeff with e^(-i*pi*k/N)
    interp1_fs = interp1_fs .* exp(-1i*pi / sz(1) * (-(sz(1)-1)/2:(sz(1)-1)/2)');
    interp2_fs = interp2_fs .* exp(-1i*pi / sz(2) * (-(sz(2)-1)/2:(sz(2)-1)/2));
end
if params.interpolation_windowing
    % Window the Fourier series of the interpolation basis function
    win1 = hann(sz(1)+2);
    win2 = hann(sz(2)+2);
    interp1_fs = interp1_fs .* win1(2:end-1);
    interp2_fs = interp2_fs .* win2(2:end-1)';
end

interp1_fs = single(interp1_fs);
interp2_fs = single(interp2_fs);

author: martin-danelljan

repo: Continuous-ConvOp


text:
clear

% problem specification
m = 400; n = 500; r = 10;
t = 1; esr = ceil(t*r);
fprintf('[m n r esr] = [%i, %i, %i, %i]\n',m,n,r,esr);

% data
% Xo = rand(m,r); Yo = rand(r,n);
Xo = abs(randn(m,r)); Yo = abs(randn(r,n));

d = ones(r,1).*(1:r)';
M = Xo*spdiags(d,0,r,r)*Yo;

% set solver options
opts.tol = 1e-5;
opts.maxit = 500;
opts.print = 1;

% call solver

sr = 0.5;
Omega = randsample(m*n, sr*m*n);
A = M(Omega);

%exact rank-estimate case
t0 = tic;
[X,Y,Out] = mc_nmf(A,Omega,esr,m,n,opts);
time = toc(t0);
relerr = norm(X*Y-M,'fro')/norm(M,'fro');
fprintf('RelErr = %5.2e, %4.2e Sec.\n',relerr,time);
%%
[numr,numc] = size(M);
I = randi([0 1],numr,numc);
Omega = find(I);
A = M(Omega);

%%
subplot(1,2,1),imagesc(M);
subplot(1,2,2),imagesc(X*Y);
%%

%over-estimate case
esr = ceil(1.5*r);
t0 = tic;
fprintf('[m n r esr] = [%i, %i, %i, %i]\n',m,n,r,esr);
[X,Y,Out] = mc_nmf(A,Omega,esr,m,n,opts);
time = toc(t0);
relerr = norm(X*Y-M,'fro')/norm(M,'fro');
fprintf('RelErr = %5.2e, %4.2e Sec.\n',relerr,time);

author: andrewssobral

repo: lrslibrary


text:
function [h, g] = ldpc_h2g(H,q)
% [h, g] = ldpc_h2g(H,q)
% converts tentative binary LDPC matrix H into a new matrix h
% (columns are permuted) and produces the generator matrix g
% H should be a sparse matrix in MATLAB format.
% q - Field base (power of 2) now only 2 4 8 1 32 64 128 and 256
%
% MEX file

%   Copyright (c) 1999 by Igor Kozintsev igor@ifp.uiuc.edu
%   $Revision: 1.1 $  $Date: 1999/08/23 $ - implementation for GFq

author: xiaoshaoning


text:
function [sFile, ChannelMat] = in_fopen_eeglab(DataFile, ImportOptions)
% IN_FOPEN_EEGLAB: Open an EEGLAB .set file (continuous recordings).
%
% FORMAT:
%     EEGLAB datasets can have two forms : one file (.SET) or two files (.SET/.DAT).
%     In both cases, .SET file is a Matlab matrix with the dataset header in 'EEG' field.
%         1) .SET      : Recordings are stored in the 'EEGDATA' field of the .SET matrix
%         2) .SET/.DAT : Recordings are stored in binary mode in a separate .DAT file,
%                        whose file name is stored in field 'EEG.datfile' of the .SET matrix.
%                        Format : [nbChan, nbTime*nbTrials] float32 binary matrix (Little-Endian)
%     Channel locations may be stored in the .SET file, in field 'EEG.chanlocs'

% @=============================================================================
% 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: Francois Tadel, 2008-2023

% ===== PARSE INPUTS =====
if (nargin < 2) || isempty(ImportOptions)
    ImportOptions = db_template('ImportOptions');
end

%% ===== READER HEADER =====
% Load .set file
if isstruct(DataFile)
    hdr = DataFile;
    DataFile = hdr.filename;
else
    hdr = load(DataFile, '-mat');
end
% If there is no EEG field
if ~isfield(hdr, 'EEG')
    if isfield(hdr, 'nbchan') && isfield(hdr, 'pnts') && isfield(hdr, 'trials')
        hdr = struct('EEG', hdr);
    else
        error('Invalid EEGLAB .set file: missing EEG structure.');
    end
end
% Add some information
hdr.isRaw = (~isfield(hdr.EEG, 'epoch') || isempty(hdr.EEG.epoch)) && (isfield(hdr.EEG, 'data') && ~isempty(hdr.EEG.data));
nChannels = hdr.EEG.nbchan;
nTime     = hdr.EEG.pnts;
nEpochs   = hdr.EEG.trials;

% === GET TIME ===
% if isfield(hdr.EEG, 'times') && ~isempty(hdr.EEG.times)   % Disabled the use of "times" because it was not always in milliseconds (2-Nov-2022)
% hdr.Time = hdr.EEG.times ./ 1000;
if isfield(hdr.EEG, 'srate') && ~isempty(hdr.EEG.srate)
    hdr.Time = hdr.EEG.xmin + (0:nTime-1) ./ hdr.EEG.srate;
else
    hdr.Time = linspace(hdr.EEG.xmin, hdr.EEG.xmax, nTime);
    if (nTime > 1)
        hdr.EEG.srate = 1 ./ (hdr.Time(2) - hdr.Time(1));
    else
        hdr.EEG.srate = 1000;
    end
end

% ===== LIST BAD TRIALS =====
if ~hdr.isRaw
    % Get accepted/rejected trials (with methods ICA, SIG, ELEC)
    iRejectedTrials = [];
    if isfield(hdr.EEG, 'reject')
        % Epochs rejected by ICA criteria
        if isfield(hdr.EEG.reject, 'icareject')
            iRejectedTrials = [iRejectedTrials, hdr.EEG.reject.icareject];
        end
        % Epochs rejected by single-channel criteria
        if isfield(hdr.EEG.reject, 'sigreject')
            iRejectedTrials = [iRejectedTrials, hdr.EEG.reject.sigreject];
        end
        % Epochs rejected by single-channel criteria
        if isfield(hdr.EEG.reject, 'elecreject')
            iRejectedTrials = [iRejectedTrials, hdr.EEG.reject.elecreject];
        end
    end
    iGoodTrials = setdiff(1:nEpochs, iRejectedTrials);

    % Remove trials that have at least one 'bad' value in its events list
    iBadTrials = [];
    if isfield(hdr.EEG.epoch(1), 'eventbad')
        for iTrial = 1:length(hdr.EEG.epoch)
            if (iscell(hdr.EEG.epoch(iTrial).eventbad) && any([hdr.EEG.epoch(iTrial).eventbad{:}])) || ...
               (~iscell(hdr.EEG.epoch(iTrial).eventbad) && hdr.EEG.epoch(iTrial).eventbad)
                iBadTrials(end+1) = iTrial;
            end
        end
        iGoodTrials = setdiff(iGoodTrials, iBadTrials);
    end
else
    iGoodTrials = 1;
end

% ===== GET RELEVANT CONDITIONS =====
epochNames = [];
isAllEmpty = 1;
if ~hdr.isRaw && isfield(hdr.EEG, 'event') && ~isempty(hdr.EEG.event)
    % Each trial is classified with many criteria
    % Need to ask the user along which creteria the trials should be classified
    % Get all fields of the 'event' structure
    listParam = fieldnames(hdr.EEG.event(1));
    % Remove entries that are not conditions of the events
    %listParam = setdiff(listParam, {'type','latency','urevent','obs','bad','badChan','epoch'});
    listParam = setdiff(listParam, {'urevent','obs','bad','badChan','epoch'});
    % Convert all the events to strings
    for iParam = 1:length(listParam)
        hdr.EEG.event = ConvertEventToStr(hdr.EEG.event, listParam{iParam});
        hdr.EEG.epoch = ConvertEventToStr(hdr.EEG.epoch, ['event' listParam{iParam}]);
    end
    % Keep only parameters for which the value varies in the valid trials
    iParam = 1;
    paramValues = {};
    while (iParam <= length(listParam))
        % Get all the unique values
        tmpValues = {hdr.EEG.event.(listParam{iParam})};
        % If char and not all the values are the same
%         if ~iscell(tmpValues{1}) && ~all(cellfun(@(c)isequal(c,tmpValues{1}), tmpValues))
        if ischar(tmpValues{1}) && ~all(cellfun(@(c)isequal(c,tmpValues{1}), tmpValues))
            % Latency: keep the native order
            if isequal(listParam{iParam}, 'latency')
                [tmp,I,J] = unique(tmpValues);
                paramValues{end + 1} = tmpValues(sort(I));
            % Else: sort in alphabetical/numerical order
            else
                paramValues{end + 1} = unique(tmpValues);
            end
            iParam = iParam + 1;
        % Else remove parameter
        else
            listParam(iParam) = [];
        end
    end

    % If there are different types of events that can be used to classify the epochs
    if (length(listParam) >= 1) && (ImportOptions.DisplayMessages || ~isempty(ImportOptions.EventsTypes))
        % Ask the user to select the parameters to be compared
        if ImportOptions.DisplayMessages
            ParamSelected = gui_show_dialog('Conditions selection', @panel_eeglab_cond, 1, [], listParam, paramValues);
        % Use the list of event types passed in input
        else
            % Get the list of parameters in input
            evtSplit = strtrim(str_split(ImportOptions.EventsTypes, ','));
            % Check that the types are available in the file
            iNotFound = find(~ismember(evtSplit, listParam));
            if ~isempty(iNotFound)
                disp(['BST> Error: Some of the event types in input where not found in the file: ', sprintf('%s ', evtSplit{iNotFound})]);
                evtSplit(iNotFound) = [];
            end
            % Get the condition name for each epoch
            ParamSelected = [];
            if ~isempty(evtSplit)
                % Find the selected event types
                [tmp, I, J] = intersect(evtSplit, listParam);
                % Get the possible combinations of values
                ParamSelected = panel_eeglab_cond('GetConditionCombinations', listParam(J), paramValues(J));
                params = listParam(J);
                % Build the condition name for each epoch (see panel_eeglab_cond.m)
                for iCond = 1:length(ParamSelected)
                    strCondName = '';
                    for iParam = 1:length(params)
                        if (iParam ~= 1)
                            strCondName = [strCondName '_'];
                        end
                        if ischar(ParamSelected(iCond).(params{iParam}))
                            tmpStr = strrep(file_standardize(ParamSelected(iCond).(params{iParam})), '_', '-');
                            strCondName  = [strCondName, params{iParam}, tmpStr];
                        else
                            strCondName  = [strCondName, params{iParam}, sprintf('%d', ParamSelected(iCond).(params{iParam}))];
                        end
                    end
                    ParamSelected(iCond).Name = file_standardize(strCondName);
                end
            end
        end
        % Build the epochs names
        epochNames = cell(1, nEpochs);
        % Process parameters selection
        if ~isempty(ParamSelected)
            paramsList = setdiff(fieldnames(ParamSelected), 'Name');
            % Create FileName in which this file should be saved
            for iTrial = 1:nEpochs
                epoch = hdr.EEG.epoch(iTrial);
                % Find in which condition should be classified this epoch
                isCondFound = 0;
                iCond = 1;
                while ~isCondFound && (iCond <= length(ParamSelected))
                    isOk = 1;
                    iParam = 1;
                    while (isOk && (iParam <= length(paramsList)))
                        % If value for targ parameter does not match
                        if isequal(paramsList{iParam}, 'latency')
                            isMatch = isequal(hdr.EEG.event(epoch.event(1)).latency, ParamSelected(iCond).(paramsList{iParam}));
                        elseif iscell(epoch.(['event', paramsList{iParam}]))
                            isMatch = any(cellfun(@(c)isequal(c,ParamSelected(iCond).(paramsList{iParam})), epoch.(['event', paramsList{iParam}])));
                        else
                            isMatch = isequal(epoch.(['event', paramsList{iParam}]), ParamSelected(iCond).(paramsList{iParam}));
                        end
                        if isMatch
                            iParam = iParam + 1;
                        else
                            isOk = 0;
                        end
                    end
                    if isOk
                        isCondFound = 1;
                    else
                        iCond = iCond + 1;
                    end
                end
                % Build new filename with the found condition
                if isCondFound
                    epochNames{iTrial} = ParamSelected(iCond).Name;
                else
                    epochNames{iTrial} = '';
                end
            end
            isAllEmpty = all(cellfun(@isempty, epochNames));
        % If no parameters were selected
        else
            epochNames = repmat({''}, 1, nEpochs);
            isAllEmpty = 1;
        end
    end
end

%% ===== GET DATA SOURCE =====
% EEG.data
if isfield(hdr.EEG, 'data') && ~isempty(hdr.EEG.data)
    if isfield(hdr, hdr.EEG.data)
        EEGDATA = hdr.(hdr.EEG.data);
    else
        EEGDATA = hdr.EEG.data;
    end
% EEGDATA
elseif isfield(hdr, 'EEGDATA') && ~isempty(hdr.EEGDATA)
    EEGDATA = hdr.EEGDATA;
% EEGDATA.datfile
elseif isfield(hdr.EEG, 'datfile') && ~isempty(hdr.EEG.datfile)
    EEGDATA = hdr.EEG.datfile;
% Default data file: use the same file name than for .SET file, with .DAT/.FDT extension
else
    [fPath, fBase, fExt] = bst_fileparts(DataFile);
    EEGDATA = bst_fullfile(fPath, [fBase, '.dat']);
    if ~file_exist(EEGDATA)
        EEGDATA = bst_fullfile(fPath, [fBase, '.fdt']);
    end
end
% In case of attached binary file
if ischar(EEGDATA)
    % Check binary file existence
    if ~file_exist(EEGDATA)
        % Try with the same name as the .set file
        if file_exist(strrep(DataFile, '.set', '.fdt'))
            EEGDATA = strrep(DataFile, '.set', '.fdt');
        else
            % Try with adding a full path
            [fPath, fBase, fExt] = bst_fileparts(EEGDATA);
            EEGDATA = bst_fullfile(bst_fileparts(DataFile), [fBase, fExt]);
            % File is not accessible
            if ~file_exist(EEGDATA)
                error(['EEGLAB binary file does not exist: ', EEGDATA]);
            end
        end
    end
    % Save correct filename
    hdr.EEG.data = EEGDATA;
    % Remove EEGDATA field
    if isfield(hdr, 'EEGDATA')
        hdr = rmfield(hdr, 'EEGDATA');
    end
else
    hdr.EEGDATA = EEGDATA;
end

%% ===== FILL STRUCTURE =====
% Initialize returned file structure
sFile = db_template('sfile');
% Add information read from header
sFile.filename   = DataFile;
sFile.fid        = [];
sFile.format     = 'EEG-EEGLAB';
sFile.device     = 'EEGLAB';
sFile.byteorder  = 'l';
% Properties of the recordings
if (nTime == 1)
    sFile.prop.times = hdr.Time(1) + [0, 1./hdr.EEG.srate];
else
    sFile.prop.times = [hdr.Time(1), hdr.Time(end)];
end
sFile.prop.sfreq = hdr.EEG.srate;
sFile.prop.nAvg  = 1;
sFile.header = hdr;
% Channel file
if ImportOptions.DisplayMessages
    isFixUnits = [];
elseif (ImportOptions.ChannelAlign >= 1)
    isFixUnits = 1;
else
    isFixUnits = 0;
end
ChannelMat = in_channel_eeglab_set(hdr, isFixUnits);

% === EPOCHS ===
for i = 1:nEpochs
    if ~isempty(epochNames) && ~isempty(epochNames{i})
        sFile.epochs(i).label = epochNames{i};
        sFile.epochs(i).select  = 1;
    elseif isAllEmpty
        sFile.epochs(i).label = sprintf('%s (#%d)', hdr.EEG.setname, i);
        sFile.epochs(i).select  = 1;
    else
        sFile.epochs(i).select  = 0;
    end
    sFile.epochs(i).times   = sFile.prop.times;
    sFile.epochs(i).nAvg    = 1;
    sFile.epochs(i).bad     = ~ismember(i, iGoodTrials);
    % Bad channels
    sFile.epochs(i).channelflag = ones(nChannels, 1);
    if ~hdr.isRaw && isfield(hdr.EEG, 'epoch') && ~isempty(hdr.EEG.epoch) && isfield(hdr.EEG.epoch(1), 'eventbadChan')
        % Get all the bad channels for that epoch
        eventbadChan = hdr.EEG.epoch(iTrial).eventbadChan;
        if iscell(eventbadChan)
            iBadChan = [];
            for j = 1:length(eventbadChan)
                if ~isempty(eventbadChan{j})
                    iBadChan = [iBadChan, str2num(eventbadChan{j})];
                end
            end
        else
            iBadChan = str2num(eventbadChan);
        end
        % Report the bad channels in the ChannelFlag array
        sFile.epochs(i).channelflag(iBadChan) = -1;
    end
end
% Global channel flag
if (nEpochs == 1)
    sFile.channelflag = sFile.epochs(1).channelflag;
    sFile.epochs = [];
else
    sFile.channelflag = ones(nChannels, 1);
end

% === EVENTS ====
if isfield(hdr.EEG, 'event') && ~isempty(hdr.EEG.event) && isfield(hdr.EEG.event, 'type') % && hdr.isRaw
    % Get event types
    intTypes = [];
    if ischar(hdr.EEG.event(1).type)
        listTypes = unique({hdr.EEG.event.type});
    elseif isnumeric(hdr.EEG.event(1).type)
        intTypes = unique([hdr.EEG.event.type]);
        listTypes = cell(1, length(intTypes));
        for iType = 1:length(intTypes)
            listTypes{iType} = num2str(intTypes(iType));
        end
    else
        return;
    end
    % Initialize structure
    events = repmat(db_template('event'), [1, length(listTypes)]);
    % Process each event type
    for iEvt = 1:length(listTypes)
        % Get all the event occurrences
        if ~isempty(intTypes)
            listOcc = find([hdr.EEG.event.type] == intTypes(iEvt));
        else
            listOcc = find(strcmpi({hdr.EEG.event.type}, listTypes{iEvt}));
        end
        % Get event label
        events(iEvt).label = listTypes{iEvt};
        % If no occurrences: skip
        if isempty(listOcc)
            continue;
        end
        % Get epochs indices
        if ~hdr.isRaw && (isfield(hdr.EEG.event(listOcc(1)), 'epoch') && ~isempty(hdr.EEG.event(listOcc(1)).epoch))
            events(iEvt).epochs = [hdr.EEG.event(listOcc).epoch];
        else
            events(iEvt).epochs = ones(1, length(listOcc));
        end
        % Get samples
        if isfield(hdr.EEG.event(listOcc(1)), 'latency')
            allSmp = {hdr.EEG.event(listOcc).latency};
        elseif isfield(hdr.EEG.event(listOcc(1)), 'sample')
            allSmp = {hdr.EEG.event(listOcc).sample};
        else
            disp(['EEGLAB> Missing fields "latency" or "sample" in event "', hdr.EEG.event(listOcc(1)).type, '".']);
        end
        % Convert to values if available as strings
        iChar = find(cellfun(@ischar, allSmp));
        if ~isempty(iChar)
            allSmp(iChar) = cellfun(@str2num, allSmp(iChar), 'UniformOutput', 0);
        end
        % Remove empty latencies
        iEmpty = find(cellfun(@isempty, allSmp));
        if ~isempty(iEmpty)
            [allSmp{iEmpty}] = deal(1);
        end
        samples = round([allSmp{:}]);
        % Add durations if there are more than one sample
        if isfield(hdr.EEG.event(listOcc), 'duration') && ~ischar(hdr.EEG.event(listOcc(1)).duration)
            allDur = [hdr.EEG.event(listOcc).duration];
            if any(allDur > 1) && (length(samples) == length(allDur))
                samples(2,:) = samples + allDur;
            end
        end
        % For epoched files: convert events to samples local to each epoch
        if ~hdr.isRaw
            nSamples = round((sFile.prop.times(2) - sFile.prop.times(1)) .* sFile.prop.sfreq) + 1;
            samples = samples - (events(iEvt).epochs - 1) * nSamples + sFile.prop.times(1) * sFile.prop.sfreq - 1;
        end
        % Compute times
        events(iEvt).times = samples ./ sFile.prop.sfreq;
        % Additional fields
        events(iEvt).channels = [];
        events(iEvt).notes    = [];
    end
    % Save structure
    sFile.events = events;
end

end

%% ===== CONVERT EVENTS TO STRING =====
function s = ConvertEventToStr(s, param)
    % For each event
    for i = 1:length(s)
        % Double values: convert to strings
        if isnumeric(s(i).(param))
            s(i).(param) = val2str(s(i).(param));
        % Cell array of doubles
        elseif iscell(s(i).(param))
            for iCell = 1:length(s(i).(param))
                if isnumeric(s(i).(param){iCell})
                    s(i).(param){iCell} = val2str(s(i).(param){iCell});
                end
            end
        end
    end
end

function str = val2str(val)
    str = '';
    for iVal = 1:length(val)
        if (str > 1)
            str = [str ' '];
        end
        str = [str, num2str(val(iVal))];
    end
end

author: brainstorm-tools


text:
function OutputFile = db_add(iTarget, InputFile, isReload)
% DB_ADD: Create a new file in a given study.
%
% USAGE:  db_add(iStudy/iSubject, InputFile, isReload)
%         db_add(iStudy/iSubject, sNew     , isReload)
%         db_add(iStudy/iSubject)             : Read a variable from Matlab workspace

% @=============================================================================
% 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: Francois Tadel, 2011-2022

%% ===== GET INPUT FILE =====
if (nargin < 3) || isempty(isReload)
    isReload = 1;
end
% USAGE:  db_add(iTarget)
if (nargin < 2) || isempty(InputFile)
    % Get variable from workspace
    [InputFile, varname] = in_matlab_var([], 'struct');
    if isempty(InputFile)
        OutputFile = [];
        return
    end
end
% USAGE:  db_add(iTarget, sMat)
if isstruct(InputFile)
    % Get structure
    sMat = InputFile;
    InputFile = [];
    % Detect file type
    fileType = file_gettype(sMat);
    % Add sub-category for timefreq files
    fileSubType = '';
    if ismember(fileType, {'timefreq', 'ptimefreq'})
        if ismember(lower(sMat.Method), {'corr', 'cohere', 'granger', 'spgranger', 'plv', 'plvt', 'aec', 'pte','henv'})
            if (length(sMat.RefRowNames) > 1) && (isequal(sMat.RowNames, sMat.RefRowNames) || isequal(sMat.RowNames, sMat.RefRowNames'))
                fileSubType = ['connectn_', lower(sMat.Method), '_'];
                if (size(sMat.TF,1) == length(sMat.RefRowNames)*length(sMat.RowNames))
                    sMat.Options.isSymmetric = 0;
                else
                    sMat.Options.isSymmetric = 1;
                end
            else
                fileSubType = ['connect1_', lower(sMat.Method), '_'];
                sMat.Options.isSymmetric = 0;
            end
        elseif ismember(lower(sMat.Method), {'pac', 'dpac', 'tpac'}) && isfield(sMat,'sPAC')
            if isfield(sMat.sPAC, 'DirectPAC') && ~isempty(sMat.sPAC.DirectPAC)
                fileSubType = 'pac_fullmaps_';
            else
                fileSubType = 'dpac_fullmaps_';
            end
        else
            fileSubType = [lower(sMat.Method), '_'];
        end
    end
    % Add stat thresholding method
    if ismember(fileType, {'pdata', 'presults', 'ptimefreq', 'pmatrix'})
        if isfield(sMat, 'Correction') && ~isempty(sMat.Correction)
            fileSubType = [fileSubType, lower(sMat.Correction), '_'];
        end
    end
    % Add zscore tag
    if isfield(sMat, 'ZScore') && ~isempty(sMat.ZScore)
        if isfield(sMat.ZScore, 'abs') && isequal(sMat.ZScore.abs, 1)
            fileSubType = [fileSubType 'abs_zscore_'];
        else
            fileSubType = [fileSubType 'zscore_'];
        end
    end
    % Some files cannot be imported
    if ismember(fileType, {'brainstormsubject', 'brainstormstudy', 'unknown'})
        bst_error('This structure cannot be imported in the database.', 'Add file to database', 0);
        return;
    end
    % Surfaces subtypes
    if ismember(fileType, {'fibers', 'fem'})
        fileType = 'tess';
        fileSubType = [fileType, '_'];
    end
    isAnatomy = ismember(fileType, {'subjectimage', 'tess'});
    % Spikes: file tag is data
    if strcmp(fileType, 'spike')
        fileType = 'data';
        fileSubType = '0ephys_';
    end
    % Create a new output filename
    c = clock;
    strTime = sprintf('%02.0f%02.0f%02.0f_%02.0f%02.0f', c(1)-2000, c(2:5));
    OutputFile = [fileType '_' fileSubType strTime '.mat'];

% USAGE:  db_add(iTarget, InputFile)
else
    % Get full file path and type
    [InputFileFull, fileType, isAnatomy] = file_fullpath(InputFile);
    % Output filename
    [fPath, fBase] = bst_fileparts(InputFile);
    OutputFile = [fBase, '.mat'];
    % Read file
    sMat = load(InputFileFull);
end

%% ===== OUTPUT FILENAME =====
% Get protocol folders
ProtocolInfo = bst_get('ProtocolInfo');
% Full output filename
if isAnatomy
    % Get destination study
    sSubject = bst_get('Subject', iTarget);
    % Build full filename
    OutputFileFull = file_unique(bst_fullfile(ProtocolInfo.SUBJECTS, bst_fileparts(sSubject.FileName), OutputFile));
    OutputFile = file_short(OutputFileFull);
else
    % Get destination study
    sStudy = bst_get('Study', iTarget);
    % Build full filename
    OutputFileFull = file_unique(bst_fullfile(ProtocolInfo.STUDIES, bst_fileparts(sStudy.FileName), OutputFile));
    OutputFile = file_short(OutputFileFull);
end

%% ===== REPLACE EXISTING FILES =====
% Check for files that can exist in only one version
if ismember(fileType, {'subjectimage', 'channel', 'noisecov', 'ndatacov'})
    % Get file to replace
    delfile = '';
    switch(fileType)
        case 'subjectimage'
%             if ~isempty(sSubject.Anatomy)
%                 delfile = bst_fullfile(ProtocolInfo.SUBJECTS, sSubject.Anatomy(1).FileName);
%             end
        case 'channel'
            if ~isempty(sStudy.Channel)
                delfile = bst_fullfile(ProtocolInfo.STUDIES, sStudy.Channel(1).FileName);
            end
        case 'noisecov'
            if ~isempty(sStudy.NoiseCov) && ~isempty(sStudy.NoiseCov(1).FileName)
                delfile = bst_fullfile(ProtocolInfo.STUDIES, sStudy.NoiseCov(1).FileName);
            end
        case 'ndatacov'
            if ~isempty(sStudy.NoiseCov) && (length(sStudy.NoiseCov) >= 2) && ~isempty(sStudy.NoiseCov(2).FileName)
                delfile = bst_fullfile(ProtocolInfo.STUDIES, sStudy.NoiseCov(2).FileName);
            end
    end
    % Replace file
    if ~isempty(delfile)
        % Ask for user confirmation
        isdel = java_dialog('confirm', ['Replace existing ' fileType ' file?'], 'Add file to database');
        if ~isdel
            OutputFile = [];
            return
        end
        % Delete existing file
        file_delete(delfile, 1);
    end
end

%% ===== ADD COMMENT TAG =====
matVer = 'v7';
if isAnatomy
    % Get subject structure
    sSubject = bst_get('Subject', iTarget);
    % Update comment with a file tag, to make it unique
    switch (fileType)
        case 'subjectimage'
            % Nothing to do: file is replaced anyway
        case {'tess', 'cortex', 'scalp', 'outerskull', 'innerskull', 'fibers', 'fem'}
            sMat.Comment = file_unique(sMat.Comment, {sSubject.Surface.Comment});
    end
else
    % Add comment if missing
    if isempty(sMat.Comment)
        [tmp__, sMat.Comment] = bst_fileparts(InputFile);
    end
    % Get study structure
    sStudy = bst_get('Study', iTarget);
    % Update comment with a file tag, to make it unique
    switch (fileType)
        case 'data'
            sMat.Comment = file_unique(sMat.Comment, {sStudy.Data.Comment});
            matVer = 'v6';
        case 'headmodel'
            sMat.Comment = file_unique(sMat.Comment, {sStudy.HeadModel.Comment});
        case 'results'
            if ~isempty(sMat.DataFile)
                iRes = find(file_compare(sMat.DataFile, {sStudy.Result.DataFile}));
                if isempty(iRes)
                    iRes = find(cellfun(@isempty, {sStudy.Result.DataFile}));
                end
                if ~isempty(iRes)
                    sMat.Comment = file_unique(sMat.Comment, {sStudy.Result(iRes).Comment});
                end
            else
                sMat.Comment = file_unique(sMat.Comment, {sStudy.Result.Comment});
            end
            matVer = 'v6';
        case 'stat'
            sMat.Comment = file_unique(sMat.Comment, {sStudy.Stat.Comment});
            matVer = 'v6';
        case 'dipoles'
            sMat.Comment = file_unique(sMat.Comment, {sStudy.Dipoles.Comment});
        case 'timefreq'
            sMat.Comment = file_unique(sMat.Comment, {sStudy.Timefreq.Comment});
            matVer = 'v6';
        case 'matrix'
            sMat.Comment = file_unique(sMat.Comment, {sStudy.Matrix.Comment});
            matVer = 'v6';
    end
end

%% ===== ADD NEW FILE TO DATABASE =====
% Progress bar
isProgress = bst_progress('isVisible');
if ~isProgress
    bst_progress('start', 'Add file to database', 'Saving new file...');
else
    bst_progress('text', 'Saving new file...');
end
% Save new file
bst_save(OutputFileFull, sMat, matVer);
% Close progress bar
if ~isProgress
    bst_progress('stop');
end
% Reload output study
if isReload
    if isAnatomy
        db_reload_subjects(iTarget);
    else
        db_reload_studies(iTarget);
    end
    panel_protocols('UpdateTree');
    panel_protocols('SelectNode', [], OutputFileFull);
end


text:
function G = curl(F)
%CURL  Numerical surface curl.
%   G = CURL(F) returns the numerical surface curl of the SPHEREFUNV F.
%   This only makes mathematical sense for a tangential vector field.
%
% See also SPHEREFUNV/DIVERGENCE, SPHEREFUN/GRADIENT, SPHEREFUNV/VORTICITY.

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

% Empty check.
if ( isempty(F) )
    G = spherefunv;
    return
end

Fc = F.components;

% The following derivatives are all tangential derivatives.

% First component: d/dy(fz) - d/dz(fy)
gx = diff(Fc{3}, 2) - diff(Fc{2}, 3);

% Second component: d/dz(fx) - d/dx(fz)
gy = diff(Fc{1}, 3) - diff(Fc{3}, 1);

% Third component: d/dx(fy) - d/dy(fx)
gz = diff(Fc{2}, 1) - diff(Fc{1}, 2);

G = spherefunv(gx,gy,gz);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function pHW=x2p_(xHW,layer,convnet)
Stride=convnet.targetStride(layer);
centerStart=convnet.targetCenter(layer);
pHW=centerStart+(xHW-1).*Stride;
end

repo: interpretableCNN


text:
function RHS = convectionTvdRHSSpherical1D(u, phi, FL)
% This function uses the upwind scheme to discretize a 1D
% convection term in the form \grad (u \phi) where u is a face vactor
% It also returns the x and y parts of the matrix of coefficient.
%
% SYNOPSIS:
%   [M, RHS] = convectionTvdTermCylindrical1D(u, phi, FL)
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% extract data from the mesh structure
Nr = u.domain.dims(1);
% DXp = u.domain.cellsize.x(2:end-1);
dx = 0.5*(u.domain.cellsize.x(1:end-1)+u.domain.cellsize.x(2:end));
% r = u.domain.cellcenters.x;
rf = u.domain.facecenters.x;
RHS = zeros(Nr+2, 1);
psi_p = zeros(Nr+1,1);
psi_m = zeros(Nr+1,1);

% extract the velocity data
% note: size(ux) = [1:m+1, 1:n] and size(uy) = [1:m, 1:n+1]
ux = u.xvalue;

% calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
% P is 3:Nr+2
% W is 2:Nr+1
% WW is 1:Nr
dphi_p = (phi.value(2:Nr+2)-phi.value(1:Nr+1))./dx;
rp = dphi_p(1:end-1)./fsign(dphi_p(2:end));
psi_p(2:Nr+1) = 0.5*FL(rp).*(phi.value(3:Nr+2)-phi.value(2:Nr+1));
psi_p(1) = 0.0; % left boundary will be handled explicitly

% calculate the upstream to downstream gradient ratios for u<0 (- ratio)
% P is 3:Nr+2
% W is 2:Nr+1
% WW is 1:Nr
rm = dphi_p(2:end)./fsign(dphi_p(1:end-1));
psi_m(1:Nr) = 0.5*FL(rm).*(phi.value(1:Nr)-phi.value(2:Nr+1));
psi_m(Nr+1) = 0.0; % right boundary will be handled explicitly

% reassign the east, west, north, and south velocity vectors for the
% code readability
ue = ux(2:Nr+1);		uw = ux(1:Nr);
re = rf(2:Nr+1);     rw = rf(1:Nr);

% find the velocity direction for the upwind scheme
ue_min = min(ue,0);	ue_max = max(ue,0);
uw_min = min(uw,0);	uw_max = max(uw,0);

% calculate the TVD correction term
RHS(2:Nr+1) = -(1./(1/3*(rf(2:Nr+1).^3-rf(1:Nr).^3))).*(re.^2.*(ue_max.*psi_p(2:Nr+1)+ue_min.*psi_m(2:Nr+1))- ...
              rw.^2.*(uw_max.*psi_p(1:Nr)+uw_min.*psi_m(1:Nr)));

end

function phi_out = fsign(phi_in)
% This function checks the value of phi_in and assigns an eps value to the
% elements that are less than or equal to zero, while keeping the signs of
% the nonzero elements
    phi_out = (abs(phi_in)>=eps).*phi_in+eps*(phi_in==0)+eps*(abs(phi_in)<eps).*sign(phi_in);
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function b_sorted = comp_b_sorted(R, bF)

[~,Rfp,~,Rpp] = seg_R(R);

% Optimized bP
bP_opt = -Rpp\Rfp'*bF;
b_sorted = [bF;bP_opt];

end

author: yorgoon

repo: minimum-snap-geometric-control


text:
function [x, infos] = snpa(V, num_col, in_options)
% Successive Nonnegative Projection Algorithm (variant with f(.) = ||.||^2)
%
%       At each step of the algorithm, the column of X maximizing ||.||_2 is
%       extracted, and X is updated with the residual of the projection of its
%       columns onto the convex hull of the columns extracted so far.
%
% Inputs:
%       matrix      V
%       num_col     number of columns to be extracted.
%       options     options
%           normalization: 1: scale the columns of X so that they sum to one,
%                           hence matrix H will satisfy the assumption above for any
%                           nonnegative separable matrix X.
%                          0: the default value for which no scaling is
%                           performed. For example, in hyperspectral imaging, this
%                           assumption is already satisfied and normalization is not
%                           necessary.
%
% Output:
%       w           solution of w
%           K        : index set of the extracted columns.
%           H        : optimal weights, that is, H = argmin_{Y >= 0} ||X-X(:,K)Y||_F
%       infos       information
%
% References:
%       N. Gillis,
%       "Successive Nonnegative Projection Algorithm for Robust Nonnegative Blind Source Separation,"
%       SIAM J. on Imaging Sciences 7 (2),
%       pp. 1420-1450, 2014.
%
%
% This file is part of NMFLibrary.
%
% This file has been ported from
%   SNPA.m at https://gitlab.com/ngillis/nmfbook/-/tree/master/algorithms
%   by Nicolas Gillis (nicolas.gillis@umons.ac.be)
%
% Change log:
%
%   June. 21, 2022 (Hiroyuki Kasai): Added initialization module.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.disp_freq = 1;
    local_options.normalize = 0;
    local_options.relerr = 1e-6;
    local_options.inner_max_epoch = 500;
    local_options.inner_nnls_alg = 'fpgm';
    local_options.special_stop_condition = @(epoch, infos, options, stop_options) spna_stop_func(epoch, infos, options, stop_options);

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize
    method_name = 'SNPA';
    i = 0;
    grad_calc_count = 0;
    stop_options = [];

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    % initialize for this algorithm
    options.max_epoch = num_col+1;
    if options.normalize == 1
        % normalize the columns of V of which colum is L1-norm = 1
        V = normalize_W(V, 1);
    end

    U = zeros(m, num_col);
    K = zeros(1, num_col);
    H = zeros(num_col, n);
    normV0 = sum(V.^2);
    nVmax = max(normV0);
    normR = normV0;
    VtUK = [];
    UKtUK = [];

    % set for nnls subsolver
    nnls_options = [];
    nnls_options.verbose = 0;
    nnls_options.inner_max_epoch = options.inner_max_epoch;
    nnls_options.algo = options.inner_nnls_alg;

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, eye(m), V, [], options, [], i, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('%s: Epoch = 0000, cost = %.16e, optgap = %.4e\n', method_name, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main loop
    i = i + 1;
    while true

        % check stop condition
        stop_options.normR = normR;
        stop_options.nVmax = nVmax;
        [stop_flag, reason, max_reached_flag] = check_stop_condition(i, infos, options, stop_options);
        if stop_flag
            display_stop_reason(i, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % select the column of the residual R with largest l2-norm
        [a, ~] = max(normR);

        % check ties up to 1e-6 precision
        b = find((a-normR)/a <= 1e-6);

        % In case of a tie, select column with largest norm of the input matrix X
        if length(b) > 1
            [~, d] = max(normX0(b));
            b = b(d);
        end

        % update the index set, and extracted column
        K(i) = b;
        U(:, i) = V(:, b);

        % update MtUJ
        VtUK = [VtUK, V' * U(:, i)];

        % update UJtUJ
        if i == 1
            UtUi = [];
        else
            UtUi = U(:, 1:i-1)' * U(:, i);
        end
        UKtUK = [UKtUK, UtUi ; UtUi', U(:, i)' * U(:, i)];

        % update residual
        if i == 1
            % Fast gradient method for min_{y in Delta} ||M(:, i)-M(:,J)y||
            [H, ~, ~, ~] = nnls_solver(V, V(:, K(1:i)), nnls_options);
        else
            H(:, K(i)) = 0;
            h = zeros(1,n); h(K(i)) = 1;
            H = [H; h];
            nnls_options.init = H;
            [H, ~, ~, ~] = nnls_solver(V, V(:, K(1:i)), nnls_options);
        end

        % update the norm of the columns of the residual without computing it explicitely.
        if i == 1
            normR = normV0 - 2 * ((VtUK') .* H) + (H .* (UKtUK*H));
        else
            normR = normV0 - 2 * sum((VtUK') .* H) + sum(H .* (UKtUK*H));
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % measure gradient calc count
        grad_calc_count = grad_calc_count + m*n;

        % update epoch
        i = i + 1;

        % store info
        infos = store_nmf_info(V, U(:, 1:i-1), H, [], options, infos, i-1, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, i-1, infos, options);

    end

    x.K = K;
    x.U = U;
    x.H = H;

end

function [stop_flag, reason, rev_infos] = spna_stop_func(epoch, infos, options, stop_options)

    stop_flag = false;
    reason = [];
    rev_infos = [];

    normR = stop_options.normR;
    nVmax = stop_options.nVmax;

    if sqrt(max(normR)/nVmax) < options.relerr
        stop_flag = true;
        reason = sprintf('precision reached: sqrt(max(normR)/nVmax) = %.4e < options.relerr = %.4e\n', sqrt(max(normR)/nVmax), options.relerr);
        return;
    end

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [align] = mamarama(in,config)

Cst = in.Cst;
Rot = in.Rot;

Octave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
if ~Octave,
  drawscene(in.Xe,Cst',Rot,41,'cloud','Graphical Output Validation: View from top or bottom (no sRt)',config.cal.cams2use);
end

% definition of the absolute world frame (in cm)

%% for hummingbird cage
%
%cam(3).C = [168, 173, 185]';
%cam(4).C = [142, 76, 89]';
%cam(1).C = [81, 72, 165]';
%cam(2).C = [22,  15, 176]';
%cam(5).C = [53.5, 175, 78]';

% definition of the absolute world frame (in mm) (Doug's cube)

% cam(1).C = [10, 390, 600]';
% cam(2).C = [460,  380, 600]';
% cam(3).C = [430, 115, 580]';
% cam(4).C = [0, 0, 610]';
% cam(5).C = [240, 150, 610]';

% definition of the absolute world frame (in mm) (windtunnel June 2005)

%cam(1).C = [1510, 180, 500+322]';
%cam(2).C = [1140,  260, 450+322]';
%cam(3).C = [1040, 25, 480+322]';
%cam(4).C = [490, 160, 390+322]';
%cam(5).C = [900, 25, 400+322]';

% the above numbers (from June 2005) transformed to acheive better
% agreement (July 3, 2005)

%cam(1).C = [1585.0266, 157.1808, 725.3687]';
%cam(2).C = [1221.0045, 230.7819, 717.1524]';
%cam(3).C = [1123.2774, -2.1424, 731.8758]';
%cam(4).C = [580.8019, 124.5745, 705.3737]';
%cam(5).C = [978.8653, 3.3878, 663.1311]';

% cal July 5, 2005

%cam(1).C = [1710, 375, 820]';
%cam(2).C = [1110, 316, 800]';
%cam(3).C = [1110, 100, 466]';
%cam(4).C = [-380, 290, 875]';
%cam(5).C = [440, 15, 680]';

% touch-up July 6, 2005 (5c)

%cam(1).C = [1853.5325, 228.1160, 971.5365]';
%cam(2).C = [1202.0386, 199.5437, 953.7752]';
%cam(3).C = [1183.3920, 41.6797, 537.2034]';
%cam(4).C = [-406.8365, 229.6378, 1073.0449]';
%cam(5).C = [460.2840, -60.0974, 778.1218]';

% touch-up Aug 22, 2005 (touch-up from A applied to B)

%cam(1).C = cam(1).C - [130.0*0.83, 0., 110.0*0.83]';
%cam(2).C = cam(2).C - [130.0*0.83, 0., 110.0*0.83]';
%cam(3).C = cam(3).C - [130.0*0.83, 0., 110.0*0.83]';
%cam(4).C = cam(4).C - [130.0*0.83, 0., 110.0*0.83]';
%cam(5).C = cam(5).C - [130.0*0.83, 0., 110.0*0.83]';

% touch-up Jan 27, 2006 ( hacked-in estimate )

%cam(1).C = cam(1).C - [130.0*0.83+20.0, 0., 110.0*0.83]';
%cam(2).C = cam(2).C - [130.0*0.83+10.0, 0., 110.0*0.83]';
%cam(3).C = cam(3).C - [130.0*0.83+10.0, 0., 110.0*0.83]';
%cam(4).C = cam(4).C - [130.0*0.83, 0., 110.0*0.83]';
%cam(5).C = cam(5).C - [130.0*0.83, 0., 110.0*0.83]';

% new Feb 14, 2006

%cam(1).C = [1760, 190, 850]';
%cam(2).C = [1060, 260, 785]';
%cam(3).C = [1070, 70, 470 ]';
%cam(4).C = [-380, 230, 873]';
%cam(5).C = [380, -10, 670]';

% moved cam 2006 03 07 19:45

%cam(1).C = [1760, 190, 850]';
%cam(2).C = [1060, 260, 785]';
%cam(3).C = [1070, 190, 510 ]';
%cam(4).C = [-380, 230, 873]';
%cam(5).C = [380, 120, 670]';

if 1,

% 2006 03 31

%cam(1).C = [640, -105, 170]';
%cam(2).C = [1220, 172, 912]';
%cam(3).C = [475, 172, 580 ]';
%cam(4).C = [-350, 172, 932]';
%cam(5).C = [320, -105, 170]';

% 2006 04 03e

%cam(3).C = [380, 240, 530 ]';

% 2006 04 04a

%cam(4).C = [-115, 220, 800]';
end

% 2006 09 14

%cam(1).C = [685, 220, 775]';
%cam(2).C = [1220, 172, 912]';
%cam(3).C = [380, 240, 530 ]';
%cam(4).C = [-115, 220, 800]';
%cam(5).C = [180, 240, 740]';

%% 2006 10 23 (from DLT)
%cam(1).C = [ 651.7963301   169.95615265  749.51356735]';
%cam(2).C = [ 1090.40810501   190.82014771   866.8980674 ]';
%cam(3).C = [ 327.41874025  237.43022637  539.80942371]';
%cam(4).C = [-241.364468    194.27934314  767.35249515]';
%cam(5).C = [ 147.20028399  171.45774441  725.08522715]';

% 2006 12 19 riverside hummingbirds
% estimates from DLT
%cam(1).C = [976.22314128 -2289.85157108   489.78104927]';
%cam(2).C = [-17278.94752856   1794.85211148  15855.51116479]';
%cam(3).C = [ 2266.45466953  2590.2964507    -74.53159994]';
%cam(4).C = [-329.01069449  796.46308431  624.04191004]';

% hand measurements
%ft2mm = 12*25.4;
%cam(1).C = [ 1350 -1580 540]';
%cam(2).C = [ -450 770 (750+5*ft2mm) ]';
%cam(3).C = [ 1500 -1600 2250]';
%cam(4).C = [-280 820 570]';

% 2006 12 01 DLT in flydra WT
%cam(1).C = [867.37446261  198.97798549  832.51824736]';
%cam(2).C = [1194.97461953 206.718389208 898.66723201]';
%cam(3).C = [459.669418616 168.068949056 713.890827]';
%cam(4).C = [-186.0833207 207.892635698 775.298571621]';
%cam(5).C = [115.230367611 202.584140956 904.723766629]';

% 2007 10 03 DLT in flydra WT
%cam(1).C = [ 985.827871     182.49267738  1134.97921258]';
%cam(2).C = [362.27817754   189.63204145  1134.23310277]';
%cam(3).C = [721.10852967  167.31210242  932.32564817]';
%cam(4).C = [-330.47764957   168.38518854  1014.81801646]';
%cam(5).C = [1774.53765317   144.79725152   928.10954671]';

% 2007 11 16 mamarama
%cam(1).C = ([20     25  30]*25.4)';
%cam(2).C = ([26   36  33]*25.4)';
%cam(3).C = ([26  -36  31]*25.4)';
%cam(4).C = ([74   2  30]*25.4)';
%cam(5).C = ([17   -25   30]*25.4)';

% 2007 12 18 mamarama
% cam(1).C = [130   450    705]';
% cam(2).C = [0     660    665]';
% cam(3).C = [1460  640    665]';
% cam(4).C = [685   1865   710]';
% cam(5).C = [1315  480    715]';

% % 2008 03 12 mamarama
% cam(1).C = [710   1350    785]';
% cam(2).C = [795   1255    785]';
% cam(3).C = [-610   1255    790]';
% cam(4).C = [-505   1355    790]';
% cam(5).C = [110   -355     785]';
% cam(6).C = [-850   305     785]'; %mama07.2
% cam(7).C = [-810   190     785]';
% cam(8).C = [1020   305     785]'; %mama08.2
% cam(9).C = [980   190     785]';

% 2008 04 17 mamarama
cam(1).C = [645   -265    760]';    %mama01
cam(2).C = [640   1250    760]';    %mama02
cam(3).C = [-455  1260    780]';   %mama03
cam(4).C = [-450  -295    780]';   %mama04
cam(5).C = [100   495     805]';   %mama05
cam(6).C = [85    1460    780]';   %mama06_0
cam(7).C = [85    -470    780]';   %mama06_1
cam(8).C = [-745  920     775]';   %mama07_0
cam(9).C = [-685  -65     775]';    %mama07_1
cam(10).C = [890   1030    775]';  %mama08_0
cam(11).C = [935   80     775]';   %mama08_1

% of the similarity computation

[align.simT.s, align.simT.R, align.simT.t]  = estsimt([Cst'],[cam(:).C]);
[align.P, align.X] = align3d(in.Pe,in.Xe,align.simT);
% save aligned data
if 1 % SAVE_STEPHI | SAVE_PGUHA
	[align.Cst,align.Rot] = savecalpar(align.P,config);
end

if ~Octave,
  drawscene(align.X,align.Cst',align.Rot,61,'cloud','Graphical Output Validation: Aligned data',config.cal.cams2use);

  set(gca,'CameraTarget',[0,0,0]);
  set(gca,'CameraPosition',[0,0,1]);

  figure(61),
  % print -depsc graphevalaligned.eps
  eval(['print -depsc ', config.paths.data, 'topview.eps'])

  drawscene(align.X,align.Cst',align.Rot,62,'cloud','Graphical Output Validation: Aligned data',config.cal.cams2use);

  set(gca,'CameraTarget',[0,0,0.9]);
  set(gca,'CameraPosition',[2,0,0.9]);

  %figure(62),
  % print -depsc graphevalaligned.eps
  %eval(['print -depsc ', config.paths.data, 'sideview.eps'])
end

return

author: strawlab

repo: MultiCamSelfCal


text:
function f = uplus(f)
%UPLUS   Unary plus of a TRIGTECH.
%   UPLUS(F) = F.
%
% See also UMINUS, PLUS.

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

% This function does not change anything.

end


text:
function y = hsv2hsl(x)

%hsv2hsl - Convert hue-saturation-value colors to hue-saturation-luminance.
%
%  USAGE
%
%    y = hsv2hsl(x)
%
%    x              Nx3 RGB matrix (values in [0..1])

% Copyright (C) 2013 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Check number of parameters
if nargin < 1,
  error('Incorrect number of parameters (type ''help hsv2hsl'' for details).');
end
if ~isdmatrix(x) || size(x,2) ~= 3 || any(x(:)<0) || any(x(:)>1),
  error('Incorrect HSL matrix (type ''help hsv2hsl'' for details).');
end

% Convert HSV to HSL
h = x(:,1);
s = x(:,2);
v = x(:,3);

H = h;
S = s .* v;
L = (2-s).*v;
in = L <= 1;
S(in) = S(in) ./ L(in);
S(in&v==0) = 0;
S(~in) = S(~in) ./ (2-L(~in));
S(~in&s==0&v==1) = 0;
L = L/2;

y = Clip([H S L],0,1);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [a]=times(b,c,varargin)
%A=B.*C
%   [A]=TIMES(B,C) Computes Hadamard (elementwise) product of two
%   TT-matrices
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
n = b.n;
m = b.m;
b = b.tt;
c = c.tt;
a = b.*c;
a = tt_matrix(a, n, m);
end

author: oseledets

repo: TT-Toolbox


text:
function name = PTKGetUserVisibleNameForContext(context)
    % PTKGetUserVisibleNameForContext.
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2014.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    context = char(context);

    if strcmp(context, 'Lungs')
        name = 'Both lungs';
    elseif strcmp(context, 'LeftLung')
        name = 'Left lung';
    elseif strcmp(context, 'RightLung')
        name = 'Right lung';
    elseif strcmp(context, 'LeftUpperLobe')
        name = 'Left upper lobe';
    elseif strcmp(context, 'LeftLowerLobe')
        name = 'Left lower lobe';
    elseif strcmp(context, 'RightUpperLobe')
        name = 'Right upper lobe';
    elseif strcmp(context, 'RightLowerLobe')
        name = 'Right lower lobe';
    elseif strcmp(context, 'RightMiddleLobe')
        name = 'Right middle lobe';
    else
        name = context;
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
%testGrowthExpMatch tests the functionality of all the components of growthExpMatch
%
%Procedure to run growthExpMatch:
%(1) obtain all input files (ie. model, CompAbr, and KEGGID are from BiGG, KEGGList is from KEGG website)
%(2) remove desired reaction from model with removeRxns, or set the model
%       on a particular Carbon or Nitrogen source
%(3) create an SUX Matrix by using the function MatricesSUX =
%       generateSUXMatrix(model,Dictionary, KEGGFilename,Compartment)
%(4) run it through growthExpMatch using [solution,b,solInt]=Smiley(MatricesSUX)
%(5) solution.importedRxns contains the solutions to all iterations (in
%       this particular test, we removed reaction ENO from model and
%       obtained 'R00658_f' which is the KEGGID for ENO
%
%   Joseph Kang 11/16/09
%   Adaptions to CI - Thomas Pfau Okt 2017

global CBTDIR

oriFolder = pwd;

%moves to testing folder that contains testGrowthExpMatch
test_folder = fileparts(which('testGrowthExpMatch.m'));
cd(test_folder);

%load Model
model = getDistributedModel('ecoli_core_model.mat');

%removes reaction ENO
disp('------------------------------------')
disp('Removing reaction, ENO, from test model:')
model = removeRxns(model, 'ENO');

%moves to folder w/ input files
% w = what('testing');
% p = w.path;
% cd(p);
d = load([test_folder filesep 'Dictionary.mat']);
KEGGFilename = 'testTest_KEGG_Reaction_List.lst';

%Test for the default solvers
solverPkgs = {'gurobi', 'glpk'};

for k = 1:length(solverPkgs)
    fprintf('   Running solveCobraLPCPLEX using %s ... ', solverPkgs{k});

    % change the COBRA solver
    solverOKLP = changeCobraSolver(solverPkgs{k}, 'LP', 0);
    solverOKMILP = changeCobraSolver(solverPkgs{k}, 'MILP', 0);

    if solverOKLP && solverOKMILP
        %runs growthExpMatch and obtains solution
        [solution]=growthExpMatch(model, KEGGFilename,'[c]', 1, d.dictionary);

        %if R00658_f is solution result, returns a positive answer, else negative
        assert(isequal(solution.importedRxns, {{'R00658_f'}}))
    end
end
%perform cleanup
if exist([test_folder filesep 'GEMLog_solution_1.mat'],'file')
    delete('GEMLog_solution_1.mat')
end
if exist([test_folder filesep 'CobraMILPSolver.log'],'file')
    delete('CobraMILPSolver.log')
end
if exist([test_folder filesep 'GEMLog.txt'],'file')
    delete('GEMLog.txt')
end
close all

cd(oriFolder);

author: opencobra

repo: cobratoolbox


text:
function proposals= calcEdgeBoxesForIm( input, ebconfig )
rmpath(genpath([pwd '/rigor']));
rmpath(genpath([pwd '/mcg/']));
if(isstr(input))
	im = im2double(imread(input));
else
        im = im2double(input); % Just to make input consistent
end
if(size(im, 3) == 1)
        im=repmat(im,[1,1,3]);
end
% load pre-trained edge detection model and set opts
if(~exist(ebconfig.modelPath))
    fprintf('Path to model does not exist. Please make sure you give a proper full path\n');
    return;
end

model = load(ebconfig.modelPath);
model = model.model;
model.opts.multiscale = 0;
model.opts.sharpen = 2;
model.opts.nThread = 4;

%Write code to set options
% call edgeBoxes() to get back options

opts = edgeBoxes();
bbs=edgeBoxes(im,model,opts);
if(isfield((ebconfig.opts),'numProposals'))
	numProposals=ebconfig.opts.numProposals;
        if(size(bbs,1)>=numProposals)
       	        bbs=bbs(1:numProposals);
       	else
               	fprintf('Only %d proposals were generated for input image.\n',size(bbs,1));
        end
end
%edges boxes produces baoxes as "[x y w, h]"
%we convert to [x y x+w y+h]==[xmin ymin xmax ymax]
boxes=bbs(:,1:4);
boxes=[boxes(:,1) boxes(:,2) boxes(:,1)+ boxes(:,3) boxes(:,2)+boxes(:,4)];
proposals.boxes= boxes;
proposals.scores = bbs(:,5);

addpath(genpath([pwd '/mcg/']));
addpath(genpath([pwd '/rigor']));
end

author: Cloud-CV

repo: object-proposals


text:
function value = fall ( x, n )

%*****************************************************************************80
%
%% FALL computes the falling factorial function [X]_N.
%
%  Discussion:
%
%    The number of "injections" or 1-to-1 mappings from
%    a set of N elements to a set of M elements is [M]_N.
%
%    The number of permutations of N objects out of M is [M}_N.
%
%    The Stirling numbers of the first kind can be used
%    to convert a falling factorial into a polynomial, as follows:
%
%      [X]_N = S^0_N + S^1_N * X + S^2_N * X^2 + ... + S^N_N X^N.
%
%  Formula:
%
%    [X]_N = X * ( X - 1 ) * ( X - 2 ) * ... * ( X - N + 1 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X, the argument of the falling factorial
%    function.
%
%    Input, integer N, the order of the falling factorial function.
%    If N = 0, FALL = 1, if N = 1, FALL = X.  Note that if N is
%    negative, a "rising" factorial will be computed.
%
%    Output, integer VALUE, the falling factorial function.
%
  value = 1;

  arg = x;

  if ( 0 < n )

    for i = 1 : n
      value = value * arg;
      arg = arg - 1;
    end

  elseif ( n < 0 )

    for i = -1 : - 1 : - n
      value = value * arg;
      arg = arg + 1;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [set_ar, varxh, klifit] = ARShat_misd(ng,xg,L,lag_max,set_ar_start);

%ARSHAT_MISD AR models of increasing order from measurements with missing data
%   [set_ar, varxh, klifit] = ARShat_misd(ng,xg,[Lmin Lmax],lag_max) estimates
%   autoregressive models from measurements containing missing data for
%   orders Lmin, ..., Lmax.
%   lag_max contains the maximum lag for the ARFIL algorithm for all orders
%   that are considered.
%
%   NOTE: The mean of the signal is NOT subtracted.
%
%   [set_ar, varxh, klifit] = ARShat_misd(ng,xg,[Lmin Lmax],set_ar_start)
%   The cell array set_ar_start contains starting models for all model
%   orders considered.
%
%   Output:
%   set_ar is a cell array containing all estimated models;
%   varxh  is the estimated signal variance;
%   klifit is the fit for each model in terms of the Kullback-Leibler Index
%          (kli) for missing data.
%
%   See also: ARHAT_MISD, SIG2AR_MISD.

%S. de Waele, August 2001.

showwb = 0; %1 = Show waitbar

nseg = length(ng);
%-------------------------------------------------------------
%Reading input arguments
%-------------------------------------------------------------
mstart = 0;
if exist('set_ar_start'),
    mstart = 1;
else
    mstart = 0;
end

Lmin = L(1);
Lmax = L(2);
nk = Lmax-Lmin+1;

%-------------------------------------------------------------
%Estimation
%-------------------------------------------------------------
sumvars = 0;
for seg = 1:nseg
   sumvars = sumvars + std(xg{seg})^2;
end
varxh = sumvars/nseg

nobsi = length(xg);
set_ar = cell(nk,1);
klifit  = zeros(nk,1);

if showwb,
	hwb = waitbar(0/nk,[int2str(Lmin) ' to ' int2str(Lmax) ' AR models - ARSirr-tan']);
end

rch_prev_calc = 0; %Currently, no previous model has been calculated
for ik = 1:nk,
   k = Lmin+ik-1;
   if k == 0,
      %White noise estimate
      klifit(ik) = ARMLfit([],ng,xg,0,lag_max(ik));
      rch = [];
   else	%if k == 0
      %Determination starting point.
      %Option 1:	white noise estimate
      rc_00 = zeros(1,k);
      % klifit_wn = fit_cL_irr(rc_start,ng,xg,fs);
      klifit_wn = ARMLfit([],ng,xg,rc_00,lag_max(ik));
      rc_start = rc_00;
      klifit_start(1) = klifit_wn;
      %Option 2:	Append zero to previous model
      rc_old0 = [rch_prev 0];
      klifit_old0 = ARMLfit([],ng,xg,rc_old0,lag_max(ik));
      rc_start(2,:) = rc_old0;
      klifit_start(2) = klifit_old0;
      if mstart,
         %Option 3:	User-provided Starting Model.
         [dummy rc_sm] = ar2arset(set_ar_start{ik});
         rc_sm = rc_sm(2:end);
         klifit_sm = ARMLfit([],ng,xg,rc_sm,lag_max(ik));
         rc_start(3,:) = rc_sm;
         klifit_start(3) = klifit_sm;
      end
      %Selection of inition conditions with best fit.
      [klifit_start, i_min] = min(klifit_start);
      rc_start = rc_start(i_min,:)
      methods = {'allzero','previous+0','User-provided'};
      disp(['Selected starting point: ' methods{i_min}])
      %Optimization starting in selected point
      [rch,fit_temp] = ARhat_misd(ng,xg,rc_start,[],lag_max(ik));
      klifit(ik) = fit_temp;
      if klifit(ik) > klifit_start
         warning(['Optimization did not yield better model for k = ' int2str(k) '.'])
         rch = rc_start;
         klifit(ik) = klifit_start;
      end
   end %if k == 0
   set_ar{ik} = rc2arset([1 rch]);
   rch_prev = rch; %For the next round
   rch_prev_calc = 1;
   if showwb,
      waitbar(ik/nk)
   end
end
if showwb,
	close(hwb)
end

author: Sable

repo: mcbench-benchmarks


text:
vectors_file_name = add_data_dir_base('GoogleNews_vectors_norm.mat');

if word2vec_sqrt
	vectors_file_name = fname_concat(vectors_file_name, '_sqrt');
end

is_sampled = true;
%is_sampled = false;

% without normalizing we had numeric problems in the LMM
normalize  = true;

if is_sampled
    vectors_file_name = fname_concat(vectors_file_name, '_sampled');
end

output(1, 'reading vectors file %s\n', vectors_file_name);
load(vectors_file_name);
output(1, 'done\n');

%TargetDim = 256;
TargetDim = 300;
%TargetDim = 100;

output(1, 'Use PCA to reduce to dimension %d\n', TargetDim);

[pc,~,latent] = princomp(vectors_');
sumlatent = cumsum(latent) ./ sum(latent);

pc = pc(:,1:TargetDim);
trans_matrix = pc';
vectors_ = trans_matrix * vectors_;

[dim, num_vectors] = size(vectors_)

if normalize
    % normalize
    output(1, 'normalizing...\n');
    for i = 1:num_vectors
      vectors_(:,i) = norma(vectors_(:,i));
    end
end

output(1, 'saving to file...\n');

post_fix = sprintf('_pca_%d', TargetDim);
pca_vectors_file_name = fname_concat(vectors_file_name, post_fix);
trans_matrix_file_name = add_data_dir_base(strcat('trans_matrix', post_fix));

save(pca_vectors_file_name, 'vectors_', '-v7.3');
output(1, 'saved to file: %s\n', pca_vectors_file_name);

save(trans_matrix_file_name, 'trans_matrix', '-v7.3');
output(1, 'saved to file: %s\n', trans_matrix_file_name);

clear vectors_;
clear trans_matrix;

author: layumi

repo: Image-Text-Embedding


text:
% Initialization
clear all; close all; clc;

a1 = [ 2;  1];
a2 = [ 2; -1];
a3 = [-1;  2];
a4 = [-1; -2];

A = [a1 a2 a3 a4]';

b = ones(4,1);
span = [-2 2 -2 2];
spx.graphics.plot.polyhedron(A, b, span);

author: indigits

repo: sparse-plex


text:
% Johann Diep (johann.diep@esa.int) - August 2021
%
% This script builds the structure for the visual dataset.

%% Conversion Constants

G =  9.80665;       % Gravity constant, m/s^2
G2MSS = G;          % g to m/s^2
MSS2G = (1/G);      % m/s^2 to g

D2R = (pi/180);     % degrees to radians
R2D = (180/pi);     % radians to degrees

KT2MS = 0.514444;   % knot to m/s
MS2KMH = 3.6;       % m/s to km/h

%% OpenVINS PLANETARY DATASET ERROR PROFILE

load('./canada-planetary-data/data//gnss_planetary_r.mat');
load('./canada-planetary-data/data//imu_planetary.mat');

visual_data_planetary_pose = readtable('./external-data/canada/pose_data.csv');
visual_data_planetary_twist = readtable('./external-data/canada/twist_data.csv');
visual_data_planetary_posecov = readtable('./external-data/canada/posecov_data.csv');
visual_data_planetary_twistcov = readtable('./external-data/canada/twistcov_data.csv');

%   visual, visual data structure
%         t: Vx1 time vector (seconds)
%       lat: Vx1 latitude (radians).
%       lon: Vx1 longitude (radians).
%         h: Vx1 altitude (m).
%       vel: Vx3 NED velocities (m/s).
%     covvm: Vx36 velocity and position covariance matrices (m^2 and m^2/s^2).
%      covv: Vx9 velocity covariance matrices (m^2/s^2).
%      larm: 3x1 lever arm from IMU to camera (x-fwd, y-right, z-down) (m).
%      freq: 1x1 sampling frequency (Hz).
%       eps: 1x1 time interval to compare current IMU time to current GNSS time vector (s).

Cutoff = 1500;

% parsing the timestamps of the measurements
Seconds = table2array(visual_data_planetary_pose(Cutoff:end,4));
Nanoseconds = table2array(visual_data_planetary_pose(Cutoff:end,5));
Timestamps = Seconds + Nanoseconds/10^(9);
visual_planetary.t = Timestamps - Timestamps(1);

% position and velocity preprocessing
Positions = table2array(visual_data_planetary_pose(Cutoff:end,6:7))-table2array(visual_data_planetary_pose(Cutoff,6:7));
Heights = table2array(visual_data_planetary_pose(Cutoff:end,8))-table2array(visual_data_planetary_pose(Cutoff,8));
RotationAngle = deg2rad(-100);
RotationMatrix = [cos(RotationAngle),-sin(RotationAngle);sin(RotationAngle),cos(RotationAngle)];
Positions_rotated = RotationMatrix * Positions';

[visual_planetary.lat, visual_planetary.lon, visual_planetary.h] = ned2geodetic(Positions_rotated(2,:)',Positions_rotated(1,:)',-Heights, ...
    gnss_planetary_r.lat(1),gnss_planetary_r.lon(1),gnss_planetary_r.h(1),wgs84Ellipsoid,'radians');

Velocities = table2array(visual_data_planetary_twist(Cutoff:end,1:2));
Velocities_rotated = RotationMatrix * Velocities';
visual_planetary.vel(:,1:2) = Velocities_rotated';
visual_planetary.vel(:,3) = -table2array(visual_data_planetary_twist(Cutoff:end,3));

% parsing the covariances
Count = 1;
for i = Cutoff:size(table2array(visual_data_planetary_posecov),1)
        PoseCovariance = table2array(visual_data_planetary_posecov(i,:));
        TwistCovariance = table2array(visual_data_planetary_twistcov(i,:));

        PoseCovariance = reshape(PoseCovariance,[6,6])';
        TwistCovariance = reshape(TwistCovariance,[6,6])';

        PositionCovariance = PoseCovariance(1:3,1:3);
        VelocityCovariance = TwistCovariance(1:3,1:3);

        % [~,Std] = corrcov(PositionCovariance);
        % Std_x(i) = Std(1);
        % Std_y(i) = Std(2);
        % Std_z(i) = Std(3);

        IntermediateCovariance = [VelocityCovariance,zeros(3,3);zeros(3,3),PositionCovariance];

        visual_planetary.covv(Count,:) = reshape(VelocityCovariance',[1,9]);
        visual_planetary.covvm(Count,:) = reshape(IntermediateCovariance',[1,36]);
        Count = Count + 1;
end

visual_planetary.larm = [0,0,0]'; % no lever arm needed in case of OpenVINS
visual_planetary.eps = mean(diff(imu_planetary.t))*20; % rule of thumb for choosing eps

% downsample
visual_planetary.t = downsample(visual_planetary.t,33);
visual_planetary.lat = downsample(visual_planetary.lat,33);
visual_planetary.lon = downsample(visual_planetary.lon,33);
visual_planetary.h = downsample(visual_planetary.h,33);
visual_planetary.vel = downsample(visual_planetary.vel,33);
visual_planetary.covv = downsample(visual_planetary.covv,33);
visual_planetary.covvm = downsample(visual_planetary.covvm,33);

visual_planetary.freq = 1/mean(diff(visual_planetary.t)); % estimating the frequency

save('./canada-planetary-data/data/visual_planetary.mat','visual_planetary');

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function [Vo, Ro, No, V] = generate_syntheticdata(F, N, K, noise_level, rho)
%
%     % Vo
%     W = randn(F,K);
%     W = max(W, 0);
%     W = min(W, 1);
%
%     H = randn(K,N);
%     H = max(H, 0);
%     H = min(H, 1);
%
%     Vo = W * H;
%
%     % noise
     No = noise_level*randn(F,N);
%
%     % add outlier
%     [V, Ro] = add_outlier(rho, F, N, Vo);
%
%     % V
%     V = V + No;
%     %V = max(V, 0);
%     %V = min(V,1);

    %
    sigma2 = 1 / sqrt(K);
    HN = makedist('Normal', 'mu', 0, 'sigma', sqrt(sigma2));
    Vo_n = random(HN, F, N);
    Vo = abs(Vo_n) ;
    Vo = min(Vo, 1);

    nu = rho;
    nu_tilda = 0.1;
    I = nu * N;
    card = nu_tilda * F;
    Ro = zeros(F,N);
    if rho > 0
        for i = 1 : N
            n_before = 0;
            if i < I
                for f = 1 : card
                    c = randi(F);
                     Ro(c,i) =  1 + (1+1)*rand(1, 1);
                     n = nnz(Ro(:,i));
                     if n_before == n
                         f = f - 1;
                     end
                     n_before = n;
                end

            end
        end
    end

    % V
    V = Vo + Ro + No;
    %V = max(V, 0);
    %V = min(V,1);

    index = find(V<0);
    V(index) = 0;

end

repo: SGDLibrary


text:
classdef ConnecCoordFromInterpAndMesh < handle

    properties (Access = public)
        coord
        connec
    end

    properties (Access = private)
        mesh
        interp
    end

    methods (Access = public)

        function obj = ConnecCoordFromInterpAndMesh(cParams)
            obj.init(cParams)
        end

        function compute(obj)
            m = obj.mesh;
            shapesInVarNodes = obj.computeShapesInVariableNodes(m);
            xPointsMesh      = m.coord;
            Tmesh            = m.connec;

            xPoints = zeros(1,obj.interp.ndime);
            Tinterp = zeros(m.nelem,obj.interp.nnode);

            inode = 1;
            for ielem = 1:m.nelem
                for inodeVar = 1:obj.interp.nnode
                    xNode = zeros(1,obj.interp.ndime);
                    for inodeMesh = 1:m.nnodeElem
                        node = Tmesh(ielem,inodeMesh);
                        shapes = shapesInVarNodes(inodeVar,inodeMesh);
                        xNode = xNode + shapes*xPointsMesh(node,:);
                    end

                    node = obj.findPointInList(xNode,xPoints);

                    if isempty(node)
                        xPoints(inode,:) = xNode;
                        Tinterp(ielem,inodeVar) = inode;
                        inode = inode+1;
                    else
                        Tinterp(ielem,inodeVar) = node;
                    end
                end
            end
            obj.coord = xPoints;
            obj.connec = Tinterp;
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.mesh   = cParams.mesh;
            obj.interp = cParams.interpolation;
        end

        function shapes = computeShapesInVariableNodes(obj,mesh)
            interpMesh = Interpolation.create(mesh,'LINEAR');
            nNodeMesh = interpMesh.nnode;
            nNodeVar  = obj.interp.nnode;
            shapes    = zeros(nNodeVar,nNodeMesh);
            nodesVar  = obj.interp.pos_nodes;
            for inodeVar = 1:obj.interp.nnode
                nodesPoints = nodesVar(inodeVar,:);
                interpMesh.computeShapeDeriv(nodesPoints')
                shapes(inodeVar,:) = interpMesh.shape;
            end
        end

    end

    methods (Access = private, Static)

        function ind = findPointInList(node,xPoints)
            match = true(size(xPoints,1),1);
            for idime = 1:size(node,2)
                match = match & xPoints(:,idime) == node(idime);
            end
            ind = find(match);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function M = euclideancomplexfactory(m, n)
% Returns a manifold struct to optimize over complex m-by-n matrices.
%
% function M = euclideancomplexfactory(m, n)
%
% Returns M, a structure describing the vector space of complex m-by-n
% matrices, as a manifold for Manopt.
%
% The complex plane is here viewed as R^2. The inner product between two
% m-by-n matrices A and B is given by: real(trace(A'*B)). This choice
% guides the proper definition of gradient and Hessian for this geometry.
% This is not the classical Euclidean inner product for complex matrices;
% it is a real inner product.
%
% See also: euclideanfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, April 7, 2015.
% Contributors:
% Change log:

    if ~exist('n', 'var') || isempty(n)
        n = 1;
    end

    M.name = @() sprintf('Vector space C^(%dx%d)', m, n);

    M.dim = @() 2*m*n;

    M.inner = @(x, d1, d2) real(d1(:)'*d2(:));

    M.norm = @(x, d) norm(d, 'fro');

    M.dist = @(x, y) norm(x-y, 'fro');

    M.typicaldist = @() sqrt(m*n);

    M.proj = @(x, d) d;

    M.egrad2rgrad = @(x, g) g;

    M.ehess2rhess = @(x, eg, eh, d) eh;

    M.tangent = M.proj;

    M.exp = @exp;
    function y = exp(x, d, t)
        if nargin == 3
            y = x + t*d;
        else
            y = x + d;
        end
    end

    M.retr = M.exp;

	M.log = @(x, y) y-x;

    M.hash = @(x) ['z' hashmd5([real(x(:)) ; imag(x(:))])];

    M.rand = @() (randn(m, n) + 1i*randn(m, n))/sqrt(2);

    M.randvec = @randvec;
    function u = randvec(x) %#ok<INUSD>
        u = randn(m, n) + 1i*randn(m, n);
        u = u / norm(u, 'fro');
    end

    M.lincomb = @matrixlincomb;

    M.zerovec = @(x) zeros(m, n);

    M.transp = @(x1, x2, d) d;

    M.pairmean = @(x1, x2) .5*(x1+x2);

    mn = m*n;
    M.vec = @(x, u_mat) [real(u_mat(:)) ; imag(u_mat(:))];
    M.mat = @(x, u_vec) reshape(u_vec(1:mn), [m, n]) + 1i*reshape(u_vec((mn+1):end), [m, n]);
    M.vecmatareisometries = @() true;

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function y = Fault_decision_wg2_t1(u)
global st10

% Using model defined by the structure st10 to make decision about type1 wg2 fault

sker=st10.x2sup+(abs(u))'*abs(u)*ones(st10.Nlsup,1)-2*st10.xsup*abs(u);
y=(st10.w)'*exp(-sker./(2*(st10.sigma).^2))+st10.b;
if y>=0
    y=1;
else
    y=0;
end
end


text:
%  Figure 10.32      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% Fig. 10.32 Script for computation of the root locus of yaw response,
% with no compensation, direct feedback

f =[-0.0558   -0.9968    0.0802    0.0415;
    0.5980   -0.1150   -0.0318         0 ;
   -3.0500    0.3880   -0.4650         0 ;
         0    0.0805    1.0000         0] ;
g =[0.0073;
   -0.4750;
    0.1530;
         0];
h = [0     1     0     0];
j =[0];
% the equations of the actuator:

na=[0 10];
da=[1 10];
[fa,ga,ha,ja]=tf2ss(na,da);

% the equations of the aircraft with actuator:

[fp,gp,hp,jp]=series(fa,ga,ha,ja,f,g,h,j);
% the uncompensated rootlocus
 hold off; clf
 rlocus(fp,-gp,hp,jp);
 axis([-2, 2, -1.5, 1.5]);
 grid;
 title('Fig. 10.32 Root locus of the aircraft with positive feedback')


text:
% Test file for singfun/flipud.m

function pass = test_flipud(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebfunpref();
end

% Generate a few random points to use as test values.
seedRNG(6178);
d = 10;
x = 2*(1-10^(-d)) * rand(100, 1) - (1-10^(-d));

% The order of the exponents:
a = 0.64;
b = -0.64;
c = 1.28;
d = -1.28;

%%
% Spot-check derivatives for a couple of functions.

% fractional root at the left endpoint
data.exponents = [a 0];
data.singType = {'root', 'none'};
f = singfun(@(x) (1+x).^a.*exp(x), data, pref);
g = flipud(f);
vals_df = feval(g, x);
flip_exact = @(x) (1-x).^a.*exp(-x);
vals_exact = feval(flip_exact, x);
err = vals_df - vals_exact;
pass(1) = (norm(err, inf) < 1e1*eps*norm(vals_exact, inf));

% fractional pole at the left endpoint
data.exponents = [d 0];
data.singType = {'sing', 'none'};
f = singfun(@(x) (1+x).^d.*sin(x), data, pref);
g = flipud(f);
vals_df = feval(g, x);
flip_exact = @(x) -(1-x).^d.*sin(x);
vals_exact = feval(flip_exact, x);
err = vals_df - vals_exact;
pass(2) = (norm(err, inf) < 10*eps*norm(vals_exact, inf));

% fractional root at the right endpoint
data.exponents = [0 c];
data.singType = {'none', 'root'};
f = singfun(@(x) (1-x).^c.*cos(x), data, pref);
g = flipud(f);
vals_df = feval(g, x);
flip_exact = @(x) (1+x).^c.*cos(x);
vals_exact = feval(flip_exact, x);
err = vals_df - vals_exact;
pass(3) = (norm(err, inf) < 1e1*eps*norm(vals_exact, inf));

% fractional pole at the right endpoint
data.exponents = [0 b];
data.singType = {'none', 'sing'};
f = singfun(@(x) (1-x).^b.*(x.^5), data, pref);
g = flipud(f);
vals_df = feval(g, x);
flip_exact = @(x) -(1+x).^b.*(x.^5);
vals_exact = feval(flip_exact, x);
err = vals_df - vals_exact;
pass(4) = (norm(err, inf) < 1e1*eps*norm(vals_exact, inf));

% a combination of fractional pole and fractional root
data.exponents = [b c];
data.singType = {'sing', 'root'};
f = singfun(@(x) (1+x).^b.*sin(x).*(1-x).^c, data, pref);
g = flipud(f);
vals_df = feval(g, x);
flip_exact = @(x) -(1-x).^b.*sin(x).*(1+x).^c;
vals_exact = feval(flip_exact, x);
err = vals_df - vals_exact;
pass(5) = (norm(err, inf) < 1e1*eps*norm(vals_exact, inf));

%%
% Verify that calling flipud() gives the reasonably accurate answer as direct
% construction.

data.exponents = [b b];
data.singType = {'sing', 'sing'};
f = singfun(@(x) (1+x).^b.*sin(2*x).*(1-x).^b, data, pref);
g = flipud(f);
vals_df = feval(g, x);
flip_exact = @(x) -(1-x).^b.*sin(2*x).*(1+x).^b;
vals_exact = feval(flip_exact, x);
err = vals_df - vals_exact;
pass(6) = (norm(err, inf) < 1e1*eps*norm(vals_exact, inf));

%%
% Check higher-order derivatives.

data.exponents = [a b];
data.singType = {'root', 'sing'};
f = singfun(@(x) (1+x).^a.*sin(x).*(1-x).^b, data, pref);
df2 = flipud(f);
vals_df2 = feval(df2, x);
df2_exact = @(x) -(1-x).^a.*sin(x).*(1+x).^b;
vals_exact = feval(df2_exact, x);
err = vals_df2 - vals_exact;
pass(7) = (norm(err, inf) < 5*eps*norm(vals_exact, inf));

end

author: chebfun

repo: chebfun


text:
clear
close all
clc;

%time
dt=0.01;

%coeffs
Cg=50.;
Cd=.3;
minSpeed=5;

name='greenpyr.tiff'
sanitizedname=[strrep(name,'.','') '_']
mkdir([sanitizedname 'outputdir'])
%green
%greenheight=double(imread('testgreen.tiff'))/256*1.7;
%greenheight=double(imread('greencone.tiff'))/256*.3;
%greenheight=double(imread('greenconeinv.tiff'))/256*.3;
greenheight=fliplr(rot90(rot90(double(imread(name))/256*.4)));
greenheightPLOT=fliplr(rot90(rot90(double(imread(name))/256*.2)));
[gradX,gradY]=gradient(greenheight(:,:,1));
forceFun=@(R,RL) (RL-R)/dt*Cd + [interp2(-gradX,R(1,:),R(2,:));interp2(-gradY,R(1,:),R(2,:))]*Cg;
%forceFun=@(R,RL) [interp2(-gradX,R(1,:),R(2,:));interp2(-gradY,R(1,:),R(2,:))];
l=size(greenheight,1);

%pin
holeLoc=size(greenheight)'/2;
holeRadius=3;

%graphics
%sph=sphere(30)
alt=30.85;
az=-39+90;
cmapgreen=flipud([181,228,138; 160,220,104; 124,216,87; 81,212,76; 52,188,67; 37,167,60; 32,154,61; 1,114,56; 0,86,19]/255);
stickX=[holeLoc(1),holeLoc(1)];
stickY=[holeLoc(2),holeLoc(2)];
stickZ=[greenheight(floor(holeLoc(1)),floor(holeLoc(1))),greenheight(floor(holeLoc(1)),floor(holeLoc(1)))+100];
flagX=[stickX;stickX+40];
flagY=[stickY;stickY];
flagZ=[stickZ(2),stickZ(2)-30;stickZ(2),stickZ(2)-30];
flagC=zeros(2,2,3);
flagC(:,:,1)=.9;

% balls
startLoc=[64;64]*3.5;

%angles=(45)*pi/180;
%speeds=400;
%angles=linspace((45-15)*pi/180,(45+15)*pi/180,10);%-pi/4;
angles=linspace((45-180)*pi/180,(45+180)*pi/180,100);%-pi/4;
speeds=linspace(minSpeed,300,100);

numParticles=length(speeds)*length(angles);

speeds=reshape(speeds,1,1,length(speeds));
angleVectors=[cos(angles);sin(angles)];
startConditions=bsxfun(@times,speeds,angleVectors);
startConditions=reshape(startConditions,2,length(speeds)*length(angles));

r=repmat(startLoc,1,numParticles);
rl=r-startConditions*dt;

rstart=r;
rlstart=rl;

%sequentialhits staging
%each ball laucnhes after 1/4 second (300 iters/4) = 80 ticks
launchTicks=80

i=0;
while sum(sum(r~=rl))
    %energy verlet + drag
    i=i+1;
    rn=2*r-rl+(forceFun(r,rl))*dt^2;
    rl=r;
    r=rn;

    %%%LAUNCH CONDITIONS
    if i<50
        r=rstart;
        rl=rlstart;
    end

    %%%STOP CONDITIONS
    %static friction
    dr=rl-r;
    s=sqrt(dr(1,:).^2+dr(2,:).^2)/dt;
    %walls
    haltBallsEdge=(r>l)|(r<1);
    %in the hole
    distToHole=bsxfun(@minus,holeLoc,rl);%use rl not r to check if the ball is in the hole so it can't escape
    distToHole=sqrt(distToHole(1,:).^2+distToHole(2,:).^2);
    %stop those in stop conditions
    holeBalls=(distToHole<holeRadius);
    rl(:,holeBalls)=repmat(holeLoc,1,sum(holeBalls));
    haltBalls=(haltBallsEdge(1,:)|haltBallsEdge(2,:))|(s<minSpeed)|(distToHole<holeRadius);
    r(:,haltBalls)=rl(:,haltBalls);

    if mod(i,10)==0
        figure(1)
        surf(greenheightPLOT,'edgecolor','none')
        colormap(cmapgreen)
        set(gca,'DataAspectRatio',[1 1 1])
        hold on
        scatter3(r(1,:),r(2,:),3+greenheightPLOT(sub2ind(size(greenheightPLOT), floor(r(2,:)),floor(r(1,:)))),'wo','filled')
        plot3(stickX,stickY,stickZ,'k-')
        surf(flagX,flagY,flagZ,flagC,'edgecolor','none')
        hold off
        az=az+.1;
        view([az alt])
%         imshow(greenheight,[min(min(greenheight)),max(max(greenheight))],'colormap',colormap('parula'))
        hold on
        %rectangle('Position',[holeLoc(1)-holeRadius holeLoc(2)-holeRadius holeRadius*2 holeRadius*2],'Curvature',[1,1],'FaceColor','black');
        %scatter(r(1,:),r(2,:),'w.')
        hold off
        xlim([0,512])
        ylim([0,512])
        %title(num2str(s(1)))
        set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
        saveas(gcf,[[sanitizedname 'outputdir'] '\' num2str(i/10,'%05i') '.png'])
    end
end

%phase space

figure(2)
maxDist=256;
maxColor=reshape([1,0,0],1,1,3);
minColor=reshape([1,1,1],1,1,3);
dists=fliplr(flipud(reshape(distToHole,length(angles),length(speeds))'));
imshow(dists,[holeRadius,maxDist],'colormap',colormap('jet'))
hold on
h=imshow(repmat(ones(size(dists)),1,1,3));
hold off
set(h,'AlphaData',dists<holeRadius);
ylabel('Putt speed')
xlabel('putt angle')

author: BrianHaidet

repo: AlphaPhoenix


text:
% pop_fmrib_qrsdetect() -  GUI for fmrib_qrsdetect:
%   Detect QRS events from ECG channel and add them to
%   EEGLAB data event structure.
%
% Usage:
%    >> [EEGOUT, COM] = pop_fmrib_qrsdetec(EEG) % for gui
% or >> [EEGOUT, COM] = pop_fmrib_qrsdetec(EEG,ecgchan,'ename','yes')
% for command line. 'no' or 'yes' can be used in the last argument:
%
% ecgchan: is the ecg channel number
% 'ename': event name to call the found peaks
% 'yes' or 'no': do you want to delete ecg channel?
%
% Inputs: EEG data structure
%
% GUI:
%   Enter ECG Channel Number: ECG Data channel to extract QRS peaks from.
%   What would you like to name the QRS events?: Event type to assign to
%       the QRS peaks detected.
%   Delete ECG channel when finished (yes/no)?: I think this is obvious!
%
%
% Author: Rami Niazy, FMRIB Centre, University of Oxford.
%
%
%
% Copyright (c) 2004 University of Oxford.
%

% Copyright (C) 2004 University of Oxford
% Author:   Rami K. Niazy, FMRIB Centre
%           rami@fmrib.ox.ac.uk
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% Last update:
% Oct 26, 2004

function [EEGOUT, command] = pop_fmrib_qrsdetect(EEG,varargin);

EEGOUT = [];
command = '';

% check inputs
%--------------------
nargchk(1,4,nargin);

if nargin>1 & nargin<4
    error('incorrect input arguments','Please see help.');
end

if isempty(EEG.data)
    errordlg2(strvcat('Error: no data'), 'Error');
    error('pop_fmrib_qrsdetect error: no data','fmrib_qrsdetect() error!'); return;
end

% check for DSP Toolbox
%----------------------
if ~exist('firls')
   error('QRS event detection requires the signal processing toolbox.','fmrib_qrsdetect() error!');
end

%Launch user input GUI and get inputs
%------------------------------------
if nargin ==1
	prompt={'Enter ECG channel number:','What would you like to name the QRS events?','Delete ECG channel when finished (yes/no)?'};
	def={'','qrs','no'};
	dlgtitle='pop_fmrib_qrsdetect()';
	helpfunc='fmrib_qrsdetect';
	answer=inputdlg2(prompt,dlgtitle,1,def,helpfunc);
    if ~isempty(answer)
		ecgchan=str2num(char(answer(1)));
		qrsevent=char(answer(2));
		ecgdel=char(answer(3));
    else
        return;
    end
else
    ecgchan=varargin{1};
    qrsevent=varargin{2};
    ecgdel=varargin{3};
end

%Detect QRS
%----------
fprintf('\nFinding QRS Peaks...\n');
Peaks=fmrib_qrsdetect(EEG,ecgchan);

%Add Peaks to events
%-------------------
fprintf('Writing QRS events to event structure...\n');
for index = 1:length(Peaks)
	EEG.event(end+1).type  = qrsevent;
	EEG.event(end).latency = Peaks(index);
    if EEG.trials > 1 | isfield(EEG.event, 'epoch');
        EEG.event(end).epoch = 1+floor((EEG.event(end).latency-1) / EEG.pnts);
    end
end

if EEG.trials > 1
    EEG = pop_editeventvals( EEG, 'sort', {  'epoch' 0 'latency', [0] } );
else
    EEG = pop_editeventvals( EEG, 'sort', { 'latency', [0] } );
end

if isfield(EEG.event, 'urevent')
    EEG.event = rmfield(EEG.event, 'urevent');
end

EEG = eeg_checkset(EEG, 'eventconsistency');
EEG = eeg_checkset(EEG, 'makeur');

% delete ECG channel
% ---------------
if strcmp(lower(ecgdel), 'yes')
	EEG = pop_select(EEG, 'nochannel', ecgchan);
end

%return command
%--------------
EEGOUT=eeg_checkset(EEG);
command = sprintf('EEG=pop_fmrib_qrsdetect(EEG,%d,''%s'',''%s'')',ecgchan,qrsevent,ecgdel);
return;

author: PatternRecognition

repo: OpenBMI


text:
function [data_len, winlen, winstep] = mrcg_extract( audio_dir, varargin )

    %UNTITLED2 이 함수의 요약 설명 위치
    %   자세한 설명 위치
    [noisy_speech, audio_sr] = audioread(audio_dir);

    if isempty(varargin)
        y_label = zeros(length(noisy_speech), 1);
    else
        label_dir = varargin{1};
        load(label_dir);
    end

    mkdir './sample_data'
    mkdir './sample_data/Labels'

    save_dir = './sample_data';
    name_mrcg = [save_dir, '/mrcg'];
    name_label = [save_dir, '/Labels/label'];

%     noisy_speech = noisy_speech(1:1000000,:);
%     y_label = y_label(1:1000000,:);

    mrcg_mat = MRCG_features(noisy_speech, audio_sr)';
    winlen             = ceil(audio_sr*25*0.001);	%window length (default : 25 ms)
    winstep            = ceil(audio_sr*10*0.001);	%window step (default : 10 ms)
    train_mean = mean(mrcg_mat,1);
    train_std = std(mrcg_mat,1,1);

    %% Read label

    framed_label = Truelabel2Trueframe( y_label, winlen, winstep );
    num = 0;

    if (length(mrcg_mat) > length(framed_label))
        binary_saver( name_mrcg, mrcg_mat(1:length(framed_label), :), num );
        binary_saver( name_label, framed_label, num );
        data_len = length(framed_label);
    else
        binary_saver( name_mrcg, mrcg_mat, num );
        binary_saver( name_label, framed_label(1:length(mrcg_mat), 1), num );
        data_len = length(mrcg_mat);
    end

    %% Save global normalization factor

    save([save_dir, '/normalize_factor'], 'train_mean', 'train_std');

    disp('MRCG extraction is successifully done.')

end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function sparse_grid_composite_test ( )

%*****************************************************************************80
%
%% SPARSE_GRID_COMPOSITE_TEST tests the SPARSE_GRID_COMPOSITE library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 August 2009
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_COMPOSITE_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the SPARSE_GRID_COMPOSITE library.\n' );
%
%  Count number of points in sparse rule from DIM_MIN to DIM_MAX, LEVEL_MAX_MAX.
%
  sparse_grid_composite_test01 ( 1, 6, 0, 5 );
  sparse_grid_composite_test01 ( 6, 10, 0, 4 );
  sparse_grid_composite_test01 ( 100, 100, 0, 1 );
%
%  Compute abstract grid indices of sparse grid points as selected from
%  product grid for DIMENSION, LEVEL_MAX.
%
  dim_num = 2;
  level_max = 3;
  sparse_grid_composite_test02 ( dim_num, level_max );

  dim_num = 2;
  level_max = 4;
  sparse_grid_composite_test02 ( dim_num, level_max );

  dim_num = 3;
  level_max = 0;
  sparse_grid_composite_test02 ( dim_num, level_max );

  dim_num = 3;
  level_max = 2;
  sparse_grid_composite_test02 ( dim_num, level_max );

  dim_num = 6;
  level_max = 2;
  sparse_grid_composite_test02 ( dim_num, level_max );
%
%  Compute sparse composite rule for DIMENSION, LEVEL_MAX.
%
  dim_num = 2;
  level_max = 3;
  sparse_grid_composite_test03 ( dim_num, level_max );

  dim_num = 3;
  level_max = 0;
  sparse_grid_composite_test03 ( dim_num, level_max );

  dim_num = 3;
  level_max = 1;
  sparse_grid_composite_test03 ( dim_num, level_max );
%
%  Test sum of weights for DIMENSION, LEVEL_MAX.
%
  sparse_grid_composite_test04 ( 2, 4 );
  sparse_grid_composite_test04 ( 3, 0 );
  sparse_grid_composite_test04 ( 3, 1 );
  sparse_grid_composite_test04 ( 3, 6 );
  sparse_grid_composite_test04 ( 10, 3 );
%
%  Test monomial exactness for DIMENSION, LEVEL_MAX, DEGREE_MAX.
%
  sparse_grid_composite_test05 ( 2, 0, 2 );
  sparse_grid_composite_test05 ( 2, 1, 2 );
  sparse_grid_composite_test05 ( 2, 2, 2 );
  sparse_grid_composite_test05 ( 2, 3, 2 );
  sparse_grid_composite_test05 ( 2, 4, 2 );
  sparse_grid_composite_test05 ( 2, 5, 2 );

  sparse_grid_composite_test05 ( 3, 0, 2 );
  sparse_grid_composite_test05 ( 3, 1, 2 );
  sparse_grid_composite_test05 ( 3, 2, 2 );
  sparse_grid_composite_test05 ( 3, 3, 2 );
%
%  Show how to write a rule to a file.
%
  dim_num = 2;
  level_max = 3;
  sparse_grid_composite_test06 ( dim_num, level_max );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_COMPOSITE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function Rx=xrot(phi)

Rx = [1 0 0; 0 cos(phi) -sin(phi);0 sin(phi) cos(phi)];

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function hmm = versCompatibilityFix(hmm)

if ~isfield(hmm.train,'orders')
    hmm.train.orders = formorders(hmm.train.order,hmm.train.orderoffset,...
        hmm.train.timelag,hmm.train.exptimelag);
end

if isfield(hmm,'cache')
    hmm = rmfield(hmm,'cache');
end

if isfield(hmm.train,'state')
    hmm.train = rmfield(hmm.train,'state');
end
if ~isfield(hmm.train,'onpower')
    hmm.train.onpower = 0;
end
if ~isfield(hmm.train,'leida')
    hmm.train.leida = 0;
end
if ~isfield(hmm.train,'embeddedlags')
    hmm.train.embeddedlags = 0;
end
if ~isfield(hmm.train,'embeddedlags_batched')
    hmm.train.embeddedlags_batched = 0;
end
if ~isfield(hmm.train,'pca')
    hmm.train.pca = 0;
end
if ~isfield(hmm.train,'filter')
    hmm.train.filter = [];
end
if ~isfield(hmm.train,'detrend')
    hmm.train.detrend = 0;
end
if ~isfield(hmm.train,'downsample')
    hmm.train.downsample = 0;
end
if ~isfield(hmm.train,'leakagecorr')
    hmm.train.leakagecorr = 0;
end
if ~isfield(hmm.train,'Pstructure')
    hmm.train.Pstructure = true(hmm.K);
end
if ~isfield(hmm.train,'Pistructure')
    hmm.train.Pistructure = true(1,hmm.K);
end
if ~isfield(hmm.train,'PriorWeightingP')
    hmm.train.PriorWeightingP = 1;
end
if ~isfield(hmm.train,'PriorWeightingPi')
    hmm.train.PriorWeightingPi = 1;
end
if ~isfield(hmm.train,'lowrank')
    hmm.train.lowrank = 0;
end
if ~isfield(hmm.train,'updateObs')
    hmm.train.updateObs = 1;
end
if ~isfield(hmm.train,'updateGamma')
    hmm.train.updateGamma = 1;
end
if ~isfield(hmm.train,'updateP')
    hmm.train.updateP = 1;
end
if ~isfield(hmm.train,'id_mixture')
    hmm.train.id_mixture = 0;
end
if ~isfield(hmm.train,'id_mixture')
    hmm.train.id_mixture = 0;
end
if ~isfield(hmm.train,'distribution')
    hmm.train.distribution = 'Gaussian';
end
if ~isfield(hmm.train,'pca_spatial')
    hmm.train.pca_spatial = 0;
end
if ~isfield(hmm.train,'Gamma_constraint')
    hmm.train.Gamma_constraint = [];
end
if ~isfield(hmm.train,'Gamma_constraint_radius')
    hmm.train.Gamma_constraint = [];
end
if ~isfield(hmm.train,'acrosstrial_constrained')
    hmm.train.acrosstrial_constrained = [];
end
if ~isfield(hmm.train,'episodic')
    hmm.train.episodic = 0;
end
end

author: OHBA-analysis

repo: HMM-MAR


text:
report_this_filefun(mfilename('fullpath'));

SA = 1;

ste = [bvg(:,2) bvg(:,1)+180  bvg(:,4) bvg(:,3)+180 bvg(:,6) bvg(:,5)+180 bvg(:,7) ];
sor = ste;
ste2 = ste;
sor(:,SA*2) = sor(:,SA*2)+90;

normlap2=ones(length(tmpgri(:,1)),1)*nan;

re3=reshape(normlap2,length(yvect),length(xvect));
s11 = re3;

normlap2(ll)= bvg(:,8);
rama=reshape(normlap2,length(yvect),length(xvect));

normlap2(ll)= bvg(:,7);
r=reshape(normlap2,length(yvect),length(xvect));

normlap2(ll)= bvg(:,1);
s11=reshape(normlap2,length(yvect),length(xvect));

normlap2(ll)= bvg(:,4);
s31=reshape(normlap2,length(yvect),length(xvect));

normlap2(ll)= bvg(:,2);
s12=reshape(normlap2,length(yvect),length(xvect));

old1 = re3;

figure_w_normalized_uicontrolunits('visible','on')
plq = quiver(newgri(:,1),newgri(:,2),-cos(sor(:,1)*pi/180),sin(sor(:,1)*pi/180),0.6,'.')
set(plq,'LineWidth',0.5,'Color','k')
hold on

%pl = quiver(newgri(:,1),newgri(:,2),-cos(s(:,2)*pi/180/3),-sin(s(:,2)*pi/180),2)
%set(pl,'LineWidth',2,'Color','r')

px = get(plq,'Xdata'); py = get(plq,'Ydata');
close
figure_w_normalized_uicontrolunits('pos',[100 100 860 600])
watchon;
%whitebg(gcf);
set(gcf,'color','w');
axes('pos',[0.12 0.12 0.8 0.8]);
hold on
n = 0;
l0 = []; l1 = []; l2 = []; l3 = []; l4 = []; l5 = [];
for i = 1:3:length(px)-1
    n = n+1;j = jet;
    col = floor(ste(n,SA*2-1)/60*62)+1;
    if col > 64 ; col = 64; end
    pl = plot(px(i:i+1),py(i:i+1),'k','Linewidth',1.,'Markersize',1,'color',[ 0 0 0  ] );
    hold on
    dx = px(i)-px(i+1);
    dy = py(i) - py(i+1);
    pl2 = plot(px(i),py(i),'ko','Markersize',0.1,'Linewidth',0.5,'color',[0 0 0] );
    l0 = pl;
    pl3 = plot([px(i) px(i)+dx],[py(i) py(i)+dy],'k','Linewidth',1.,'color',[0 0 0] );

  if ste(n,1) < 40  && ste(n,3) > 45  && ste(n,5) < 20 ; set([pl pl3],'color',[0.2 0.8 0.2]); set(pl2,'color',[0.2 0.8 0.2]); l3 = pl; end
   if ste(n,1) > 52  &&                 ste(n,5) < 35 ;                 set([pl pl3],'color','r'); set(pl2,'color','r'); l1 = pl; end
  if ste(n,1) < 20  && ste(n,3) > 45  && ste(n,5) < 40 ;  set([pl pl3],'color',[0.2 0.8 0.2]); set(pl2,'color',[0.2 0.8 0.2]);l3 = pl; end
  if ste(n,1) > 40  &&                 ste(n,1) < 52  && ste(n,5) < 20 ;  set([pl pl3],'color','m'); set(pl2,'color','m'); l2 = pl; end
  if ste(n,1) < 20  &&                 ste(n,5) > 40  &&  ste(n,5) <  52 ; set([pl pl3],'color','c'); set(pl2,'color','c');l4 = pl; end
   if ste(n,1) < 35  && ste(n,5) > 52  ; set([pl pl3],'color','b'); set(pl2,'color','b');l5 = pl;  end

end

if isempty(l1) == 1; pl2 = plot(px(i),py(i),'kx','Linewidth',1.,'color','r'); l1 = pl2; set(l1,'visible','off'); end
if isempty(l2) == 1; pl2 = plot(px(i),py(i),'kx','Linewidth',1.,'color','m'); l2 = pl2; set(l2,'visible','off'); end
if isempty(l3) == 1; pl2 = plot(px(i),py(i),'kx','Linewidth',1.,'color',[0.2 0.8 0.2] ); l3 = pl2; set(l3,'visible','off'); end
if isempty(l4) == 1; pl2 = plot(px(i),py(i),'kx','Linewidth',1.,'color','c' ); l4 = pl2; set(l4,'visible','off'); end
if isempty(l5) == 1; pl2 = plot(px(i),py(i),'kx','Linewidth',1.,'color','b' ); l5 = pl2; set(l5,'visible','off'); end
if isempty(l0) == 1; l0 = plot(px(i),py(i),'kx','Linewidth',1.,'color',[0 0 0 ] );  set(l0,'visible','off'); end

legend([l1 l2 l3 l4 l5 l0],'NF','NS','SS','TS','TF','U');

hold on
axis('equal')
axis('tight')

title2([name ';  '   num2str(t0b) ' to ' num2str(teb) ],'FontSize',ZmapGlobal.Data.fontsz.s,...
    'Color','k','FontWeight','normal')

xlabel('Distance along projection [km] ','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
ylabel('Depth [km] ','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)

set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','normal',...
    'FontWeight','normal','LineWidth',1.,...
    'Box','on','TickDir','out','Ticklength',[0.01 0.01])

matdraw

watchoff;

% view the variance map
re3 = r;sha = 'in';
%view_varmap

author: CelsoReyes