text:
function [anybad,wh] = iimg_check_volinfo(maskInfo,imgInfo)
% Checks a series of image .mat files and dims against a reference
% (maskInfo)
%
% :Usage:
% ::
%
%     anybad = iimg_check_volinfo(maskInfo,imgInfo)
%
% :Inputs:
%
%   **maskInfo** and **volInfo :**
%        are spm-style volume info structures
%
% see spm_vol.m

wh = [];

n = length(imgInfo);
notok = zeros(1,n);

tol = .01;

for i=1:n
    chk = abs(maskInfo.mat - imgInfo(i).mat) > tol;

    chk = any(diag(chk(1:3, 1:3)));

% %     chk = chk(:);
% %     chk = chk(1:end-1);             % eliminate SPM scale factor and translation
% %     chk1 = any(chk);

    chk2 = any(maskInfo.dim(1:3) - imgInfo(i).dim(1:3));

    notok(i) = chk | chk2;
end

anybad = any(notok);

if anybad
    wh = find(notok);

    disp('The following images'' mat files or dims differed from the first:')
    disp(num2str(wh));

    disp('First mat:');
    disp(maskInfo.mat);
    disp(maskInfo.dim);

    disp('Bad mats:');
    for i = wh
        disp(imgInfo(i).fname);
        disp(imgInfo(i).mat);
        disp(imgInfo(i).dim);
    end
end

return


text:
%ISCOLOR Test for color image
%
% ISCOLOR(IM) is true (1) if IM is a color image, that is, it its third
% dimension is equal to three.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function s = iscolor(im)
    % WxH is mono
    % WxHx3 is color
    % WxHxN is mono sequence
    % WxHx3xN is color sequence
    s = isnumeric(im) && size(im,1) > 1 && size(im,2) > 1 && size(im,3) == 3;

author: petercorke

repo: machinevision-toolbox-matlab


text:
function n = nrow(x)

%	function n = nrow(x)
%	Return number of rows in x.

	n = size(x); n = n(1);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function varargout = process_evt_groupname( varargin )
% PROCESS_EVT_GROUPNAME: Combine different categories of events into one (by name)
%
% USAGE:  OutputFiles = process_evt_groupname('Run', sProcess, sInputs)

% @=============================================================================
% 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, 2013-2023

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Group by name';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Events';
    sProcess.Index       = 51;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/EventMarkers#Other_menus';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data', 'raw', 'matrix'};
    sProcess.OutputTypes = {'data', 'raw', 'matrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Event name
    sProcess.options.combine.Comment = ['Example: We have three events (A,B,C) and want to create new combinations:
' ...
                                        '    E: Event A and B occurring at the same time
' ...
                                        '    F: Event A and C occurring at the same time
' ...
                                        'For that, use the following classification:
' ...
                                        '    E = A,B
' ...
                                        '    F = A,C
' ...
                                        'You may add as many combinations as needed, one per line.
' ...
                                        'You can rename events with the following syntax: E=A.
' ...
                                        'You can delete or keep the original events (A,B,C) with the checkbox below.

'];
    sProcess.options.combine.Type    = 'textarea';
    sProcess.options.combine.Value   = '';
    % Maximum duration between simulateous events
    sProcess.options.dt.Comment = 'Maximum delay between simultaneous events: ';
    sProcess.options.dt.Type    = 'value';
    sProcess.options.dt.Value   = {0, 'ms', 0};
    % Maximum duration between simulateous events
    sProcess.options.order.Comment = {'First', 'Last', 'Event to select (if not strictly aligned): '; ...
                                     'first', 'last', ''};
    sProcess.options.order.Type    = 'radio_linelabel';
    sProcess.options.order.Value   = 'first';
    % Delete original events
    sProcess.options.delete.Comment = 'Delete the original events';
    sProcess.options.delete.Type    = 'checkbox';
    sProcess.options.delete.Value   = 0;
end

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

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Return all the input files
    OutputFiles = {};

    % ===== GET OPTIONS =====
    % Get the options
    isDelete = sProcess.options.delete.Value;
    Order = sProcess.options.order.Value;
    dt = sProcess.options.dt.Value{1};
    % Combination string
    combineStr = strtrim(sProcess.options.combine.Value);
    % Split in lines
    combine_lines = str_split(combineStr, [10 13]);
    % Split each line
    combineCell = {};
    for i = 1:length(combine_lines)
        % No information on this line: skip
        combine_lines{i} = strtrim(combine_lines{i});
        if isempty(combine_lines{i})
            continue;
        end
        % Split with "="
        lineCell = str_split(combine_lines{i}, '=');
        if (length(lineCell) ~= 2)
            continue;
        end
        % Split with ";,"
        eventsCell = str_split(lineCell{2}, ';,');
        % Add combination entry
        iComb = size(combineCell,1) + 1;
        combineCell{iComb,1} = strtrim(lineCell{1});
        combineCell{iComb,2} = cellfun(@strtrim, eventsCell, 'UniformOutput', 0);
    end
    % If no combination available
    if isempty(combineCell)
        bst_report('Error', sProcess, [], 'Invalid combinations format.');
        return;
    end

    % ===== PROCESS ALL FILES =====
    % For each file
    for iFile = 1:length(sInputs)
        % ===== GET FILE DESCRIPTOR =====
        % Load the raw file descriptor
        isRaw = strcmpi(sInputs(iFile).FileType, 'raw');
        if isRaw
            DataMat = in_bst_data(sInputs(iFile).FileName, 'F');
            sEvents = DataMat.F.events;
            sFreq = DataMat.F.prop.sfreq;
        else
            DataMat = in_bst_data(sInputs(iFile).FileName, 'Events', 'Time');
            sEvents = DataMat.Events;
            sFreq = 1 ./ (DataMat.Time(2) - DataMat.Time(1));
        end
        % If no markers are present in this file
        if isempty(sEvents)
            bst_report('Error', sProcess, sInputs(iFile), 'This file does not contain any event. Skipping File...');
            continue;
        end
        % Convert the distance in time to distance in samples
        ds = round(dt .* sFreq);
        % Call the grouping function
        [sEvents, isModified] = Compute(sInputs(iFile), sEvents, combineCell, ds, isDelete, sFreq, Order);

        % ===== SAVE RESULT =====
        % Only save changes if something was change
        if isModified
            % Report changes in .mat structure
            if isRaw
                DataMat.F.events = sEvents;
            else
                DataMat.Events = sEvents;
            end
            % Save file definition
            bst_save(file_fullpath(sInputs(iFile).FileName), DataMat, 'v6', 1);
        end
        % Return all the input files
        OutputFiles{end+1} = sInputs(iFile).FileName;
    end
end

%% ===== GROUP EVENTS =====
function [eventsNew, isModified] = Compute(sInput, events, combineCell, ds, isDelete, sfreq, Order)
    % Parse inputs
    if (nargin < 7) || isempty(Order)
        Order = 'first';
    end
    % No modification
    isModified = 0;
    eventsNew = events;
    removeEvt = zeros(0,2);
    % Loop on the different combinations
    for iComb = 1:size(combineCell,1)
        AllEvt = zeros(2,0);
        iEvtList = [];
        isSkip = 0;
        % Use the output of the previous combination as the input to this combination
        events = eventsNew;
        % Remove duplicated events
        evtList = unique(combineCell{iComb,2});
        % Get events for this combination
        for iCombEvt = 1:length(evtList)
            % Find event in the list
            evtLabel = evtList{iCombEvt};
            iEvt = find(strcmpi({events.label}, evtLabel));
            % If events doesn't exist: skip
            if isempty(iEvt)
                bst_report('Warning', 'process_evt_groupname', sInput, ['Event "' evtLabel '" does not exist. Skipping group...']);
                isSkip = 1;
                break;
            end
            % If renaming event: simply update the label and move on
            if (length(combineCell{iComb,2}) == 1)
                eventsNew(iEvt).label = combineCell{iComb,1};
                isSkip = 1;
                isModified = 1;
                break;
            end
            % If events are extended events: skip
            if (size(events(iEvt).times,1) > 1)
                bst_report('Error', 'process_evt_groupname', sInput, 'Cannot process extended events. Skipping group...');
                isSkip = 1;
                break;
            end
            % Add to the list of all the processes
            iEvtList(end+1) = iEvt;
            AllEvt = [AllEvt, [round(events(iEvt).times .* sfreq); repmat(iEvt, size(events(iEvt).times))]];
        end
        % Skip combination if one of the events is not found or not a simple event
        if isSkip
            continue;
        end

        % Fix the events times according to the maximum allowed distance between events
        if (ds > 0)
            N = size(AllEvt,2);
            % Calculate distance between all the pairs of events:
            % [alldist(i,j)=5] means that event #i is 5 samples before event #j
            alldist = repmat(AllEvt(1,:),N,1) - repmat(AllEvt(1,:)',1,N);
            % Get the distances for the different event types only: we don't want to collapse two events of the same category
            if (length(iEvtList) > 1)
                diffmask = (repmat(AllEvt(2,:),N,1) ~= repmat(AllEvt(2,:)',1,N));
            % Unless there is only one event in input
            else
                diffmask = (eye(N) == 0);
            end
            % Find the events that can be collapsed
            collapse = ((alldist > 0) & (alldist <= ds)) .* diffmask;
            collapse(:,sum(collapse,1) > 1) = 0;
            collapse(sum(collapse,2) > 1,:) = 0;
            [iBefore, iAfter] = find(collapse);
            % Remove cycles
            iCycles = find(ismember(iBefore, iAfter));
            iBefore(iCycles) = [];
            iAfter(iCycles) = [];
            % Replace the time of the last events with the time of the first event
            if ~isempty(iBefore)
                for iEvt = 1:length(iAfter)
                    smpBefore = round(events(AllEvt(2,iBefore(iEvt))).times .* sfreq);
                    smpAfter = round(events(AllEvt(2,iAfter(iEvt))).times .* sfreq);
                    iOccBefore = find(smpBefore == AllEvt(1,iBefore(iEvt)));
                    iOccAfter  = find(smpAfter == AllEvt(1,iAfter(iEvt)));
                    switch (Order)
                        case 'first'
                            events(AllEvt(2,iAfter(iEvt))).times(iOccAfter) = events(AllEvt(2,iBefore(iEvt))).times(iOccBefore);
                        case 'last'
                            events(AllEvt(2,iAfter(iEvt))).times(iOccBefore) = events(AllEvt(2,iBefore(iEvt))).times(iOccAfter);
                    end
                end
                switch (Order)
                    case 'first'
                        AllEvt(1,iAfter) = AllEvt(1,iBefore);
                    case 'last'
                        AllEvt(1,iBefore) = AllEvt(1,iAfter);
                end
            end
        end

        % Process each unique time value
        uniqueSamples = unique(AllEvt(1,:));
        for iSmp = 1:length(uniqueSamples)
            % Look for all the events happening at this time
            iEvts = AllEvt(2, (AllEvt(1,:) == uniqueSamples(iSmp)));
            % If only one occurrence: skip to the next time
            if (length(iEvtList) > 1) && (length(iEvts) < length(iEvtList))   % Two different events
                continue;
            elseif (length(iEvtList) == 1) && (length(iEvts) < 2)    % One event with itself
                continue;
            end
            % Remove occurrence from each event type (and build new event name)
            for i = 1:length(iEvts)
                % Find the occurrence indice
                iOcc = find(round(events(iEvts(i)).times .* sfreq) == uniqueSamples(iSmp));
                % Get the values
                if (i == 1)
                    newTime     = events(iEvts(i)).times(iOcc);
                    newEpoch    = events(iEvts(i)).epochs(iOcc);
                    if ~isempty(events(iEvts(i)).channels)
                        newChannels = events(iEvts(i)).channels(iOcc);
                    else
                        newChannels = [];
                    end
                    if ~isempty(events(iEvts(i)).notes)
                        newNotes = events(iEvts(i)).notes(iOcc);
                    else
                        newNotes = [];
                    end
                end
                % Remove this occurrence
                if isDelete
                    removeEvt = [removeEvt; repmat(iEvts(i),length(iOcc),1), reshape(iOcc,[],1)];
                end
            end
            % New event name
            newLabel = combineCell{iComb,1};
            % Find this event in the list
            iNewEvt = find(strcmpi(newLabel, {eventsNew.label}));
            % Create event category if does not exist yet
            if isempty(iNewEvt)
                % Initialize new event
                iNewEvt = length(eventsNew) + 1;
                sEvent = db_template('event');
                sEvent.label = newLabel;
                % Re-use color of renamed event
                if (length(iEvtList) == 1) && isDelete
                    sEvent.color = events(iEvtList(1)).color;
                % Create a new color
                else
                    sEvent.color = panel_record('GetNewEventColor', iNewEvt, eventsNew);
                end
                % Add new event to list
                eventsNew(iNewEvt) = sEvent;
            end
            % Add occurrences
            eventsNew(iNewEvt).times    = [eventsNew(iNewEvt).times,    newTime];
            eventsNew(iNewEvt).epochs   = [eventsNew(iNewEvt).epochs,   newEpoch];
            % Sort
            [eventsNew(iNewEvt).times, indSort] = unique(eventsNew(iNewEvt).times);
            eventsNew(iNewEvt).epochs   = eventsNew(iNewEvt).epochs(indSort);
            % Channels
            if ~isempty(newChannels)
                if isempty(eventsNew(iNewEvt).channels)
                    eventsNew(iNewEvt).channels = cell(1, size(eventsNew(iNewEvt).times,2) - size(newTime,2));
                end
                eventsNew(iNewEvt).channels = [eventsNew(iNewEvt).channels, newChannels];
                eventsNew(iNewEvt).channels = eventsNew(iNewEvt).channels(indSort);
            end
            % Notes
            if ~isempty(newNotes)
                if isempty(eventsNew(iNewEvt).notes)
                    eventsNew(iNewEvt).notes = cell(1, size(eventsNew(iNewEvt).times,2) - size(newTime,2));
                end
                eventsNew(iNewEvt).notes = [eventsNew(iNewEvt).notes, newNotes];
                eventsNew(iNewEvt).notes = eventsNew(iNewEvt).notes(indSort);
            end
            isModified = 1;
        end
    end
    % Remove events that were tagged
    if ~isempty(removeEvt)
        uniqueEvt = unique(removeEvt(:,1));
        for i = 1:length(uniqueEvt)
            iEvt = uniqueEvt(i);
            iOcc = removeEvt(removeEvt(:,1)==iEvt,2)';
            eventsNew(iEvt).times(iOcc)    = [];
            eventsNew(iEvt).epochs(iOcc)   = [];
            if ~isempty(eventsNew(iEvt).channels)
                eventsNew(iEvt).channels(iOcc) = [];
            end
            if ~isempty(eventsNew(iEvt).notes)
                eventsNew(iEvt).notes(iOcc) = [];
            end
        end
        % Remove empty categories
        iEmptyEvt = find(cellfun(@isempty, {eventsNew.times}));
        if ~isempty(iEmptyEvt)
            eventsNew(iEmptyEvt) = [];
        end
    end
end

author: brainstorm-tools


text:
%% ----------------------------
% Input: Work_mode: Mode of working condition 1 --> BBDST, 2 --> constant current
%                   SOC_est_init: The initial value of estimated SOC
%% ----------------------------
function main(Work_mode, SoC_est_init)
    if nargin == 0  % Set parameter by default
        Work_mode = 1;
        SoC_est_init = 1;
    elseif nargin == 1
        SoC_est_init = 1;
    end
    if Work_mode == 1
        sim BBDST_workingcondition;
        I = -(current.data)' * 1.5 / 50;
    elseif Work_mode == 2
        N = 60001;
        I = 1.5 * ones(1, N);
        I(ceil(N / 5) : ceil(N * 3 / 9)) = 0;
        I(ceil(N * 5 / 9) : ceil(N * 4 / 5)) = 0;
    else
        disp("Input error!");
        disp("Work_mode: Mode of working condition");
        disp("           1 --> BBDST, 2 --> constant current ");
        disp("SOC_est_init : The initial value of estimated SOC");
        return;
    end
    tic;  % start time
    [avr_err_EKF, std_err_EKF, avr_err_UKF, std_err_UKF] = EKF_UKF_Thev(SoC_est_init, I);
    toc;  % end time
    fprintf('Initial SOC : %f\nWorking Mode: %d\n', SoC_est_init, Work_mode);
    fprintf("avr_err_EKF --> %f\n", avr_err_EKF);
    fprintf("standard_err_EKF --> %f\n", std_err_EKF);
    fprintf("avr_err_UKF --> %f\n", avr_err_UKF);
    fprintf("standard_err_UKF --> %f\n", std_err_UKF);
end

author: AlterWL


text:
function results = test_sscor(eeg, model, is_ensemble)
% Test phase of the sum of squared correlations (SSCOR)-based
% steady-state visual evoked potentials (SSVEPs) detection [1].
% * This script is the exact same as test_trca.m except for the function
% name.
%
% function results = test_sscor(eeg, model, is_ensemble)
%
% Input:
%   eeg             : Input eeg data
%                     (# of targets, # of channels, Data length [sample])
%   model           : Learning model for tesing phase of the ensemble
%                     SSCOR-based method
%   is_ensemble     : 0 -> SSCOR-based method,
%                     1 -> Ensemble SSCOR-based method (defult: 1)
%
% Output:
%   results         : The target estimated by this method
%
% See also:
%   train_sscor.m
%
% Reference:
%   [1] G. R. Kumar and M. R. Reddy,
%       "Designing a Sum of Squared Correlations Framework for Enhancing SSVEP
%        Based BCIs",
%       IEEE Trans. Neural Syst. Rehabil. Eng., vol. 27, pp. 2044-2050, 2019.
%
% Kuan-Jung Chiang and Masaki Nakanishi, 25-Nov-2019
% Swartz Center for Computational Neuroscience, Institute for Neural
% Computation, University of California San Diego
% E-mail: masaki@sccn.ucsd.edu

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

if ~exist('model', 'var')
    error('Training model based on SSCOR is required. See train_sscor().');
end

fb_coefs = [1:model.num_fbs].^(-1.25)+0.25;

for targ_i = 1:1:model.num_targs
    test_tmp = squeeze(eeg(targ_i, :, :));
    for fb_i = 1:1:model.num_fbs
        testdata = filterbank(test_tmp, model.fs, fb_i);
        for class_i = 1:1:model.num_targs
            traindata =  squeeze(model.trains(class_i, fb_i, :, :));
            if ~is_ensemble
                w = squeeze(model.W(fb_i, class_i, :));
            else
                w = squeeze(model.W(fb_i, :, :))';
            end
            r_tmp = corrcoef(testdata'*w, traindata'*w);
            r(fb_i,class_i) = r_tmp(1,2);
        end % class_i
    end % fb_i
    rho = fb_coefs*r;
    [~, tau] = max(rho);
    results(targ_i) = tau;
end % targ_i

author: mnakanishi

repo: TRCA-SSVEP


text:
function eframes = vl_frame2oell(frames)
% FRAMES2OELL   Convert generic feature frames to oriented ellipses
%   EFRAMES = VL_FRAME2OELL(FRAMES) converts the specified FRAMES to
%   the oriented ellipses EFRAMES.
%
%   A frame is either a point, disc, oriented disc, ellipse, or
%   oriented ellipse. These are represened respecively by
%   2, 3, 4, 5 and 6 parameters each, as described in VL_PLOTFRAME().
%
%   An oriented ellipse is the most general frame. When an unoriented
%   frame is converted to an oriented ellipse, the rotation is selected
%   so that the positive Y direction is unchanged.
%
%   See: VL_PLOTFRAME(), VL_HELP().

% Author: Andrea Vedaldi

% Copyright (C) 2013 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

[D,K] = size(frames) ;
eframes = zeros(6,K) ;

switch D
  case 2
    eframes(1:2,:) = frames(1:2,:) ;

  case 3
    eframes(1:2,:) = frames(1:2,:) ;
    eframes(3,:)   = frames(3,:) ;
    eframes(6,:)   = frames(3,:) ;

  case 4
    r = frames(3,:) ;
    c = r.*cos(frames(4,:)) ;
    s = r.*sin(frames(4,:)) ;

    eframes(1:2,:) = frames(1:2,:) ;
    eframes(3:6,:) = [c ; s ; -s ; c] ;

  case 5
    eframes(1:2,:) = frames(1:2,:) ;
    eframes(3:6,:) = mapFromS(frames(3:5,:)) ;

  case 6
    eframes = frames ;

  otherwise
     error('FRAMES format is unknown.') ;
end

% --------------------------------------------------------------------
function A = mapFromS(S)
% --------------------------------------------------------------------
% Returns the (stacking of the) 2x2 matrix A that maps the unit circle
% into the ellipses satisfying the equation x' inv(S) x = 1. Here S
% is a stacked covariance matrix, with elements S11, S12 and S22.
%
% The goal is to find A such that AA' = S. In order to let the Y
% direction unaffected (upright feature), the assumption is taht
% A = [a b ; 0 c]. Hence
%
%  AA' = [a^2, ab ; ab, b^2+c^2] = S.

A = zeros(4,size(S,2)) ;
a = sqrt(S(1,:));
b = S(2,:) ./ max(a, 1e-18) ;

A(1,:) = a ;
A(2,:) = b ;
A(4,:) = sqrt(max(S(3,:) - b.*b, 0)) ;

author: jianxiongxiao

repo: ProfXkit


text:
function laplacian_test01 ( )

%*****************************************************************************80
%
%% LAPLACIAN_TEST01 tests L1DD and similar routines.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAPLACIAN_TEST01\n' );
  fprintf ( 1, '  A full-storage matrix is returned by:\n' );
  fprintf ( 1, '  L1DD: Dirichlet/Dirichlet BC;\n' );
  fprintf ( 1, '  L1DN: Dirichlet/Neumann BC;\n' );
  fprintf ( 1, '  L1ND: Neumann/Dirichlet BC;\n' );
  fprintf ( 1, '  L1NN: Neumann/Neumann BC;\n' );
  fprintf ( 1, '  L1PP: Periodic BC;\n' );

  n = 5;

  for test = 1 : 2

    if ( test == 1 )
      h = 1.0;
    else
      h = 1.0 / ( n + 1 );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using spacing H = %g\n', h );

    l = l1dd ( n, h );
    r8mat_print ( n, n, l, '  L1DD:' );

    l = l1dn ( n, h );
    r8mat_print ( n, n, l, '  L1DN:' );

    l = l1nd ( n, h );
    r8mat_print ( n, n, l, '  L1ND:' );

    l = l1nn ( n, h );
    r8mat_print ( n, n, l, '  L1NN:' );

    l = l1pp ( n, h );
    r8mat_print ( n, n, l, '  L1PP:' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function volcrop = crop_volume(vol, crp)
if ndims(vol) == 3
    volcrop = vol(crp(1,1):crp(1,2), crp(2,1):crp(2,2), crp(3,1):crp(3,2));
elseif ndims(vol) == 2
    volcrop = vol(crp(1,1):crp(1,2), crp(2,1):crp(2,2));
end
end

repo: pTVreg


text:
function set_rotation_pointer(h)
    %set rotation pointer rotate
    zz=load('rotate.mat');
    val = sum(zz.cdata,3)+1;
    h.Pointer='custom';
    h.PointerShapeCData=[nan(1,16);val];
end

author: CelsoReyes


text:
function simPrices = simulatePrices(expReturns, expCovariances, ...
    price0, allTimes)
% SIMULATEPRICES: calculate Monte Carlo prices for a set of correlated
% assets
%
% INPUTS:
% expReturn: a 1-by-nAssets vector of expected returns
% expCovariances: a nAssets-by-nAssets matrix of expected covariances
% price0: a 1-by-nAssets of initial asset prices
% allTimes: an nTimes-by-1 vector of integer times for which we need to
%   calculate the prices (ex. [30 90 180]' for 1-, 3-, and 6-month prices)
% nSims: a scalar integer for the number of Monte Carlo paths generated
%
% OUTPUT:
% simPrices: an nTimes-by-nAssets matrix of the simulated prices at the
% requested times.
%
% The variables expReturns, expCovariances, and allTimes share the same
% units-- if allTimes is in units of days, then the returns and covariances
% must be _daily_ moments.

% Calculates correlated, multivariate, normally-distributed numbers.  Uses
% Statistics Toolbox.
simReturns = mvnrnd(expReturns, expCovariances, max(allTimes));

% Convert simple returns to prices
simPrices = cumprod([price0; simReturns+1]);

% Downsample to only the times requested.  The first row is t==0 (i.e.:
% price0)
simPrices = simPrices(allTimes+1, :);

author: Sable

repo: mcbench-benchmarks


text:
function chanind = selectchannels(this, channels)
% Method for getting channel indices based on labels and/or types
% FORMAT  res = selectchannels(this, label)
% this       - MEEG object
% channels   - string or cell array of labels that may also include
%              'all', or types ('EEG', 'MEG' etc.)
%
% res        - vector of channel indices matching labels
%__________________________________________________________________________
% Copyright (C) 2010-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: selectchannels.m 7253 2018-02-04 17:20:57Z vladimir $

if ischar(channels)
    channels = {channels};
end

chanind = [];

for i = 1:numel(channels)
    if strncmpi('regexp_', channels{i}, 7)
        re        = channels{i}(8:end);
        match     = regexp(chanlabels(this), re);
        chanind   = [chanind find(~cellfun('isempty', match))];
    else
        cind    = indchannel(this, channels{i});
        if ~isempty(cind)
            chanind = [chanind cind];
        elseif ismember(upper(channels{i}), ...
                {'ALL','MEG', 'MEGPLANAR', 'MEGMAG', 'MEGGRAD', 'MEGCOMB','EEG',...
                'EOG', 'ECG', 'EMG', 'LFP', 'SRC', 'PHYS', 'ILAM', 'OTHER', 'REF', 'REFMAG', 'REFGRAD'})
            chanind = [chanind indchantype(this, upper(channels{i}))];
        end
    end

    if any(size(chanind) == 0)
        chanind = [];
    end
end

chanind = unique(chanind);

author: spm


text:
% --------------------------------------------------------------------------------------------------------
function [im_patch, im_patch_original] = get_subwindow_tracking(im, pos, model_sz, original_sz, avg_chans)
%GET_SUBWINDOW_TRACKING Obtain image sub-window, padding with avg channel if area goes outside of border
% -------------------------------------------------------------------------------------------------
    if isempty(original_sz)
        original_sz = model_sz;
    end
    sz = original_sz;
    im_sz = size(im);
    %make sure the size is not too small
    assert(all(im_sz(1:2) > 2));
    c = (sz+1) / 2;

    %check out-of-bounds coordinates, and set them to black
    context_xmin = round(pos(2) - c(2)); % floor(pos(2) - sz(2)/2);
    context_xmax = context_xmin + sz(2) - 1;
    context_ymin = round(pos(1) - c(1)); % floor(pos(1) - sz(1)/2);
    context_ymax = context_ymin + sz(1) - 1;
    left_pad = max(0, 1-context_xmin);
    top_pad = max(0, 1-context_ymin);
    right_pad = max(0, context_xmax - im_sz(2));
    bottom_pad = max(0, context_ymax - im_sz(1));

    context_xmin = context_xmin + left_pad;
    context_xmax = context_xmax + left_pad;
    context_ymin = context_ymin + top_pad;
    context_ymax = context_ymax + top_pad;

    if top_pad || left_pad
        R = padarray(im(:,:,1), [top_pad left_pad], avg_chans(1), 'pre');
        G = padarray(im(:,:,2), [top_pad left_pad], avg_chans(2), 'pre');
        B = padarray(im(:,:,3), [top_pad left_pad], avg_chans(3), 'pre');
        im = cat(3, R, G, B);
    end

    if bottom_pad || right_pad
        R = padarray(im(:,:,1), [bottom_pad right_pad], avg_chans(1), 'post');
        G = padarray(im(:,:,2), [bottom_pad right_pad], avg_chans(2), 'post');
        B = padarray(im(:,:,3), [bottom_pad right_pad], avg_chans(3), 'post');
        im = cat(3, R, G, B);
    end

    xs = context_xmin : context_xmax;
    ys = context_ymin : context_ymax;

    im_patch_original = im(ys, xs, :);
    if ~isequal(model_sz, original_sz)
        im_patch = imresize(im_patch_original, model_sz(1)/original_sz(1));
        % Strangely, sometimes model_sz/original_sz doesn't output model_sz and output of rescaling has to be forced to model_sz
        % (gpu version of imresize does not allow array as model size, only ratio
        if size(im_patch,1)~=model_sz(1)
            % WARNING: camera ready used mexResize and not imresize, this could cause minor differences in results
%             im_patch = gpuArray(mexResize(gather(im_patch_original), model_sz, 'auto'));
            im_patch = gpuArray(imresize(gather(im_patch_original), model_sz));
        end
    else
        im_patch = im_patch_original;
    end
end

author: bertinetto

repo: siamese-fc


text:
% Displays a matrix in a pop-up window for convenient editing.
% This control presents an m-by-n table of matrix entries.

% Copyright 2008-2010 Levente Hunyadi
classdef MatrixEditor < UIControl
    properties
        Type = PropertyType.empty(1,0);
    end
    properties (Dependent)
        Control;
        Item;
        % Whether to show transpose of the matrix instead of actual matrix.
        Transpose;
        % Whether matrix is read-only.
        ReadOnly;
    end
    properties (Access = protected)
        MatrixTranspose = false;
        Matrix;
        Table;
        JideTableModel;
        JideTable;
    end
    properties (Access = private)
        Menu;
    end
    methods
        function self = MatrixEditor(varargin)
            self = self@UIControl(varargin{:});
        end

        function self = Instantiate(self, parent)
            if nargin < 2
                parent = figure;
            end

            % container panel
            panel = uipanel(parent, ...
                'Units', 'normalized', ...
                'Position', [0 0 1 1], ...
                'Tag', '__MatrixEditor__', ...
                'UserData', self);

            % initialization
            com.mathworks.mwswing.MJUtilities.initJIDE;
            com.jidesoft.grid.CellRendererManager.initDefaultRenderer();
            com.jidesoft.grid.CellEditorManager.initDefaultEditor();

            % table model
            self.JideTableModel = handle(javax.swing.table.DefaultTableModel(), 'CallbackProperties');
            set(self.JideTableModel, 'TableChangedCallback', @MatrixEditor.OnTableChanged);

            % JIDE table grid control
            self.JideTable = objectEDT('com.jidesoft.grid.ContextSensitiveTable', self.JideTableModel);
            self.JideTable.setCellSelectionEnabled(true);

            pixelpos = getpixelposition(panel);
            [control,self.Table] = javacomponent(self.JideTable, [10 10 pixelpos(3)-20 pixelpos(4)-20], panel); %#ok<ASGLU>
            set(self.Table, ...
                'Units', 'normalized');

            % menu
            menu = uimenu(ancestor(parent, 'figure'), ...
                'Label', 'Matrix', ...
                'Callback', @self.OnMatrixMenuClicked);
            self.Menu = struct( ...
                'ShowTranspose', uimenu(menu, 'Label', 'Show transpose', 'Callback', @self.OnShowTranspose), ...
                'AddColumnBefore', uimenu(menu, 'Label', 'Add column before', 'Callback', @self.OnAddColumnBefore), ...
                'AddColumnAfter', uimenu(menu, 'Label', 'Add column after', 'Callback', @self.OnAddColumnAfter), ...
                'AddRowBefore', uimenu(menu, 'Label', 'Add row before', 'Callback', @self.OnAddRowBefore), ...
                'AddRowAfter', uimenu(menu, 'Label', 'Add row after', 'Callback', @self.OnAddRowAfter), ...
                'RemoveColumn', uimenu(menu, 'Label', 'Remove column', 'Callback', @self.OnRemoveColumn), ...
                'RemoveRow', uimenu(menu, 'Label', 'Remove row', 'Callback', @self.OnRemoveRow));
            self.OnMatrixMenuClicked();
        end

        function control = get.Control(self)
            control = self.Table;
        end

        function t = get.Transpose(self)
            t = self.MatrixTranspose;
        end

        function self = set.Transpose(self, t)
            validateattributes(t, {'logical'}, {'scalar'});
            item = self.Item;
            self.MatrixTranspose = t;
            self.Item = item;
        end

        function r = get.ReadOnly(self)
            r = ~self.JideTable.isEnabled();
        end

        function self = set.ReadOnly(self, r)
            validateattributes(r, {'logical'}, {'scalar'});
            self.JideTable.setEnabled(~r);
        end

        function matrix = get.Item(self)
            matrix = self.Matrix;
            if self.MatrixTranspose
                matrix = transpose(matrix);
            end
        end

        function self = set.Item(self, matrix)
            validateattributes(matrix, {'numeric','logical'}, {'2d'});
            if self.MatrixTranspose
                matrix = transpose(matrix);
            end
            self.Matrix = matrix;
            if isempty(self.Type)
                self.Type = PropertyType.AutoDiscover(matrix);
            end
            data = self.Type.GetJavaVectorOfVectors(matrix);
            if size(matrix,2) > 0
                list = java.util.Arrays.asList(javaArray('java.lang.String', size(matrix,2)));  % no labels
                labels = java.util.Vector(list);
            else
                labels = java.util.Vector();
            end
            self.JideTableModel.setDataVector(data, labels);
            self.JideTable.doLayout();
        end

        function self = set.Type(self, type)
            validateattributes(type, {'PropertyType'}, {'scalar'});
            self.Type = type;
        end

        function item = GetItemAsString(self)
            validateattributes(self, {'MatrixEditor'}, {'scalar'});
            matrix = self.Matrix;
            if self.MatrixTranspose
                matrix = transpose(matrix);
            end
            item = mat2str(matrix);
        end
    end
    methods (Static, Access = private)
        function self = FindEditor(obj)
        % Finds the editor object to which a table model object belongs.
        %
        % Input arguments:
        % obj:
        %    a javax.swing.table.DefaultTableModel instance
            h = findobjuser(@(userdata) userdata.JideTableModel == obj, '__MatrixEditor__');  % find which HGO contains the object for which the callback is executing
            self = get(h, 'UserData');
        end
    end
    methods (Static)
        function OnTableChanged(source, event)
        % Fired when a matrix entry has been changed.
            javarow = event.getFirstRow();
            javacol = event.getColumn();
            if javarow < 0 || javacol < 0 || event.getType() ~= javax.swing.event.TableModelEvent.UPDATE
                return;  % no selection
            end
            self = MatrixEditor.FindEditor(source);
            row = javarow + 1;
            col = javacol + 1;

            entryvalue = source.getValueAt(javarow, javacol);
            value = entryvalue;  % represented as-is by default
            if ischar(entryvalue)  % represented as text
                if isnumeric(self.Type)
                    value = str2double(entryvalue);
                    if isnan(value)
                        value = [];
                    end
                elseif islogical(self.Type)
                    switch entryvalue
                        case 'true'
                            value = true;
                        case 'false'
                            value = false;
                        otherwise
                            value = str2double(entryvalue);
                            if isnan(value)
                                value = [];
                            end
                    end
                end
            end
            if ~isempty(value)
                try
                    value = self.Type.ConvertFromMatLab(value);
                    self.Matrix(row,col) = value;  % persist changes
                catch %#ok<CTCH>
                    % no action
                end
            end
            if ischar(entryvalue)
                value = self.Matrix(row,col);
                javavalue = self.Type.GetPrimitiveJavaValue(value);
                if javavalue ~= source.getValueAt(javarow, javacol)  % do not trigger change event if value to set is equivalent to value in cell
                    source.setValueAt(javavalue, javarow, javacol);
                end
            end
        end
    end
    methods (Access = private)
        function self = OnMatrixMenuClicked(self, source, event) %#ok<INUSD>
        % Fired when the Matrix menu (or the context menu) is to be shown.
            if ~isempty(self.Type) && ~strcmp(self.Type.Shape, 'column')
                addcolstate = 'on';
            else
                addcolstate = 'off';
            end
            if ~isempty(self.Type) && ~strcmp(self.Type.Shape, 'row')
                addrowstate = 'on';
            else
                addrowstate = 'off';
            end
            if ~isempty(self.Type) && ~strcmp(self.Type.Shape, 'column') ...  % cannot remove single column for column vector
                    && size(self.Matrix,2) > 0 ...     % check if there is anything to remove
                    && ~isempty(self.GetSelectedColumn())
                remcolstate = 'on';
            else
                remcolstate = 'off';
            end
            if ~isempty(self.Type) && ~strcmp(self.Type.Shape, 'row') ...  % cannot remove single row for row vector
                    && size(self.Matrix,1) > 0 ...  % check if there is anything to remove
                    && ~isempty(self.GetSelectedRow())
                remrowstate = 'on';
            else
                remrowstate = 'off';
            end
            set(self.Menu.AddColumnBefore, 'Enable', addcolstate);
            set(self.Menu.AddColumnAfter, 'Enable', addcolstate);
            set(self.Menu.AddRowBefore, 'Enable', addrowstate);
            set(self.Menu.AddRowAfter, 'Enable', addrowstate);
            set(self.Menu.RemoveColumn, 'Enable', remcolstate);
            set(self.Menu.RemoveRow, 'Enable', remrowstate);
        end

        function self = OnShowTranspose(self, source, event) %#ok<INUSD>
            self.Transpose = ~self.Transpose;  % call setter method
        end

        function self = OnAddColumnBefore(self, source, event) %#ok<INUSD>
            index = self.GetSelectedColumn();
            if isempty(index)
                index = 0;
            end
            self.AddColumnAt(index);  % index ranges from 0 to column count-1
        end

        function self = OnAddColumnAfter(self, source, event) %#ok<INUSD>
            index = self.GetSelectedColumn();
            if ~isempty(index)
                self.AddColumnAt(index+1);  % index+1 ranges from 1 to column count
            else
                self.AddColumnAt(size(self.Matrix,2));
            end
        end

        function self = OnAddRowBefore(self, source, event) %#ok<INUSD>
            index = self.GetSelectedRow();
            if isempty(index)
                index = 0;
            end
            self.AddRowAt(index);  % index ranges from 0 to row count-1
        end

        function self = OnAddRowAfter(self, source, event) %#ok<INUSD>
            index = self.GetSelectedRow();
            if ~isempty(index)
                self.AddRowAt(index+1);  % index+1 ranges from 1 to row count
            else
                self.AddRowAt(size(self.Matrix,1));
            end
        end

        function self = OnRemoveColumn(self, source, event) %#ok<INUSD>
            index = self.GetSelectedColumn();  % index ranges from 0 to column count-1
            if isempty(index)
                return;
            end
            matrix = self.Matrix;
            matrix(:,index+1) = [];
            self.Matrix = matrix;
            data = self.JideTableModel.getDataVector();
            for k = 0 : size(matrix,1)-1
                row = data.elementAt(k);
                row.removeElementAt(index);
            end
            self.JideTableModel.setColumnCount(size(matrix,2));
            self.JideTable.repaint();
        end

        function self = OnRemoveRow(self, source, event) %#ok<INUSD>
            index = self.GetSelectedRow();  % index ranges from 0 to row count-1
            if isempty(index)
                return;
            end
            matrix = self.Matrix;
            matrix(index+1,:) = [];
            self.Matrix = matrix;
            self.JideTableModel.removeRow(index);
            self.JideTable.repaint();
        end

        function index = GetSelectedColumn(self)
        % Returns the currently selected column.
        %
        % Output arguments:
        % index:
        %    a zero-based index
            index = self.JideTable.getSelectedColumn();
            if index < 0  % no row is selected
                index = [];
            end
        end

        function index = GetSelectedRow(self)
        % Returns the currently selected row.
        %
        % Output arguments:
        % index:
        %    a zero-based index
            index = self.JideTable.getSelectedRow();
            if index < 0  % no row is selected
                index = [];
            end
        end

        function self = AddColumnAt(self, index)
        % Add new column at the specified index.
        % Items whose index is greater than the specified index are shifted
        % to the right.
        %
        % Input argument:
        % index:
        %    an insertion index between 0 (to insert as first column) and
        %    size(self.Matrix,2) (to insert as last column)
            matrix = self.Matrix;
            if isempty(matrix)
                ix1 = [];
                ix2 = [];
            else
                ix1 = 1 : index;
                ix2 = index+1 : size(matrix,2);
            end
            matrix = [ matrix(:,ix1) , zeros(size(matrix,1),1) , matrix(:,ix2) ];
            self.Matrix = matrix;
            data = self.JideTableModel.getDataVector();
            for k = 0 : size(matrix,1)-1
                row = data.elementAt(k);
                row.insertElementAt(self.Type.GetPrimitiveJavaValue(0), index);
            end
            self.JideTableModel.setColumnCount(size(matrix,2));
            self.JideTable.repaint();
        end

        function self = AddRowAt(self, index)
        % Add new row at the specified index.
        % Items whose index is greater than the specified index are shifted
        % downwards.
        %
        % Input argument:
        % index:
        %    an insertion index between 0 (to insert as first row) and
        %    size(self.Matrix,1) (to insert as last row)
            matrix = self.Matrix;
            if isempty(matrix)
                ix1 = [];
                ix2 = [];
            else
                ix1 = 1 : index;
                ix2 = index+1 : size(matrix,1);
            end
            matrix = [ matrix(ix1,:) ; zeros(1,size(matrix,2)) ; matrix(ix2,:) ];
            self.Matrix = matrix;
            data = self.Type.GetJavaVector(zeros(1,size(matrix,2)));
            self.JideTableModel.insertRow(index, data);
            self.JideTable.repaint();
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function Agal = mg_diff_setup(x,y)
%mg_diff_setup   GMG diffusion problem on square domain
%   Agal=mg_diff_setup(x,y)
%   input
%          x       x coordinate vector for coarse grid
%          y       y coordinate vector for coarse grid
%   output
%          Agal    discrete diffusion operator
%
%   IFISS function: AR; 19 November, 2001.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
n=length(x)-1; np=n/2; nq=n/4;
nvtx=(n+1)*(n+1);
[X,Y]=meshgrid(x,y);
xx=reshape(X',nvtx,1);
yy=reshape(Y',nvtx,1);
xy=[xx(:),yy(:)];
% assembly process
kx = 1;
ky = 1;
mel=0;
for j=1:np
   for i=1:np
      mref=(n+1)*(ky-1)+kx;
      pref=(np+1)*(j-1)+i;
      mel=mel+1;
      nvv(1) = mref;
      nvv(2) = mref+2;
      nvv(3) = mref+2*n+4;
      nvv(4) = mref+2*n+2;
      nvv(5) = mref+1;
      nvv(6) = mref+n+3;
      nvv(7) = mref+2*n+3;
      nvv(8)=  mref+n+1;
      nvv(9)=  mref+n+2;
      npp(1) = pref;
      npp(2) = pref+1;
      npp(3) = pref+np+2;
      npp(4) = pref+np+1;
      mv(mel,1:9)=nvv(1:9);
      mp(mel,1:4)=npp(1:4);
      kx = kx + 2;
   end
   ky = ky + 2;
   kx = 1;
end
%
% compute boundary vertices and edges
% four boundary edges
k1=find( xy(:,2)==-1 );
e1=[]; for k=1:mel, if any(mv(k,5)==k1), e1=[e1,k]; end, end
ef1=ones(size(e1));
%
k2=find( xy(:,1)==1  & xy(:,2)<1   & xy(:,2) >-1);
e2=[]; for k=1:mel, if any(mv(k,6)==k2), e2=[e2,k]; end, end
ef2=2*ones(size(e2));
%
k3=find( xy(:,2)==1 );
e3=[]; for k=1:mel, if any(mv(k,7)==k3), e3=[e3,k]; end, end
ef3=3*ones(size(e3));
%
k4=find( xy(:,1)==-1 & xy(:,2)<1   & xy(:,2) >-1 );
e4=[]; for k=1:mel, if any(mv(k,8)==k4), e4=[e4,k]; end, end
ef4=4*ones(size(e4));
%
bound=sort([k1;k2;k3;k4]);
mbound=[e1',ef1';e2',ef2';e3',ef3';e4',ef4'];
%
% set up matrices for Q1 approximation
[ev,ebound]=mg_q1grid(x,y,xy,mv,bound,mbound);
[A,M,fdummy] = mg_q1diff(xy,ev);
%
% impose zero boundary conditions
Agal = mg_zerobc(A,xy,bound);


text:
function net = load_cnn(fparams, im_size)

net = load(['networks/' fparams.nn_name]);
net = vl_simplenn_tidy(net);
net.layers = net.layers(1:max(fparams.output_layer));

if strcmpi(fparams.input_size_mode, 'cnn_default')
    base_input_sz = net.meta.normalization.imageSize(1:2);
elseif strcmpi(fparams.input_size_mode, 'adaptive')
    base_input_sz = im_size(1:2);
else
    error('Unknown input_size_mode');
end

net.meta.normalization.imageSize(1:2) = round(base_input_sz .* fparams.input_size_scale);
net.meta.normalization.averageImageOrig = net.meta.normalization.averageImage;

if isfield(net.meta,'inputSize')
    net.meta.inputSize = base_input_sz;
end

if size(net.meta.normalization.averageImage,1) > 1 || size(net.meta.normalization.averageImage,2) > 1
    net.meta.normalization.averageImage = imresize(single(net.meta.normalization.averageImage), net.meta.normalization.imageSize(1:2));
end

net.info = vl_simplenn_display(net);
end

author: martin-danelljan

repo: ECO


text:
% Copyright (C) 1995-2009, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function b = ilabeltest
    b = [
    0 0 0 1 1 0 0 0 0 0
    0 0 0 1 1 0 1 1 1 0
    0 1 0 0 0 0 0 1 0 0
    0 1 0 1 0 1 0 1 0 0
    0 1 1 1 1 1 1 1 1 0
    0 0 0 0 0 1 0 0 1 0
    0 1 1 0 1 1 1 1 1 0
    0 0 0 0 0 0 0 0 0 0];


text:
function value = iif(condition, value1, value2)
% SYNTAX:
%   value = iif(condition, value1, value2)
% DESCRIPTION:
%   return value1 if condition is true otherwise value2

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:       Andrea Gatti
%  Contributors:     Andrea Gatti ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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.
%
%   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, see <http://www.gnu.org/licenses/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

    if condition
        value = value1;
    else
        value = value2;
    end
end

author: goGPS-Project


text:
function double_c_plot ( x, y, c )

%*****************************************************************************80
%
%% DOUBLE_C_PLOT displays the "double C" data.
%
%  Discussion:
%
%    The data comprises a "C" shape and its reverse, which do not intersect,
%    but which nestle together fairly closely.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%   Input, real X(N), Y(N), the coordinates of the data points.  The points
%   have been shuffled so that data belonging to the two components is
%   interleaved.
%
%   Input, integer C(N), is 1 or 2, depending on which set the corresponding
%   data point belongs to.
%
  i1 = ( c == 1 );
  i2 = ( c == 2 );

  n1 = sum ( i1 );
  n2 = sum ( i2 );

  plot ( x(i1), y(i1), 'r.', ...
         x(i2), y(i2), 'b.' )

  xlabel ( '<-- X -->' );
  ylabel ( '<-- Y -->' );
  title ( sprintf ( 'Double C Data, N1 (red) = %d, N2 (blue) = %d', n1, n2 ) );
  grid on
  axis equal

  return
end


text:
function res = le(a,b)
%LE           Implements  a <= b  elementwise for intervals a and b
%
%  if true,  a  is definitely less than or equal to  b
%

% written  10/16/98     S.M. Rump
% modified 11/30/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if ~isa(a,'intval')
    a = intval(a);
  end
  if ~isa(b,'intval')
    b = intval(b);
  end

  if a.complex | b.complex
    res = real(sup(a)) <= real(inf(b)) & imag(sup(a)) <= imag(inf(b)) ;
  else
    res = sup(a) <= inf(b) ;
  end

repo: OpenMAS


text:
function [ t, wts ] = h_quadrature_rule ( nt )

%*****************************************************************************80
%
%% H_QUADRATURE_RULE: quadrature for H(i,x).
%
%  Discussion:
%
%    H(i,x) is the physicist's Hermite polynomial of degree I.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NT, the order of the rule.
%
%    Output, real T(NT,1), WTS(NT,1), the points and weights of the rule.
%
  aj = zeros ( nt, 1 );
  bj = sqrt ( ( 1 : nt )' / 2.0 );
  wts = zeros ( nt, 1 );
  wts(1,1) = sqrt ( sqrt ( pi ) );

  [ t, wts ] = imtqlx ( nt, aj, bj, wts );

  wts(1:nt,1) = wts(1:nt,1).^2;

  return
end


text:
function [hrf, fit, e, param] = Fit_sFIR(tc, TR, Run, T, mode)
% Fits FIR and smooth FIR model
%
% :Usage:
% ::
%
%     function [hrf, fit, e, param] = Fit_sFIR(tc,TR,Runs,T,mode)
%
% :Inputs:
%
%   **tc:**
%        time course
%
%   **TR:**
%        time resolution
%
%   **Runs:**
%        expermental design
%
%   **T:**
%        length of estimated HRF
%
%   **mode:**
%        FIR or smooth FIR
%
%        Options:
%           0 - standard FIR
%
%           1 - smooth FIR
%
% :Outputs:
%
%   **hrf:**
%        estimated hemodynamic response function
%
%   **fit:**
%        estimated time course
%
%   **e:**
%        residual time course
%
%   **param:**
%        estimated amplitude, height and width
%
% ..
%    Created by Martin Lindquist on 10/02/09
%    Last edited: 05/17/13 (ML)
% ..

numstim = length(Run);
len = length(Run{1});
t=1:TR:T;
tlen = length(t);

Runs = zeros(len,numstim);
for i=1:numstim,
    Runs(:,i) = Run{i};
end;

[DX] = tor_make_deconv_mtx3(Runs,tlen,1);

if mode == 1

    C=(1:tlen)'*(ones(1,tlen));
    h = sqrt(1/(7/TR));                       % 7 seconds smoothing - ref. Goutte

    v = 0.1;
    sig = 1;

    R = v*exp(-h/2*(C-C').^2);
    RI = inv(R);
    MRI = zeros(numstim*tlen+1);
    for i=1:numstim,
        MRI(((i-1)*tlen+1):(i*tlen),((i-1)*tlen+1):(i*tlen)) = RI;
    end;

    b = inv(DX'*DX+sig^2*MRI)*DX'*tc;
    fit = DX*b;
    e = tc - DX*b;

elseif mode == 0

    b = pinv(DX)*tc;
    fit = DX*b;
    e = tc - DX*b;

end

hrf =zeros(tlen,numstim);
param = zeros(3,numstim);

for i=1:numstim,
    hrf(:,i) = b(((i-1)*tlen+1):(i*tlen))';
    param(:,i) = get_parameters2(hrf(:,i),(1:tlen));
end;

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

function MU_funcWaterfall(Temp,Event,handles)
handles = guidata(handles.MU_matrix_display);

figure;
waterfall(double(handles.BMatrix));
colormap(handles.V.Color_map);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
clearvars;
close all;
% clc;

rng default;

M = 300;
N = 1000;
K = 60;

A = randn(M, N);
x0 = zeros(N,1);

permutation = randperm(N);
indices = permutation(1:K);
% non-zero entries
x0(indices) = randn(K,1);

b0 = A * x0;
% noise variance
sigma = 0.1;
e = sigma * randn(M, 1);
% add noise
b = b0 + e;

% orthonormalize the rows of A
% perform QR decomposition of rows of A
[Q, R] = qr(A', 0);
% keep the orthogonal rows
A = Q';
% change b accordingly to the new coordinates
b = R'\b;

options.verbose = 0;
options.max_iterations = 200;
options.tolerance = 5e-3;
solver = spx.pursuit.single.L1_ADMM_YZ(A, options);
mu = sigma;
x = solver.solve_bpdn_l2(b, mu);

r = x - x0;
max_diff = max(abs(r));
rel_error = norm(x-x0) /norm(x0);
iterations = solver.details.iterations(1);
elapsed_time = solver.details.elapsed_times(1);
fprintf('Iterations: %d, Relative error: %e, max diff: %.4f, time: %.4f seconds\n', ...
    iterations, rel_error, max_diff, elapsed_time);

print = 1;
if print
    subplot(211);
    stem(x0, '.');
    subplot(212);
    stem(x, '.');
    figure;
    iterations = solver.details.iterations(1);
    iters = 1:iterations;
    primal_objectives = solver.details.primal_objectives(iters, 1);
    dual_objectives = solver.details.dual_objectives(iters, 1);
    plot(iters, primal_objectives);
    hold on;
    plot(iters, dual_objectives);
    xlabel('Iterations');
    ylabel('Objective Value');
    legend({'Primal Objective', 'Dual Objective'});
    grid on;
end

cvx = 0;
if cvx
% time to compare with a CVX implementation
cvx_solver = spx.pursuit.single.BasisPursuit(A, b);
tstart = tic;
x = cvx_solver.solve_lasso(1/(2*mu));
elapsed_time = toc(tstart);
r = x - x0;
max_diff = max(abs(r));
rel_error = norm(x-x0) /norm(x0);
fprintf('Relative error: %e, max diff: %.4f, time: %.4f seconds\n', ...
    rel_error, max_diff, elapsed_time);
end

author: indigits

repo: sparse-plex


text:
function plot_q4_mesh ( node_num, element_num, node_xy, element_node, ...
  node_show, element_show, output_filename )

%*****************************************************************************80
%
%% TPLOT_Q4_MESH plots a Q4 mesh.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the nodes.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(4,ELEMENT_NUM), the nodes
%    that form the elements.
%
%    Input, logical NODE_SHOW:
%    0, do not show the nodes.
%    1, show the nodes.
%    2, show the nodes, and label them.
%
%    Input, logical ELEMENT_SHOW,
%    0, do not show the elements.
%    1, show the elements.
%    2, show the elements, and label them.
%
%    Input, string OUTPUT_FILENAME, the name of the output file.
%
  element_order = 4;

  x_ps_max = 576;
  x_ps_max_clip = 594;
  x_ps_min = 36;
  x_ps_min_clip = 18;
  y_ps_max = 666;
  y_ps_max_clip = 684;
  y_ps_min = 126;
  y_ps_min_clip = 108;

  date_time = timestring ( );
%
%  We need to do some figuring here, so that we can determine
%  the range of the data, and hence the height and width
%  of the piece of paper.
%
  x_max = max ( node_xy(1,1:node_num) );
  x_min = min ( node_xy(1,1:node_num) );
  x_scale = x_max - x_min;

  x_max = x_max + 0.05 * x_scale;
  x_min = x_min - 0.05 * x_scale;
  x_scale = x_max - x_min;

  y_max = max ( node_xy(2,1:node_num) );
  y_min = min ( node_xy(2,1:node_num) );
  y_scale = y_max - y_min;

  y_max = y_max + 0.05 * y_scale;
  y_min = y_min - 0.05 * y_scale;
  y_scale = y_max - y_min;

  if ( x_scale < y_scale )

    delta = round ( ( x_ps_max - x_ps_min ) ...
      * ( y_scale - x_scale ) / ( 2.0 * y_scale ) );

    x_ps_max = x_ps_max - delta;
    x_ps_min = x_ps_min + delta;

    x_ps_max_clip = x_ps_max_clip - delta;
    x_ps_min_clip = x_ps_min_clip + delta;

    x_scale = y_scale;

  elseif ( y_scale < x_scale )

    delta = round ( ( y_ps_max - y_ps_min ) ...
      * ( x_scale - y_scale ) / ( 2.0 * x_scale ) );

    y_ps_max      = y_ps_max - delta;
    y_ps_min      = y_ps_min + delta;

    y_ps_max_clip = y_ps_max_clip - delta;
    y_ps_min_clip = y_ps_min_clip + delta;

    y_scale = x_scale;

  end
%
%  Plot the mesh.
%
  file_unit = fopen ( output_filename, 'wt' );

  if ( file_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PLOT_Q4_MESH - Fatal error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'PLOT_Q4_MESH - Fatal error!' );
  end

  fprintf ( file_unit, '%%!PS-Adobe-3.0 EPSF-3.0\n' );
  fprintf ( file_unit, '%%%%Creator: plot_q4_mesh.m\n' );
  fprintf ( file_unit, '%%%%Title: %s\n', output_filename );
  fprintf ( file_unit, '%%%%CreationDate: %s\n', date_time );
  fprintf ( file_unit, '%%%%Pages: 1\n' );
  fprintf ( file_unit, '%%%%BoundingBox:  %d  %d  %d  %d\n', ...
    x_ps_min, y_ps_min, x_ps_max, y_ps_max );
  fprintf ( file_unit, '%%%%Document-Fonts: Times-Roman\n' );
  fprintf ( file_unit, '%%%%LanguageLevel: 1\n' );
  fprintf ( file_unit, '%%%%EndComments\n' );
  fprintf ( file_unit, '%%%%BeginProlog\n' );
  fprintf ( file_unit, '/inch {72 mul} def\n' );
  fprintf ( file_unit, '%%%%EndProlog\n' );
  fprintf ( file_unit, '%%%%Page: 1 1\n' );
  fprintf ( file_unit, 'save\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%  Set the RGB color to very light gray.\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '0.900  0.900  0.900 setrgbcolor\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%  Draw a gray border around the page.\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, 'newpath\n' );
  fprintf ( file_unit, '  %d  %d  moveto\n', x_ps_min, y_ps_min );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_max, y_ps_min );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_max, y_ps_max );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_min, y_ps_max );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_min, y_ps_min );
  fprintf ( file_unit, 'stroke\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%  Set the RGB color to black.\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '0.000  0.000  0.000 setrgbcolor\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%  Set the font and its size.\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '/Times-Roman findfont\n' );
  fprintf ( file_unit, '0.50 inch scalefont\n' );
  fprintf ( file_unit, 'setfont\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%  Print a title.\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%210  702  moveto\n' );
  fprintf ( file_unit, '%%(Triangulation)  show\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%  Define a clipping polygon.\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, 'newpath\n' );
  fprintf ( file_unit, '  %d  %d  moveto\n', x_ps_min_clip, y_ps_min_clip );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_max_clip, y_ps_min_clip );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_max_clip, y_ps_max_clip );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_min_clip, y_ps_max_clip );
  fprintf ( file_unit, '  %d  %d  lineto\n', x_ps_min_clip, y_ps_min_clip );
  fprintf ( file_unit, 'clip newpath\n' );
%
%  Draw the nodes.
%
  if ( node_num <= 200 )
    circle_size = 5;
  elseif ( node_num <= 500 )
    circle_size = 4;
  elseif ( node_num <= 1000 )
    circle_size = 3;
  elseif ( node_num <= 5000 )
    circle_size = 2;
  else
    circle_size = 1;
  end

  if ( 1 <= node_show )

    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Draw filled dots at the nodes.\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Set the RGB color to blue.\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '0.000  0.150  0.750 setrgbcolor\n' );
    fprintf ( file_unit, '%%\n' );

    for node = 1 : node_num

      x_ps = floor ( ...
        ( ( x_max - node_xy(1,node)         ) * x_ps_min ...
        + (         node_xy(1,node) - x_min ) * x_ps_max ) ...
        / ( x_max                   - x_min ) );

      y_ps = floor ( ...
        ( ( y_max - node_xy(2,node)         ) * y_ps_min   ...
        + (         node_xy(2,node) - y_min ) * y_ps_max ) ...
        / ( y_max                   - y_min ) );

      fprintf ( file_unit, ...
        '  newpath  %d  %d  %d 0 360 arc closepath fill\n', ...
        x_ps, y_ps, circle_size );

    end

  end
%
%  Label the nodes.
%
  if ( 2 <= node_show )

    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Label the nodes.\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Set the RGB color to darker blue.\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '0.000  0.250  0.850 setrgbcolor\n' );
    fprintf ( file_unit, '/Times-Roman findfont\n' );
    fprintf ( file_unit, '0.20 inch scalefont\n' );
    fprintf ( file_unit, 'setfont\n' );
    fprintf ( file_unit, '%%\n' );

    for node = 1 : node_num

      x_ps = floor ( ...
        ( ( x_max - node_xy(1,node)         ) * x_ps_min ...
        + (         node_xy(1,node) - x_min ) * x_ps_max ) ...
        / ( x_max                   - x_min ) );

      y_ps = floor ( ...
        ( ( y_max - node_xy(2,node)         ) * y_ps_min   ...
        + (         node_xy(2,node) - y_min ) * y_ps_max ) ...
        / ( y_max                   - y_min ) );

      fprintf ( file_unit, '  %d  %d  moveto (%d) show\n', ...
        x_ps, y_ps+5, node );

    end

  end
%
%  Draw the elements.
%
  if ( 1 <= element_show )

    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Set the RGB color to red.\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '0.900  0.200  0.100 setrgbcolor\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Draw the elements.\n' );
    fprintf ( file_unit, '%%\n' );

    for element = 1 : element_num

      fprintf ( file_unit, 'newpath\n' );

      for i = 1 : element_order + 1

        e = i4_wrap ( i, 1, element_order )

        node = element_node(e,element);

        x_ps = floor ( ...
          ( ( x_max - node_xy(1,node)         ) * x_ps_min ...
          + (         node_xy(1,node) - x_min ) * x_ps_max ) ...
          / ( x_max                   - x_min ) );

        y_ps = floor ( ...
          ( ( y_max - node_xy(2,node)         ) * y_ps_min ...
          + (         node_xy(2,node) - y_min ) * y_ps_max ) ...
          / ( y_max                   - y_min ) );

        if ( i == 1 )
          fprintf ( file_unit, '  %d  %d  moveto\n', x_ps, y_ps );
        else
          fprintf ( file_unit, '  %d  %d  lineto\n', x_ps, y_ps );
        end

      end

      fprintf ( file_unit, 'stroke\n' );

    end

  end
%
%  Label the elements.
%
  if ( 2 <= element_show )

    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Label the elements.\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '%%  Set the RGB color to darker red.\n' );
    fprintf ( file_unit, '%%\n' );
    fprintf ( file_unit, '0.950  0.250  0.150 setrgbcolor\n' );
    fprintf ( file_unit, '/Times-Roman findfont\n' );
    fprintf ( file_unit, '0.20 inch scalefont\n' );
    fprintf ( file_unit, 'setfont\n' );
    fprintf ( file_unit, '%%\n' );

    for element = 1 : element_num

      ave_x = 0.0;
      ave_y = 0.0;

      for i = 1 : element_order

        node = element_node(i,element);
        ave_x = ave_x + node_xy(1,node);
        ave_y = ave_y + node_xy(2,node);
      end

      ave_x = ave_x / element_order;
      ave_y = ave_y / element_order;

      x_ps = floor ( ...
        ( ( x_max - ave_x         ) * x_ps_min ...
        + (         ave_x - x_min ) * x_ps_max ) ...
        / ( x_max         - x_min ) );

      y_ps = floor ( ...
        ( ( y_max - ave_y         ) * y_ps_min ...
        + (         ave_y - y_min ) * y_ps_max ) ...
        / ( y_max         - y_min ) );

      fprintf ( file_unit, '  %d  %d  moveto (%d) show\n', ...
        x_ps, y_ps, element );

    end
  end

  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, 'restore  showpage\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%  End of page.\n' );
  fprintf ( file_unit, '%%\n' );
  fprintf ( file_unit, '%%%%Trailer\n' );
  fprintf ( file_unit, '%%%%EOF\n' );

  fclose ( file_unit );

  return
end


text:
function coeffs=string2MultiDimPoly(theString)
%%STRING2MULTIDIMPOLY Given a string representing a multidimensional
%           polynomial equation, convert it into a hyoermatrix
%           representation of the polynomial, as can be used in function
%           such as polyValMultiDim and polyDerMultiDim.
%
%INPUTS: theString The string representing the multivariate polynomial.
%           All coefficients come before variables. Variables are all x
%           followed by a number, such as x3. All numbers for variables
%           must be >0. The format of this string is the same as that
%           returned by multiDimPolyMat2String. For example,
%           '60*x2^2-4*x2^3+72*x1*x2+12' is a valid string.
%
%OUTPUTS: coeffs A hypermatrix of the coefficients for the multivariate
%                polynomial. These are arranged such that
%                coeffs(a1,a2,a3...an) corresponds to the coefficient of an
%                x1^(a1-1)*x2^(a2-1)*x3^(a3-1)...xn^(an-1) term.  Thus, the
%                number of indices coeffs takes is equal to the
%                dimensionality of x (not counting singleton dimensions at
%                the end of coeffs). Note that this ordering is the reverse
%                that used in the 1D polyval function that is built into
%                Matlab. The number of elements for each index in coeffs is
%                the maximum order of that dimension +1.
%
%This function just calls string2Terms and terms2MultiDimPolyMat.
%
%January 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

coeffs=terms2MultiDimPolyMat(string2Terms(theString));

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function frameloc = ufmf_write_frame(fp,stamp,x0,y0,w,h,val)

frameloc = ftell(fp);
% write the chunk id (points=1)
fwrite(fp,1,'uchar');
% write timestamp
fwrite(fp,stamp,'double');
% write number of points
npts = numel(x0);
fwrite(fp,npts,'uint32');

for i = 1:npts,
  % write x, y, width, height
  fwrite(fp,[x0(i),y0(i),w(i),h(i)],'ushort');
  % write region intensities
  fwrite(fp,val{i},'uint8');
end

author: kristinbranson

repo: JAABA


text:
function gpplot(x, y, e, varargin)
% function gpplot(x, y, e, flipaxis)
%
% x = inputs (e.g., time)
% y = expectated values of the outputs
% e = errors of the outputs

warning('Deprecated, use errorplot')

opts = struct('pseudoinputs', []);
[ opts, errmsg, wrnmsg ] = argschk( opts, varargin{:} );
if ~isempty(errmsg), error( errmsg ), end
if ~isempty(wrnmsg), warning( wrnmsg ), end

% $$$ if nargin < 4
% $$$   flipaxis = false;
% $$$ end
x = x(:);
y = y(:);
%e = sqrt(diag(Cov));
X = [x(1:end); x(end:-1:1)];
Y = [y(1:end)+e; y(end:-1:1)-e(end:-1:1)];
C = [.75,.75,.75];
% $$$ if flipaxis
% $$$   fill(Y,X,C,'EdgeColor',C);
% $$$   hold on
% $$$   plot(y,x,'k');
% $$$ else
fill(X,Y,C,'EdgeColor',C);
hold on
plot(x,y,'k');
if ~isempty(opts.pseudoinputs)
  yl = ylim;
  py = yl(1) + 0.8*(yl(2)-yl(1));
  plot(opts.pseudoinputs, py, 'k+')
end
% $$$ end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function unzipAll(varargin)

% function unzipAll(pathName,pathNameOutput)
% ------------------------------------------------------------------------
%
% This function unzips all zip files in the folder pathName. The first
% input is the folder path name, the following inputs are possible inputs
% of the unzip function.
%
% See also: dicomanon
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2018/02/09: Created
%------------------------------------------------------------------------

%%
switch nargin
    case 1
        pathName=varargin{1};
        varOutput=pathName;
    case 2
        pathName=varargin{1};
        varOutput=varargin{2};
end
%%

if ischar(varOutput) %Process only input folder
    files = dir(fullfile(pathName,'*.zip'));
    files={files(1:end).name};
    files=sort(files(:));
    if ~isempty(files)
        if numel(files)>1
            hw = waitbar(0,'Unzipping files...');
            for q=1:1:numel(files)
                fileName=fullfile(pathName,files{q});
                unzip(fileName,varOutput);
                waitbar(q/numel(files),hw,['Unzipping files...',num2str(round(100.*q/numel(files))),'%']);
            end
            close(hw)
        else
            fileName=fullfile(pathName,files{1});
            unzip(fileName,varOutput);
        end
    end
else %Process input folder and all subfolders
    %Process current
    unzipAll(pathName,pathName);

    %Process sub-folders
    [pathNames]=getSubPaths(pathName);
    if ~isempty(pathNames)
        for q=1:1:numel(pathNames)
            pathNameNow=pathNames{q};
            unzipAll(pathNameNow,pathNameNow);
        end
    end
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% 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.
%
% 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, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [idx,fstr] = PVcipmatch(clist,str)
% function [Idx,Fstr] = PVcipmatch(Clist,Str)
%
% DESCRIPTION
%   Case insensitive partial match.
%
% INPUTS
%   Clist:
%   Str:
%
% OUTPUTS
%   Idx:
%   Fstr:
%
% SYNTAX
%   [Idx,Fstr] = PVcipmatch(Clist,Str)

% 6/7/2002: Written by Packard

lstr = length(str);
tf = find(strncmpi(clist,str,lstr));
if length(tf)==1
    idx = tf;
    fstr = clist{idx};
elseif isempty(tf)
    idx = [];
    fstr = [];
else
    cnt = 0;
    for i=1:length(tf)
        if length(clist{tf(i)})==lstr
            tfuse = tf(i);
            cnt = cnt + 1;
        end
    end
    if cnt==1
        idx = tfuse;
        fstr = clist{idx};
    else
        tf = find(strncmp(clist,str,lstr));
        if length(tf)==1
            idx = tf;
            fstr = clist{idx};
        elseif isempty(tf)
            idx = [];
            fstr = [];
        else
            cnt = 0;
            for i=1:length(tf)
                if length(clist{tf(i)})==lstr
                    tfuse = tf(i);
                    cnt = cnt + 1;
                end
            end
            if cnt==1
                idx = tfuse;
                fstr = clist{idx};
            else
                error(['Ambiguous property: ''' str '''.']);
            end
        end
    end
end

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function pts = critical_points(V,F,u,tol)
%REVERSE_SUBDIVISION Find out which face contains the critical points of a
%function
%
% pts = critical_points(V,F,u,tol);
%
% Inputs:
%  V,F  the coarse input mesh
%  u  the function for which to find critical points
%  tol  the tolerance for writical points (if the gradient norm is below this,
%                                          we are at a critical point)
% Outputs:
%  pts  a list of face indices into the rows of F that tell us which faces
%       contain the critical points.
%

G = grad(V,F);
pts = ...

end

author: odedstein

repo: sgi-introduction-course


text:
%   TTeMPS Toolbox.
%   Michael Steinlechner, 2013-2016
%   Questions and contact: michael.steinlechner@epfl.ch
%   BSD 2-clause license, see LICENSE.txt
function [X, residuum, cost, times] = alsLinsolve_rankOne( L, F, X, opts )

t_start = tic();
% set default opts
if ~exist( 'opts', 'var');       opts = struct();       end
if ~isfield( opts, 'nSweeps');   opts.nSweeps = 4;      end
if ~isfield( opts, 'solver');    opts.solver = 'pcg';   end

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

normF = norm(F);
g = apply(L, X) - F;
cost = cost_function_res( X, g );
residuum = norm( g ) / normF;
times = toc(t_start);

X = orthogonalize(X, 1);
for sweep = 1:opts.nSweeps
    % ====================================================================
    % LEFT-TO-RIGHT SWEEP
    % ====================================================================
    disp( ['STARTING SWEEP ', num2str(sweep), ' from left to right'] )
    disp( '===========================================================')
    for idx = 1:d-1
        disp( ['Current core: ', num2str(idx)] )

        Fi = contract( X, F, idx );
        sz = [X.rank(idx), X.size(idx), X.rank(idx+1)];

        if strcmpi( opts.solver, 'direct' )
            % if system very small
            Li = contract( X, apply(L,X), idx );
            Ui = Li \ Fi(:);
            X.U{idx} = reshape( Ui, sz );

        elseif strcmpi( opts.solver, 'pcg' )

            [left, right] = Afun_prepare( L, X, idx );
            B1 =  prepare_precond( L.A{1}, X, idx );

            Ui = pcg( @(y) Afun( L, y, idx, sz, left, right), ...
                     Fi(:), ...
                     1e-10, 1000, ...
                     @(y) apply_precond( L.A{1}, B1, y, sz ), [],...
                     X.U{idx}(:) );

            X.U{idx} = reshape( Ui, sz );

        elseif strcmpi( opts.solver, 'diag' )
            X.U{idx} = solve_inner( L.L0, X, Fi, idx );

        else
            error( 'Unknown opts.solver type. Use either ''direct'', ''pcg'' (default) or ''diag''.' )
        end

        X = orth_at( X, idx, 'left', true );

        g = apply(L, X) - F;
        residuum = [residuum; norm( g ) / normF];
        cost = [cost; cost_function_res( X, g )];
        times = [times; toc(t_start)];
    end

    % ====================================================================
    % RIGHT-TO-LEFT
    % ====================================================================
    disp( 'Starting right-to-left half-sweep:')
    for idx = d:-1:2
        disp( ['Current core: ', num2str(idx)] )

        Fi = contract( X, F, idx );
        sz = [X.rank(idx), X.size(idx), X.rank(idx+1)];

        if strcmpi( opts.solver, 'direct' )
            % if system very small
            Li = contract( X, apply(L, X), idx );
            Ui = Li \ Fi(:);
            X.U{idx} = reshape( Ui, sz );

        elseif strcmpi( opts.solver, 'pcg' )

            [left, right] = Afun_prepare( L, X, idx );
            B1 =  prepare_precond( L.A{1}, X, idx );

            Ui = pcg( @(y) Afun( L, y, idx, sz, left, right), ...
                     Fi(:), ...
                     1e-10, 1000, ...
                     @(y) apply_precond( L.A{1}, B1, y, sz ), [],...
                     X.U{idx}(:) );

            X.U{idx} = reshape( Ui, sz );

        elseif strcmpi( opts.solver, 'diag' )
            X.U{idx} = solve_inner( L.L0, X, Fi, idx );

        else
            error( 'Unknown opts.solver type. Use either ''direct'', ''pcg'' (default) or ''diag''.' )
        end

        X = orth_at( X, idx, 'right', true );

        g = apply(L, X) - F;
        residuum = [residuum; norm( g ) / normF];
        cost = [cost; cost_function_res( X, g )];
        times = [times; toc(t_start)];
    end

end

end

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

function res = cost_function_res( X, res )
res = 0.5*innerprod( X, res );
end

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

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

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

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

    res = tmp(:);
end

function B1 = prepare_precond( L0, X, idx )

    if idx == 1
        B1 = [];
        return
    end

    n = size(L0, 1);
    r = X.rank;

    X1 = matricize( X.U{1}, 2);
    Y = X;
    Y.U{1} = tensorize( L0*X1, 2, [r(1), n(1), r(2)] );
    B1 = innerprod( X, Y, 'LR', idx-1);
end

function res = apply_precond( L0, B1, rhs, sz )

    n = size(L0, 1);
    rhs = reshape( rhs, sz );
    if isempty(B1) %idx == 1
        res = L0 \ unfold( rhs, 'left' );
        res = reshape( res, sz );
    else
        res = B1 \ unfold(rhs, 'right');
        res = reshape( res, sz );
    end
    res = res(:);
end

author: NicolasBoumal

repo: manopt


text:
function kern = componentKernReadParamsFromFID(kern, FID, version)

% COMPONENTKERNREADPARAMSFROMFID Read a component based kernel from a C++ file.
% FORMAT
% DESC reads the components fo a kernel from a file written by C++ code.
% ARG kern : the base kernel to add components to.
% ARG FID : the input file stream to add from.
% RETURN kern : the kernel with the components added in.
%
% SEEALSO : modelReadFromFID, kernReadFromFID
%
% COPYRIGHT : Neil D. Lawrence, 2005, 2006, 2008

% KERN

kern.inputDimension = readIntFromFID(FID, 'inputDim');
numParams = readIntFromFID(FID, 'numParams');
numKerns = readIntFromFID(FID, 'numKerns');

for i=1:numKerns
  if version > 0.11
    kern.comp{i} = modelReadFromFID(FID);
  else
    kern.comp{i} = kernReadFromFID(FID, version);
  end
end

for i = 1:length(kern.comp)
  kern.nParams = kern.nParams + kern.comp{i}.nParams;
  kern.comp{i}.index = [];
end
kern.paramGroups = speye(kern.nParams);

if strcmp(kern.type, 'cmpnd')
  % Summarise the total white variance in the field whiteVariance.
  kern.whiteVariance = 0;
  for i = 1:length(kern.comp)
    if strcmp(kern.comp{i}.type, 'white')
      kern.whiteVariance = kern.whiteVariance + kern.comp{i}.variance;
    else
      if(isfield(kern.comp{i}, 'whiteVariance'))
        kern.whiteVariance = kern.whiteVariance + ...
            kern.comp{i}.whiteVariance;
      end
    end
  end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ torque ] = getExternalTorqueAtJoint( tKuka,k )
%% This function is used to get the external torque at some joint, for the KUKA iiwa robot.

%% Syntax:
% [ torque ] = getExternalTorqueAtJoint( tKuka,k )

%% About:
% This function is used to get external torque (due to external forces) at
% some joint of the robot.

%% Arreguments:
% tKuka: is the TCP/IP connection
% k: is joint index. A number from 1 to 7.
% torque: is the value of the torque measured at joint (k) due to external
% force/moment acting on the robot. Unit (Newton.meter)

% Copyright, Mohammad SAFEEA, 11th of May 2017

	if(k>7)
	    disp('Error KUKA has only 7 joints');
	    return;
	end

	if(k<1)
	disp('Error, joint number shall be more than zero and less than eight');
	    return;
	end

	[ torques ] = getJointsExternalTorques( tKuka );
	torque=torques{k};
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function R = unidrndKPM(min, max, nr, nc)

if nargin < 3
  nr = 1; nc = 1;
end

R = unidrnd(max-min+1, nr, nc) + (min-1);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef (Sealed) PTKUtils < handle
    % PTKUtils. A script to update the PTK codebase via git
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2015.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    methods (Static)
        function RunScript(script_name, varargin)
            % Runs a PTKScript from the Scripts folder

            m = PTKMain();
            m.RunScript(script_name, varargin{:});
        end

        function Recompile()
            % Forces recompilation of mex files

            m = PTKMain();
            m.Recompile();
        end

        function RebuildDatabase()
            % Forces a rebuild of the image database

            m = PTKMain();
            m.RebuildDatabase();
        end

        function ClearMemory()
            % Clears out memory caches

            close all;
            clear all classes;
            path(pathdef);
        end

        function ClearDisk()
            % Clears out disk and memory caches for temporary files

            m = PTKMain();
            m.DeleteCacheForAllDatasets();
        end

        function ReloadPlugins()
            % Forces reload of plugin information. Useful if you have modified properties of the plugin

            m = PTKMain();
            m.ReloadPlugins();
        end

        function Clear()
            PTKUtils.ClearMemory();
            PTKAddPaths;
        end

    end

    methods (Access = private)
        function obj = PTKUtils()
            disp('PTKUtils contains method to help maintain your PTK');
            disp('PTKUtils.RunScript(script_name, parameters) will run a PTKScript with the specified parameters');
            disp('PTKUtils.RebuildDatabase() will reimport your imaging files and correct problems with the database');
            disp('PTKUtils.ReloadPlugins() will reload plugins to fetch property changes');
            disp('PTKUtils.Recompile() will recompile your mex files');
            disp('PTKUtils.ClearMemory() will clear all the memory caches');
            disp('PTKUtils.ClearDisk() will clear the disk caches to retrieve disk space');
            ClearDisk
        end
    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
function permutations = dtiBootGetPermutations(permuteMatrix, numPermutations, bRepetition)
%
% permutations = dtiBootGetPermutations(permuteMatrix, numPermutations, [bRepetition=0])
%
%
% HISTORY:
% 2007.06.22 RFD wrote it

if ieNotDefined('bRepetition')
    bRepetition = 0;
end

if(bRepetition==0)
    numVols = permuteMatrix;
    permutations = zeros(numVols,numPermutations);
    for(ii=1:numVols)
        %permutations(ii,:) = randsample(1:numVols,numPermutations,true);
        permutations(ii,:) = ceil(numVols*rand(1,numPermutations));
    end
else
    numVols = size(permuteMatrix,1);
    permutations = zeros(numVols,numPermutations);
    for(ii=1:numVols)
        permutations(ii,:) = permuteMatrix{ii}(ceil(length(permuteMatrix{ii}).*rand(1,numPermutations)));
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function y = im2jpeg2k(x,n,q)
%IM2JPEG2K Compresses an image using a JPEG 2000 approximation.
%   Y = IM2JPEG2K(X,N,Q) compresses image X using an N-scale JPEG 2K
%   wavelet transform, implicit or explicit coefficient quantization,
%   and Huffman symbol coding augmented by zero run-length coding. If
%   quantization vector Q contains two elements, they are assumed to be
%   implicit quantization parameters; else, it is assumed to contain
%   explicit subband step sizes.  Y is an encoding structure containing
%   Huffman-encoded data and additional parameters needed by JPEG2K2IM
%   for decoding.
%
%   See also JPEG2K2IM.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

narginchk(3,3);         % Check input arguments

if ~ismatrix(x) || ~isreal(x) || ~isnumeric(x) || ~isa(x,'uint8')
   error('The input must be a UINT8 image.');
end

if length(q) ~= 2 && length(q) ~= 3 * n + 1
   error('The quantization step size vector is bad.');
end

% Level shift the input and compute its wavelet transform.
x = double(x) - 128;
[c, s] = wavefast(x,n,'jpeg9.7');

% Quantize the wavelet coefficients.
q = stepsize(n,q);
sgn = sign(c);     sgn(sgn == 0) = 1;     c = abs(c);
for k = 1:n
   qi = 3 * k - 2;
   c = wavepaste('h',c,s,k,wavecopy('h',c,s,k) / q(qi));
   c = wavepaste('v',c,s,k,wavecopy('v',c,s,k) / q(qi + 1));
   c = wavepaste('d',c,s,k,wavecopy('d',c,s,k) / q(qi + 2));
end
c = wavepaste('a',c,s,k,wavecopy('a',c,s,k) / q(qi + 3));
c = floor(c);       c = c .* sgn;

% Run-length code zero runs of more than 10. Begin by creating
% a special code for 0 runs ('zrc') and end-of-code ('eoc') and
% making a run-length table.
zrc = min(c(:)) - 1;      eoc = zrc - 1;

% The RUNS variable is shared with the RUNS variable in the nested
% function runcode. See Chapter 3 of DIPUM3E for a discussion of nested
% functions.
RUNS = 65535;

% Find the run transition points: 'plus' contains the index of the
% start of a zero run; the corresponding 'minus' is its end + 1.
z = c == 0;                 z = z - [0 z(1:end - 1)];
plus = find(z == 1);        minus = find(z == -1);

% Remove any terminating zero run from 'c'.
if length(plus) ~= length(minus)
   c(plus(end):end) = [];      c = [c eoc];
end

% Remove all other zero runs (based on 'plus' and 'minus') from 'c'.
for i = length(minus):-1:1
   run = minus(i) - plus(i);
   if run > 10
      ovrflo = floor(run / 65535);    run = run - ovrflo * 65535;
      c = [c(1:plus(i) - 1) repmat([zrc 1],1,ovrflo) zrc ...
         runcode(run) c(minus(i):end)];
   end
end

% Huffman encode and add misc. information for decoding.
y.runs    = uint16(RUNS);
y.s       = uint16(s(:));
y.zrc     = uint16(-zrc);
y.q       = uint16(100 * q');
y.n       = uint16(n);
y.huffman = mat2huff(c);

   %-------------------------------------------------------------------%
   function y = runcode(x)
      % Find a zero run in the run-length table. If not found, create a
      % new entry in the table. Return the index of the run.

      y = find(RUNS == x);
      if length(y) ~= 1
         RUNS = [RUNS; x];
         y = length(RUNS);
      end
   end
end

%----------------------------------------------------------------------%
function q = stepsize(n,p)
% Create a subband quantization array of step sizes ordered by
% decomposition (first to last) and subband (horizontal, vertical,
% diagonal, and for final decomposition the approximation subband).

if length(p) == 2               % Implicit Quantization
   q = [];
   qn = 2 ^ (8 - p(2) + n) * (1 + p(1) / 2 ^ 11);
   for k = 1:n
      qk = 2 ^ -k * qn;
      q = [q (2 * qk) (2 * qk) (4 * qk)];
   end
   q = [q qk];
else                            % Explicit Quantization
   q = p;
end

q = round(q * 100) / 100;       % Round to 1/100th place
if any(100 * q > 65535)
   error('The quantizing steps are not UINT16 representable.');
end
if any(q == 0)
   error('A quantizing step of 0 is not allowed.');
end
end

author: dipum

repo: dipum-toolbox

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


text:
function [results] = test_all_complements(models,final,grid,friend_classes,savegtids,saverecs)
VOCinit;

keep_thresh = 0.1;
max_os = 1.0;

for exid = 1:length(models)
  fprintf(1,'!');
  Iex = im2double(imread(sprintf(VOCopts.imgpath, ...
                                 models{exid}.curid)));
  models{exid}.sizeI = size(Iex);
end

%compute all friends

% classes = VOCopts.classes;

% countvec = zeros(1,20);
% for i = 1:length(models)
%   fprintf(1,'.');
%   recs = PASreadrecord(sprintf(VOCopts.annopath,models{i}.curid));
%   bbs = cat(1,recs.objects.bbox);
%   osmat = getosmatrix_bb(bbs(models{i}.objectid,:), bbs);
%   c = {recs.objects.class};
%   goodname = ismember(c,friend_classes);
%   targets = find(osmat>keep_thresh & osmat <=max_os & goodname);

%   if length(targets) > 1
%     fprintf(1,'got multiple d, keeping best\n',length(targets));
%     [aa,bb] = sort(osmat(targets),'descend');
%     targets = targets(bb);
%     targets = targets(bb(1));
%   end

%   targetc = c(targets);

%   models{i}.friendbb = bbs(targets,:);
%   models{i}.friendc = targetc;

%   [aa,bb] = ismember(targetc,classes);
%   countvec(bb) = countvec(bb)+1;
% end

% frac_covered = countvec ./ length(models)*100;
%models_with_friends = models;

%friend_classes = find(countvec > 10);
%friend_classes = classes(friend_classes);

for z = 1:length(friend_classes)
  targetcls = friend_classes{z};

  models = cellfun2(@(x)initialize_with_friends(x,targetcls), ...
                    models);

  % for aaa = 1:length(models_with_friends)
  %   oks = find(ismember(models_with_friends{z}.friendc,targetcls));
  %   models{z}.friendbb = models_with_friends{z}.friendbb(oks,:);
  %   models{z}.friendc = models_with_friends{z}.friendc(oks);
  % end

bboxes = final.unclipped_boxes;

minoverlap = .5;
if 1 %~exist('minoverlap','var')
  fprintf(1,'transferring friends\n');
  %oldbboxes = bboxes;
  bboxes = transfer_friends(models, bboxes);
end

if 0
maxdet = cellfun2(@(x)max(x(:,end)),bboxes);
lens = cellfun(@(x)length(x),maxdet);
goods = find(lens>0);
[alpha,beta] = sort([maxdet{goods}],'descend');
for i = 1:length(beta)
  curid = grid{goods(beta(i))}.curid;
  Ibase = imread(sprintf(VOCopts.imgpath,curid));
  figure(1)
  clf
  imagesc(Ibase)
  plot_bbox(bboxes{goods(beta(i))});
  pause
end
end

for i = 1:length(bboxes)
  bboxes{i} = clip_to_image(bboxes{i},grid{i}.imbb);
end

filer=sprintf(VOCopts.detrespath,'comp3',targetcls);
fprintf(1,'Writing File %s\n',filer);
fid = fopen(filer,'w');

for i = 1:length(bboxes)
  curid = grid{i}.curid;
  for q = 1:size(bboxes{i},1)
    fprintf(fid,'%s %f %f %f %f %f\n',curid,...
            bboxes{i}(q,end), bboxes{i}(q,1:4));
  end
end
fclose(fid);

if exist('minoverlap','var')
  VOCopts.minoverlap = minoverlap;
end

[recall,prec,ap,apold] = VOCevaldet(VOCopts,'comp3',targetcls, ...
                                    true,savegtids,saverecs);
results{z}.recall = recall;
results{z}.prc = prec;
results{z}.ap = ap;
results{z}.apold = apold;
results{z}.targetcls = targetcls;

figure;
plot(recall,prec)
title(sprintf('%s to predict %s: apold=%.3f',models{1}.cls,targetcls,apold));
targetcls
drawnow
end

author: quantombone

repo: exemplarsvm


text:
function r=randSequence(maxN,N)
% r=randSequence(maxN,N)
% produces a random  increasing sequence
% that includes N integers 1:maxN
% (N<maxN)

if N>=maxN, error('N is not less than maxN!!!'); end;

fact=maxN/((maxN-N)*2);
ind=[];

while length(ind)~=N,
  ind=find(round(rand(1,maxN)*fact)>0);
end;
r=1:maxN;
r=r(ind);

repo: CanlabCore


text:
1     % GMRES/Bicgstab(2)
1e-6  % tolerance
100   % maximum number of iterations
3     % lst sqrs commutator preconditioner
1     % ideal
0     % Print to (existing) Figure
5     %                     4
4     % Color

% Data file for it_solve for NS2

author: johannesgerer

repo: jburkardt-m


text:
%gampdf(2,3,4)
p = 0.018954;
q = wishpdf(2,3,1/4);
if abs(q - p) > 1e-6
	error('wishpdf failed');
end

x = [2 1; 1 3];
a = 5;
p = 0.0038062;
q = wishpdf(x,a);
if abs(q - p) > 1e-6
	error('wishpdf failed');
end

author: tminka

repo: lightspeed

save_path: github-repos/MATLAB/tminka-lightspeed


text:
function f=shah(L,a);
%SHAH  Discrete Shah-distribution
%   Usage: f=shah(L,a);
%
%   `shah(L,a)` computes the discrete, normalized Shah-distribution of
%   length *L* with a distance of *a* between the spikes.
%
%   The Shah distribution is defined by
%
%   .. f(n*a+1)=1/sqrt(L/a)
%
%   .. math:: f(n\cdot a+1)=\frac{1}{\sqrt(L/a)}
%
%   for integer *n*, otherwise *f* is zero.
%
%   This is also known as an impulse train or as the comb function, because
%   the shape of the function resembles a comb. It is the sum of unit
%   impulses ('diracs') with the distance *a*.
%
%   If *a* divides *L*, then the |dft| of `shah(L,a)` is `shah(L,L/a)`.
%
%   The Shah function has an extremely bad time-frequency localization.
%   It does not generate a Gabor frame for any *L* and *a*.
%
%   Examples:
%   ---------
%
%   A simple spectrogram of the Shah function (includes the negative
%   frequencies to display the whole TF-plane):::
%
%     sgram(shah(256,16),'dynrange',80,'nf')
%

%   AUTHOR : Peter L. Søndergaard
%   TESTING: OK
%   REFERENCE: OK

if nargin~=2
  error('Wrong number of input parameters.');
end;

%if mod(L,a)~=0
%  error('a must divide L.');
%end;

f=zeros(L,1);

f(1:a:L)=1/sqrt(L/a);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [divI, I_dx, I_dy] = computeDivergence(I_gx, I_gy)
%
%
%       [divI, I_dx, I_dy] = computeDivergence(I_gx, I_gy)
%
%       Input:
%           -I_gx: an input image
%           -I_gx
%
%       Output:
%           -I_dx:
%           -I_dy:
%           -divI:
%
%     Copyright (C) 2011-15  Francesco Banterle
%
%     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.
%
%     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, see <http://www.gnu.org/licenses/>.
%

if(~exist('I_gx', 'var') || ~exist('I_gy', 'var'))
    error('Gradients are needed to compute divergence');
end

kernelX = [0 0 0; -1 1 0; 0  0 0];
kernelY = [0 0 0;  0 1 0; 0 -1 0];

I_dx = imfilter(I_gx, kernelX, 'same');
I_dy = imfilter(I_gy, kernelY, 'same');

divI = RemoveSpecials(I_dx + I_dy);

end

author: banterle