text:
% Construct various DBNs and examine their clique structure.
% This was used to generate various figures in chap 3-4 of my thesis.

% Examine the cliques in the unrolled mildew net

%dbn = mk_mildew_dbn;
dbn = mk_chmm(4);
ss = dbn.nnodes_per_slice;
T = 7;
N = ss*T;
bnet = dbn_to_bnet(dbn, T);

constrained = 0;
if constrained
  stages = num2cell(unroll_set(1:ss, ss, T), 1);
else
  stages = { 1:N; };
end
clusters = {};
%[jtree, root, cliques, B, w, elim_order, moral_edges, fill_in_edges] = ...
%    dag_to_jtree(bnet, bnet.observed, stages, clusters);
[jtree, root, cliques] =  graph_to_jtree(moralize(bnet.dag), ones(1,N), stages, clusters);

flip=1;
clf;[dummyx, dummyy, h] = draw_dbn(dbn.intra, dbn.inter, flip, T, -1);
dir = '/home/eecs/murphyk/WP/Thesis/Figures/Inf/MildewUnrolled';
mk_ps_from_clqs(dbn, T, cliques, [])
%mk_collage_from_clqs(dir, cliques)

% Examine the cliques in the cascade DBN

% A-A
%  \
% B B
%  \
% C C
%  \
% D D
ss = 4;
intra = zeros(ss);
inter = zeros(ss);
inter(1, [1 2])=1;
for i=2:ss-1
  inter(i,i+1)=1;
end

% 2 coupled HMMs 1,3  and 2,4
ss = 4;
intra = zeros(ss);
inter = zeros(ss); % no persistent edges
%inter = diag(ones(ss,1)); % persitence edges
inter(1,3)=1; inter(3,1)=1;
inter(2,4)=1; inter(4,2)=1;

%bnet = mk_fhmm(3);
bnet = mk_chmm(4);
intra = bnet.intra;
inter = bnet.inter;

clqs = compute_minimal_interface(intra, inter);
celldisp(clqs)

% A A
%  \
% B B
%  \
% C C
%  \
% D-D
ss = 4;
intra = zeros(ss);
inter = zeros(ss);
for i=1:ss-1
  inter(i,i+1)=1;
end
inter(4,4)=1;

ns = 2*ones(1,ss);
dbn = mk_dbn(intra, inter, ns);
for i=2*ss
  dbn.CPD{i} = tabular_CPD(bnet, i);
end

T = 4;
N = ss*T;
bnet = dbn_to_bnet(dbn, T);

constrained = 1;
if constrained
  % elim first 3 slices first in any order
  stages = {1:12, 13:16};
  %stages = num2cell(unroll_set(1:ss, ss, T), 1);
else
  stages = { 1:N; };
end
clusters = {};
%[jtree, root, cliques, B, w, elim_order, moral_edges, fill_in_edges] = ...
%    dag_to_jtree(bnet, bnet.observed, stages, clusters);
[jtree, root, cliques] =  graph_to_jtree(moralize(bnet.dag), ones(1,N), stages, clusters);

% Examine the cliques in the 1.5 slice DBN

%dbn = mk_mildew_dbn;
dbn = mk_water_dbn;
%dbn = mk_bat_dbn;
ss = dbn.nnodes_per_slice;
int = compute_fwd_interface(dbn);
bnet15 = mk_slice_and_half_dbn(dbn, int);
N = length(bnet15.dag);
stages = {1:N};

% bat
%cl1 = [16 17 19 7 14];
%cl2 = [27 25 21 23 20];
%clusters = {cl1, cl2, cl1+ss, cl2+ss};

% water
%cl1 = 1:2; cl2 = 3:6; cl3 = 7:8;
%clusters = {cl1, cl2, cl3, cl1+ss, cl2+ss, cl3+ss};

%clusters = {};
clusters = {int, int+ss};
%[jtree, root, cliques, B, w, elim_order, moral_edges, fill_in_edges] = ...
%    dag_to_jtree(bnet15, bnet.observed, stages, clusters);
[jtree, root, cliques] =  graph_to_jtree(moralize(bnet15.dag), ones(1,N), stages, clusters);

clq_len = [];
for c=1:length(cliques)
  clq_len(c) = length(cliques{c});
end
hist(clq_len, 1:max(clq_len));
h=hist(clq_len, 1:max(clq_len));
axis([1 max(clq_len)+1 0 max(h)+1])
xlabel('clique size','fontsize',16)
ylabel('number','fontsize',16)

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [C, sigma] = dataset3Params(X, y, Xval, yval)
%EX6PARAMS returns your choice of C and sigma for Part 3 of the exercise
%where you select the optimal (C, sigma) learning parameters to use for SVM
%with RBF kernel
%   [C, sigma] = EX6PARAMS(X, y, Xval, yval) returns your choice of C and
%   sigma. You should complete this function to return the optimal C and
%   sigma based on a cross-validation set.
%

% You need to return the following variables correctly.
C = 1;
sigma = 0.1;

param = [0.01 , 0.03, 0.1, 0.3, 1, 3, 10, 30];

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the optimal C and sigma
%               learning parameters found using the cross validation set.
%               You can use svmPredict to predict the labels on the cross
%               validation set. For example,
%                   predictions = svmPredict(model, Xval);
%               will return the predictions on the cross validation set.
%
%  Note: You can compute the prediction error using
%        mean(double(predictions ~= yval))
%

minError = 10000.0;

%for CVal = param,
%	for sigmaVal = param,
%		model =  svmTrain(X, y, CVal, @(x1, x2) gaussianKernel(x1, x2, sigmaVal));
%		predictions = svmPredict(model , Xval);
%		error = mean(double(predictions ~= yval));
%		if minError > error,
%			minError = error;
%			C = CVal;
%			sigma = sigmaVal;
%		end
%	end
%end

% =========================================================================

end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
function Z = ge(X,Y)
%GE Greater than or equal (>=) for tensors.
%
%   See also TENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

Z = tenfun(@ge,X,Y);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
classdef prtRvMvb < prtRv & prtRvMemebershipModel
    % prtRvMvv  Multivariate binary random variable
    %
    %   RV = prtRvMvv creates a prtRvMvv object with empty probability
    %   vector. The probability vector must be set either directly, or by
    %   calling the MLE method.
    %
    %   RV = prtRvMvn(PROPERTY1, VALUE1,...) creates a prtRvMvb object RV
    %   with properties as specified by PROPERTY/VALUE pairs.
    %
    %   A prtRvMvb object inherits all properties from the prtRv class. In
    %   addition, it has the following properties:
    %
    %   probabilities - 1 x nDimensions vector with values between 0 and 1
    %                   specifying the probability of each dimension taking
    %                   the value true.
    %
    %  Example:
    %
    %   See also: prtRv, prtRvGmm, prtRvMultinomial, prtRvUniform,
    %   prtRvUniformImproper, prtRvVq, prtRvDiscrete

    properties (SetAccess = private)
        name = 'Multi-Variate Binary';
        nameAbbreviation = 'RVMVB';
    end

    properties (SetAccess = protected)
        isSupervised = false;
        isCrossValidateValid = true;
    end

    properties (Dependent)
        probabilities
    end
    properties (Dependent, SetAccess = 'private')
        logProbabilities
        logOneMinusProbabilities
    end
    properties (Hidden, Dependent)
        minProbability
        maxProbability
    end
    properties (Hidden, SetAccess='private', GetAccess='private')
        minProbabilityDepHelper = 0;
        maxProbabilityDepHelper = 1;
        probabilitiesDepHelper
    end

    properties (Hidden = true, Dependent = true)
        nDimensions
    end

    methods
        function R = prtRvMvb(varargin)
            R.name = 'Multi-Variate Binary';
            R = constructorInputParse(R,varargin{:});
        end

        function R = mle(R,X)
           % MLE Compute the maximum likelihood estimate
            %
            % RV = RV.mle(X) computes the maximum likelihood estimate based
            % the data X. X should be nObservations x nDimensions.

            X = R.dataInputParse(X); % Basic error checking etc

            R = weightedMle(R,X,ones(size(X,1),1));
        end

        function vals = pdf(R,X)
            % PDF Output the pdf of the random variable evaluated at the points specified
            %
            % pdf = RV.pdf(X) returns  the pdf of the prtRv
            % object evaluated at X. X must be an N x nDims matrix, where
            % N is the number of locations to evaluate the pdf, and nDims
            % is the same as the number of dimensions, nDimensions, of the
            % prtRv object RV.

            [isValid, reasonStr] = R.isValid;
            assert(isValid,'PDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            X = R.dataInputParse(X); % Basic error checking etc
            assert(size(X,2) == R.nDimensions,'Data, RV dimensionality missmatch. Input data, X, has dimensionality %d and this RV has dimensionality %d.', size(X,2), R.nDimensions)
            assert((isnumeric(X) || islogical(X)) && ndims(X)==2,'X must be a 2D numeric or logical array.');

            vals = prod(bsxfun(@times,X,R.probabilities) + bsxfun(@times,~X,1-R.probabilities),2);
        end

        function vals = logPdf(R,X)
            % LOGPDF Output the log pdf of the random variable evaluated at the points specified
            %
            % logpdf = RV.logpdf(X) returns the logarithm of value of the
            % pdf of the prtRv object evaluated at X. X must be an N x
            % nDims matrix, where N is the number of locations to evaluate
            % the pdf, and nDims is the same as the number of dimensions,
            % nDimensions, of the prtRv object RV.

            [isValid, reasonStr] = R.isValid;
            assert(isValid,'LOGPDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            X = R.dataInputParse(X); % Basic error checking etc
            assert(size(X,2) == R.nDimensions,'Data, RV dimensionality missmatch. Input data, X, has dimensionality %d and this RV has dimensionality %d.', size(X,2), R.nDimensions)
            assert((isnumeric(X) || islogical(X)) && ndims(X)==2,'X must be a 2D numeric or logical array.');

            term1 = bsxfun(@times,X,R.logProbabilities); term1(isnan(term1)) = 0; %false*-inf = NaN;
            term2 = bsxfun(@times,~X,R.logOneMinusProbabilities); term2(isnan(term2)) = 0; %false*-inf = NaN;
            vals = sum(term1 + term2,2);
        end

        function vals = draw(R,N)
            % DRAW  Draw random samples from the distribution described by the prtRv object
            %
            % VAL = RV.draw(N) generates N random samples drawn from the
            % distribution described by the prtRv object RV. VAL will be a
            % N x nDimensions vector, where nDimensions is the number of
            % dimensions of RV.

            assert(numel(N)==1 && N==floor(N) && N > 0,'N must be a positive integer scalar.')
            [isValid, reasonStr] = R.isValid;
            assert(isValid,'DRAW cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            vals = bsxfun(@ge , rand(N,R.nDimensions), R.probabilities);
        end
    end

    methods (Hidden=true)
        function [val, reasonStr] = isValid(R)

            if numel(R) > 1
                val = false(size(R));
                for iR = 1:numel(R)
                    [val(iR), reasonStr] = isValid(R(iR));
                end
                return
            end

            val = ~isempty(R.probabilities);

            if val
                reasonStr = '';
            else
                reasonStr = 'because probabilities has not been set';
            end
        end

        function R = weightedMle(R,X,weights)

            assert(numel(weights)==size(X,1),'The number of weights must mach the number of observations.');

            if ~islogical(X)
                assert(isnumeric(X) && all(all(X >= 0 & X<=1)),'prt:prtRvMvb:mle','input data to prtRvMvb.mle must contain values between 0 and 1')
            end

            weights = weights(:);

            R.probabilities = sum(bsxfun(@times,X,weights))/sum(weights);

        end

        function initMembershipMat = initializeMixtureMembership(Rs,X)

            learningInitialMembershipFactor = 0.9;

            [classMeans,kmMembership] = prtUtilKmeans(X,length(Rs),'handleEmptyClusters','random','distanceMetricFn',@prtDistanceHamming,'maxIterations',100,'logicalMeans',true); %#ok<ASGLU>

            initMembershipMat = zeros(size(X,1),length(Rs));
            for iComp = 1:length(Rs)
                initMembershipMat(kmMembership == iComp, iComp) = learningInitialMembershipFactor;
            end
            initMembershipMat(initMembershipMat==0) = (1-learningInitialMembershipFactor)./(length(Rs)-1);

            % We should normalize this just in case the
            % learningInitialMembershipFactor was set poorly
            initMembershipMat = bsxfun(@rdivide,initMembershipMat,sum(initMembershipMat,2));
        end
    end

    % Get methods
    methods
        function val = get.nDimensions(R)
            val = getNumDimensions(R);
        end
        function val = get.probabilities(R)
            val = R.probabilitiesDepHelper;
        end
        function val = get.minProbability(R)
            val = R.minProbabilityDepHelper;
        end
        function val = get.maxProbability(R)
            val = R.maxProbabilityDepHelper;
        end
        function val = get.logProbabilities(R)
            val = log(R.probabilities);
        end
        function val = get.logOneMinusProbabilities(R)
            val = log(1 - R.probabilities);
        end
    end

    methods (Access = 'protected')
        function val = getNumDimensions(R)
            if ~isempty(R.probabilities)
                val = length(R.probabilities);
            else
                val = [];
            end
        end
    end
    % Set Methods
    methods
        function R = set.probabilities(R,val)
            assert((isnumeric(val) || islogical(val)) && isvector(val),'probabilities must be a vector of values between 0 and 1 (inclusively');
            assert(all(val<=1 & val>=0),'prt:prtRvMvb:badProbabilities','probabilities must have values between 0 and 1 (inclusively)')

            val = val(:)';

            val(val>R.maxProbability) = R.maxProbability;
            val(val<R.minProbability) = R.minProbability;

            R.probabilitiesDepHelper = val;
        end

        function R = set.minProbability(R,val)
            assert(prtUtilIsPositiveScalar(val) && val<=1 && val>=0,'minProbability must be a value between 0 and 1 (inclusively)');
            R.minProbabilityDepHelper = val;
        end
        function R = set.maxProbability(R,val)
            assert(prtUtilIsPositiveScalar(val) && val<=1 && val>=0,'maxProbability must be a value between 0 and 1 (inclusively)');
            R.maxProbabilityDepHelper = val;
        end

    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function P = calcP(j)
global uLINK

if j == 0
   P = [0 0 0]';
else
   c1 = uLINK(j).R * uLINK(j).c;
   P = uLINK(j).m * (uLINK(j).v + cross(uLINK(j).w, c1) );
   P = P + calcP(uLINK(j).sister) + calcP(uLINK(j).child);
end

author: s-kajita

repo: IntroductionToHumanoidRobotics


text:
function varargout = panel_eeglab_cond(varargin)
% PANEL_EEGLAB_COND: Selection of conditions from EegLab .set files.
%
% USAGE:  bstPanelNew = panel_eeglab_cond('CreatePanel', panelName)
%           selParams = panel_eeglab_cond('GetSelectedParameters')
%        selectedCond = panel_eeglab_cond('GetPanelContents')

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

eval(macro_method);
end

%% ===== CREATE PANEL =====
function [bstPanelNew, panelName] = CreatePanel(paramNames, paramValues) %#ok<DEFNU>
    % Java initializations
    import java.awt.*;
    import javax.swing.*;
    import javax.swing.table.*;
    import org.brainstorm.table.*;
    panelName = 'EeglabConditions';

    % Create panel
    jPanelNew = gui_component('Panel');

    % ===== PARAMETERS LIST =====
    jPanelParams = gui_component('Panel');
    jBorder = java_scaled('titledborder', 'Parameters');
    jPanelParams.setBorder(jBorder);

    jPanelParams.setPreferredSize(java_scaled('dimension', 110, 200));
        % JList
        jListParams = JList(paramNames);
        jListParams.setFont(bst_get('Font'));
        java_setcb(jListParams, 'ValueChangedCallback', @ParamSelectionChanged_Callback);
        % Create scroll panel
        jScrollPanelParams = JScrollPane(jListParams);
        jScrollPanelParams.setBorder([]);
        jPanelParams.add(jScrollPanelParams, BorderLayout.CENTER);
    jPanelNew.add(jPanelParams, BorderLayout.WEST);

    % ===== CONDITIONS NAMES =====
    jPanelCond = gui_component('Panel');
    jBorder = java_scaled('titledborder', 'Conditions');
    jPanelCond.setBorder(jBorder);
    jPanelCond.setPreferredSize(java_scaled('dimension', 400, 200));
        % JTable
        jTableCond = JTable(DefaultTableModel({'Parameters', 'Condition name'}, 0));
        jTableCond.setFont(bst_get('Font'));
        jTableCond.getColumnModel.getColumn(0).setCellEditor(DisabledCellEditor());
        % Create scroll panel
        jScrollPanelCond = JScrollPane(jTableCond);
        jScrollPanelCond.setBorder([]);
        jPanelCond.add(jScrollPanelCond, BorderLayout.CENTER);
    jPanelNew.add(jPanelCond, BorderLayout.CENTER);

    % ===== VALIDATION BUTTONS =====
    jPanelValidation = gui_river([5,5], [0,5,12,5]);
        % Cancel
        jButtonCancel = JButton('Cancel');
        java_setcb(jButtonCancel, 'ActionPerformedCallback', @ButtonCancel_Callback);
        jPanelValidation.add('br right', jButtonCancel);
        % Ok
        jButtonOk = JButton('OK');
        java_setcb(jButtonOk, 'ActionPerformedCallback', @ButtonOk_Callback);
        jPanelValidation.add(jButtonOk);
    jPanelNew.add(jPanelValidation, BorderLayout.SOUTH);

    % Return a mutex to wait for panel close
    bst_mutex('create', panelName);

    % Create the BstPanel object that is returned by the function
    % => constructor BstPanel(jHandle, panelName, sControls)
    bstPanelNew = BstPanel(panelName, ...
                           jPanelNew, ...
                           struct('jListParams', jListParams, ...
                                  'jTableCond',  jTableCond));

%% =================================================================================
%  === LOCAL CALLBACKS =============================================================
%  =================================================================================
%% ===== CANCEL =====
    function ButtonCancel_Callback(varargin)
        % Close panel
        gui_hide(panelName);
    end

%% ===== OK =====
    function ButtonOk_Callback(varargin)
        % Close cell editor (validate last edited string)
        cellEditor = jTableCond.getCellEditor();
        if ~isempty(cellEditor)
        	cellEditor.stopCellEditing();
        end
        % Release mutex and keep the panel opened
        bst_mutex('release', panelName);
    end

%% ===== SELECTION CHANGED =====
    function ParamSelectionChanged_Callback(hObject, ev)
        global combinCond;
        if ~ev.getValueIsAdjusting()
            % Get selected items
            iSel = double(jListParams.getSelectedIndices());
            % Update JList conditions list
            tableModel = jTableCond.getModel();
            tableModel.setRowCount(0);
            % Check number of final combinations
            nbCond = prod(cellfun(@length, paramValues(iSel+1)));
%             % Too many combinations
%             if (nbCond > 500)
%                 combinCond = [];
%                 tableModel.addRow({'(Error: too many combinations)', ''});
%             else
                % Compute a structure that describes all the combinations
                combinCond = GetConditionCombinations(paramNames(iSel+1), paramValues(iSel+1));
                % Process all the combinations
                params = paramNames(iSel+1);
                for iCond = 1:length(combinCond)
                    strCondName = '';
                    strParamList = '';
                    for iParam = 1:length(params)
                        if (iParam ~= 1)
                            strCondName = [strCondName '_'];
                            strParamList = [strParamList, ', '];
                        end
                        if ischar(combinCond(iCond).(params{iParam}))
                            tmpStr = strrep(file_standardize(combinCond(iCond).(params{iParam})), '_', '-');
                            strCondName  = [strCondName, params{iParam}, tmpStr];
                            strParamList = [strParamList, params{iParam}, '=', combinCond(iCond).(params{iParam})];
                        else
                            strCondName  = [strCondName, params{iParam}, sprintf('%d', combinCond(iCond).(params{iParam}))];
                            strParamList = [strParamList, params{iParam}, sprintf('=%d', combinCond(iCond).(params{iParam}))];
                        end
                    end
                    tableModel.addRow({strParamList, strCondName});
                end
%             end
        end
    end
end

%% =================================================================================
%  === EXTERNAL CALLBACKS ==========================================================
%  =================================================================================
%% ===== GET SELECTED PARAMETERS =====
function selParams = GetSelectedParameters() %#ok<DEFNU>
    ctrl = bst_get('PanelControls', 'EeglabConditions');
    selObj = ctrl.jListParams.getSelectedValues();
    selParams = cell(1, length(selObj));
    for i = 1:length(selObj)
        selParams{i} = char(selObj(i));
    end
end

%% ===== GET PANEL RESULTS =====
function selectedCond = GetPanelContents(varargin) %#ok<DEFNU>
    % Get panel controls
    ctrl = bst_get('PanelControls', 'EeglabConditions');
    if isempty(ctrl)
        selectedCond = [];
        return
    end
    % Get current combination description
    global combinCond;
    selectedCond = combinCond;
    combinCond = [];
    clear global combinCond;
    % If nothing is available
    if (length(selectedCond) < 2)
        selectedCond = [];
    else
        % Get selected indices
        iCondSel = ctrl.jTableCond.getSelectedRows();
        if ~isempty(iCondSel)
            iCondSel = double(iCondSel(:))' + 1;
        else
            iCondSel = 1:length(selectedCond);
        end
        % Read new conditions names
        for iCond = iCondSel
            selectedCond(iCond).Name = file_standardize(char(ctrl.jTableCond.getModel.getValueAt(iCond-1, 1)));
        end
        selectedCond = selectedCond(iCondSel);
    end
end

%% =================================================================================
%  === HELPERS =====================================================================
%  =================================================================================
function combinCond = GetConditionCombinations(paramName, paramValues)
    % Get final number of conditions
    nbValues = cellfun(@length, paramValues);
    nbCond = prod(nbValues);
    % Build structure returned for each conditions
    structCond = struct('Name', '');
    for iField = 1:length(paramName)
        structCond.(paramName{iField}) = [];
    end
    combinCond = repmat(structCond, 1, nbCond);
    % Process each condition
    for iCond = 1:nbCond
        iValues = my_ind2sub(nbValues, iCond);
        % Process each parameter
        for iParam = 1:length(paramValues)
            % Get the value for parameter
            if iscell(paramValues{iParam})
                combinCond(iCond).(paramName{iParam}) = paramValues{iParam}{(iValues(iParam))};
            else
                combinCond(iCond).(paramName{iParam}) = paramValues{iParam}((iValues(iParam)));
            end
        end
    end
end

%% ===== CONVERT INDICES TO SUBINDICES =====
function loc = my_ind2sub(siz, ndx)
    n = length(siz);
    k = [1 cumprod(siz(1:end-1))];
    loc = zeros(size(siz));
    for i = n:-1:1,
      vi = rem(ndx-1, k(i)) + 1;
      vj = (ndx - vi)/k(i) + 1;
      loc(i) = vj;
      ndx = vi;
    end
end

author: brainstorm-tools


text:
%% Sample execution for Baseline-conv5 (improved Siam-FC)
% hyper-parameters reported in Supp.material for CVPR'17, Table 2 for arXiv version
tracker_par.join.method = 'xcorr';
tracker_par.net = 'baseline-conv5_e55.mat';
tracker_par.net_gray = 'baseline-conv5_gray_e100.mat';
tracker_par.scaleStep = 1.0470;
tracker_par.scalePenalty = 0.9825;
tracker_par.scaleLR = 0.68;
tracker_par.wInfluence = 0.175;
tracker_par.zLR = 0.0102;

[~,~,dist,overlap,~,~,~,~] = run_tracker_evaluation('all', tracker_par);

author: bertinetto

repo: cfnet

save_path: github-repos/MATLAB/bertinetto-cfnet


text:
function ierror = ksubset_colex_check ( k, n, t )

%*****************************************************************************80
%
%% KSUBSET_COLEX_CHECK checks a K subset in colex form.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 January 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Donald Kreher, Douglas Simpson,
%    Combinatorial Algorithms,
%    CRC Press, 1998,
%    ISBN: 0-8493-3988-X,
%    LC: QA164.K73.
%
%  Parameters:
%
%    Input, integer K, the number of elements each K subset must
%    have. 1 <= K <= N.
%
%    Input, integer N, the number of elements in the master set.
%    N must be positive.
%
%    Input, integer T(K), describes a K subset.  T(I) is the I-th
%    element of the K subset.  The elements must be listed in
%    DESCENDING order.
%
%    Output, integer IERROR, error flag.
%    0, no error.
%    -1, N is not positive.
%    -2, K is not positive.
%    I, entry I is illegal.
%
  ierror = 0;

  if ( n < 1 )
    ierror = -1;
    return
  end

  if ( k < 1 || n < k )
    ierror = -2;
    return
  end

  tmax = n + 1;

  for i = 1 : k

    if ( t(i) <= 0 || tmax <= t(i) )
      ierror = i;
      return
    end

    tmax = t(i);

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% LOGNORMAL_RAND
%
%   X = LOGNORMAL_RAND(MU, SIGMA)
%
% MU is the location parameter in (-INF, INF)
% SIGMA is the scale parameter in (0, INF)
%
% Function calls
%   X = LOGNORMAL_RAND(MU,SIGMA,M,N,...)
%   X = LOGNORMAL_RAND(MU,SIGMA,[M,N,...])
% return an M-by-N-by-... array.
%
% See also LOGNORMAL_LOGPDF, NORMRND.

% Last modified 2010-11-12
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function x = lognormal_rand(varargin)

x = exp(normrnd(varargin{:}));

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [varargout]=quiverVec(varargin)

% function [h]=quiverVec(P,V,vecSize,colorSpec,edgeColorOpt,quiverStyleOpt,alphaLevel)
% ------------------------------------------------------------------------
% Plots the vectors V with origins P as patch data arrows. The vectors size
% can be specified as vecSize allong with the colors in colorSpec. Other
% options include the edge color, quiver style, and transparency level.
%
%
% Kevin Mattheus Moerman
%
% Change log:
% 2020/07/14 Added single origin support
% ------------------------------------------------------------------------

%% Parse input

switch nargin
    case 2
        P=varargin{1};
        V=varargin{2};
        vecSize=[];
        colorSpec=[];
        edgeColorOpt='none';
        quiverStyleOpt=1;
        alphaLevel=1;
    case 3
        P=varargin{1};
        V=varargin{2};
        vecSize=varargin{3};
        colorSpec=[];
        edgeColorOpt='none';
        quiverStyleOpt=1;
        alphaLevel=1;
    case 4
        P=varargin{1};
        V=varargin{2};
        vecSize=varargin{3};
        colorSpec=varargin{4};
        edgeColorOpt='none';
        quiverStyleOpt=1;
        alphaLevel=1;
    case 5
        P=varargin{1};
        V=varargin{2};
        vecSize=varargin{3};
        colorSpec=varargin{4};
        edgeColorOpt=varargin{5};
        quiverStyleOpt=1;
        alphaLevel=1;
    case 6
        P=varargin{1};
        V=varargin{2};
        vecSize=varargin{3};
        colorSpec=varargin{4};
        edgeColorOpt=varargin{5};
        quiverStyleOpt=varargin{6};
        alphaLevel=1;
    case 7
        P=varargin{1};
        V=varargin{2};
        vecSize=varargin{3};
        colorSpec=varargin{4};
        edgeColorOpt=varargin{5};
        quiverStyleOpt=varargin{6};
        alphaLevel=varargin{7};
end

if size(P,1)==1 %Expand to match V if single origin is given
    P=P(ones(size(V,1),1),:);
end

%%
if ~isempty(P) || ~isempty(V)

    if isempty(edgeColorOpt)
        edgeColorOpt='none';
    end

    if isempty(quiverStyleOpt)
        quiverStyleOpt=1;
    end

    switch quiverStyleOpt
        case 1 %Depart from origin
            %Keep as is
        case 2 %Arrive at origin
            P=P-V;
        case 3 %Pass through origin
            P=P-(V/2);
        case 4 %Two-sided
            P=[P;P];
            V=[V;-V];
            if ~ischar(colorSpec) && size(colorSpec,1)>1
                colorSpec=[colorSpec;colorSpec];
            end
    end

    if size(P,2)==2
        P(:,3)=0;
    end

    if size(V,2)==2
        V(:,3)=0;
    end

    if numel(vecSize)==1
        vecSize=vecSize*ones(1,2);
    end

    if ischar(colorSpec)
        [F,P,~]=quiver3Dpatch(P(:,1),P(:,2),P(:,3),V(:,1),V(:,2),V(:,3),[],vecSize);
        C=colorSpec;
    else
        if size(colorSpec,1)==1 %If only 1 color is provided
            colorSpec=colorSpec(ones(size(P,1),1),:); %copy for all vectors
        end

        [F,P,C]=quiver3Dpatch(P(:,1),P(:,2),P(:,3),V(:,1),V(:,2),V(:,3),colorSpec,vecSize);
    end

    h=gpatch(F,P,C,edgeColorOpt,alphaLevel);
else
    h=[];
end

if nargout>0
    varargout{1}=h;
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


% -------------------------------------------------------------------------
%   Description:
%       Demo script to calculate the Kendall coefficient of agreement
%       This script reproduces the results of Figure 4 in our paper.
%
%   Citation:
%       A Comparative Study for Single Image Blind Deblurring
%       Wei-Sheng Lai, Jia-Bin Huang, Zhe Hu, Narendra Ahuja, and Ming-Hsuan Yang
%       IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2016
%
%   Contact:
%       Wei-Sheng Lai
%       wlai24@ucmerced.edu
%       University of California, Merced
% -------------------------------------------------------------------------

%% input dataset and attributes
dataset = 'real';
% dataset = 'uniform';
% dataset = 'nonuniform';

attribute = {};
attribute{end+1} = 'manmade';
attribute{end+1} = 'natural';
attribute{end+1} = 'people';
attribute{end+1} = 'saturated';
attribute{end+1} = 'text';
attribute{end+1} = 'all';

num_method = 14; % total number of evaluated methods

fprintf(' %10s | %10s | %8s\n', 'dataset', 'attribute', 'kendall');
fprintf('------------------------------------\n')

for a = 1:length(attribute)

    %% load votes
    vote_filename = fullfile('votes', sprintf('votes_%s_balance_%s.csv', dataset, attribute{a}));
    M = csvread(vote_filename, 1, 0); % offset the first row to skip header

    %% convert M to winning matrix
    C = construct_winning_matrix(M, num_method);

    %% compute kendall coefficient of agreement
    k = coefficient_of_agreement(C);
    fprintf(' %10s | %10s | %f\n', dataset, attribute{a}, k);

end


text:
classdef LSMOP8 < PROBLEM
% <multi/many> <real> <large/none>
% Large-scale benchmark MOP
% nk --- 5 --- Number of subcomponents in each variable group

%------------------------------- Reference --------------------------------
% R. Cheng, Y. Jin, and M. Olhofer, Test problems for large-scale
% multiobjective and many-objective optimization, IEEE Transactions on
% Cybernetics, 2017, 47(12): 4108-4121.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties(Access = private)
        nk = 5; % Number of subcomponents in each variable group
        sublen;	% Number of variables in each subcomponent
        len;    % Cumulative sum of lengths of variable groups
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            % Parameter setting
            obj.nk = obj.ParameterSet(5);
            if isempty(obj.M); obj.M = 3; end
            if isempty(obj.D); obj.D = 100*obj.M; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = [ones(1,obj.M-1),10.*ones(1,obj.D-obj.M+1)];
            obj.encoding = ones(1,obj.D);
            % Calculate the number of variables in each subcomponent
            c = 3.8*0.1*(1-0.1);
            for i = 1 : obj.M-1
                c = [c,3.8.*c(end).*(1-c(end))];
            end
            obj.sublen = floor(c./sum(c).*(obj.D-obj.M+1)/obj.nk);
            obj.len    = [0,cumsum(obj.sublen*obj.nk)];
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            [N,D] = size(PopDec);
            M     = obj.M;
            PopDec(:,M:D) = (1+repmat(cos((M:D)./D*pi/2),N,1)).*PopDec(:,M:D) - repmat(PopDec(:,1)*10,1,D-M+1);
            G = zeros(N,M);
            for i = 1 : 2 : M
                for j = 1 : obj.nk
                    G(:,i) = G(:,i) + Griewank(PopDec(:,obj.len(i)+M-1+(j-1)*obj.sublen(i)+1:obj.len(i)+M-1+j*obj.sublen(i)));
                end
            end
            for i = 2 : 2 : M
                for j = 1 : obj.nk
                    G(:,i) = G(:,i) + Sphere(PopDec(:,obj.len(i)+M-1+(j-1)*obj.sublen(i)+1:obj.len(i)+M-1+j*obj.sublen(i)));
                end
            end
            G      = G./repmat(obj.sublen,N,1)./obj.nk;
            PopObj = (1+G+[G(:,2:end),zeros(N,1)]).*fliplr(cumprod([ones(N,1),cos(PopDec(:,1:M-1)*pi/2)],2)).*[ones(N,1),sin(PopDec(:,M-1:-1:1)*pi/2)];
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,obj.M);
            R = R./repmat(sqrt(sum(R.^2,2)),1,obj.M);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                R = obj.GetOptimum(100);
            elseif obj.M == 3
                a = linspace(0,pi/2,10)';
                R = {sin(a)*cos(a'),sin(a)*sin(a'),cos(a)*ones(size(a'))};
            else
                R = [];
            end
        end
    end
end

function f = Griewank(x)
    f = sum(x.^2,2)./4000 - prod(cos(x./repmat(sqrt(1:size(x,2)),size(x,1),1)),2) + 1;
end

function f = Sphere(x)
    f = sum(x.^2,2);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function xd = lagrange_interp_nd_grid ( m, n_1d, a, b, nd )

%*****************************************************************************80
%
%% LAGRANGE_INTERP_ND_GRID sets an M-dimensional Lagrange interpolant.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N_1D(M), the order of the 1D rule to be used
%    in each dimension.
%
%    Input, real A(M), B(M), the lower and upper interval endpoints in
%    each dimension.
%
%    Input, integer ND, the number of points in the product grid.
%
%    Output, real XD(M,ND), the points at which data was sampled.
%
    xd = zeros ( m, nd );
    for i = 1 : m
      n = n_1d(i);
      x_1d(1:n) = cc_compute_points ( n );
      x_1d(1:n) = 0.5 * ( ( 1.0 - x_1d(1:n) ) * a(i) ...
                        + ( 1.0 + x_1d(1:n) ) * b(i) );
      xd = r8vec_direct_product ( i, n, x_1d, m, nd, xd );
    end

  return
end


text:
function toolboxUpdateHeader
% Update the headers of all the files.
%
% USAGE
%  toolboxUpdateHeader
%
% INPUTS
%
% OUTPUTS
%
% EXAMPLE
%
% See also
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.40
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

header={
  'Piotr''s Computer Vision Matlab Toolbox      Version 3.40'; ...
  'Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]'; ...
  'Licensed under the Simplified BSD License [see external/bsd.txt]'};

root=fileparts(fileparts(mfilename('fullpath')));
ds=dir(root); ds=ds([ds.isdir]); ds={ds.name};
ds=ds(3:end); ds=setdiff(ds,{'.git','doc'});
subds = { '/', '/private/' };
exts = {'m','c','cpp','h','hpp'};
omit = {'Contents.m','fibheap.h','fibheap.cpp'};

for i=1:length(ds)
  for j=1:length(subds)
    for k=1:length(exts)
      d=[root '/' ds{i} subds{j}];
      if(k==1), comment='%'; else comment='*'; end
      fs=dir([d '*.' exts{k}]); fs={fs.name}; fs=setdiff(fs,omit);
      n=length(fs); for f=1:n, fs{f}=[d fs{f}]; end
      for f=1:n, toolboxUpdateHeader1(fs{f},header,comment); end
    end
  end
end

end

function toolboxUpdateHeader1( fName, header, comment )

% set appropriate comment symbol in header
m=length(header); for i=1:m, header{i}=[comment ' ' header{i}]; end

% read in file and find header
disp(fName); lines=readFile(fName);
loc = find(not(cellfun('isempty',strfind(lines,header{1}(1:40)))));
if(isempty(loc)), error('NO HEADER: %s\n',fName); end; loc=loc(1);

% check that header is properly formed, return if up to date
for i=1:m; assert(isequal(lines{loc+i-1}(1:10),header{i}(1:10))); end
if(~any(strfind(lines{loc},'NEW'))); return; end

% update copyright year and overwrite rest of header
lines{loc+1}(13:16)=header{2}(13:16);
for i=[1 3:m]; lines{loc+i-1}=header{i}; end
writeFile( fName, lines );

end

function lines = readFile( fName )
fid = fopen( fName, 'rt' ); assert(fid~=-1);
lines=cell(10000,1); n=0;
while( 1 )
  n=n+1; lines{n}=fgetl(fid);
  if( ~ischar(lines{n}) ), break; end
end
fclose(fid); n=n-1; lines=lines(1:n);
end

function writeFile( fName, lines )
fid = fopen( fName, 'w' );
for i=1:length(lines); fprintf( fid, '%s\n', lines{i} ); end
fclose(fid);
end

author: MengyangPu

repo: EDTER

save_path: github-repos/MATLAB/MengyangPu-EDTER


text:
function c8lib_test ( )

%*****************************************************************************80
%
%% C8LIB_TEST tests the C8LIB library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 March 2015
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8LIB_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the C8LIB library.\n' );

  c8_abs_test ( );
  c8_acos_test ( );
  c8_acosh_test ( );
  c8_add_test ( );
  c8_arg_test ( );
  c8_asin_test ( );
  c8_asinh_test ( );
  c8_atan_test ( );
  c8_atanh_test ( );
  c8_conj_test ( );
  c8_cos_test ( );
  c8_cosh_test ( );
  c8_cube_root_test ( );
  c8_div_test ( );
  c8_div_r8_test ( );
  c8_exp_test ( );
  c8_i_test ( );
  c8_imag_test ( );
  c8_inv_test ( );
  c8_le_l1_test ( );
  c8_le_l2_test ( );
  c8_le_li_test ( );
  c8_log_test ( );
  c8_mag_test ( );
  c8_mul_test ( );
  c8_nint_test ( );
  c8_norm_l1_test ( );
  c8_norm_l2_test ( );
  c8_norm_li_test ( );
  c8_normal_01_test ( );
  c8_one_test ( );
  c8_print_test ( );
  c8_real_test ( );
  c8_sin_test ( );
  c8_sinh_test ( );
  c8_sqrt_test ( );
  c8_sub_test ( );
  c8_tan_test ( );
  c8_tanh_test ( );
  c8_to_cartesian_test ( );
  c8_to_polar_test ( );
  c8_uniform_01_test ( );
  c8_zero_test ( );

  c8mat_identity_test ( );
  c8mat_indicator_test ( );
  c8mat_norm_fro_test ( );
  c8mat_norm_l1_test ( );
  c8mat_norm_li_test ( );
  c8mat_uniform_01_test ( );

  c8vec_indicator_test ( );
  c8vec_nint_test ( );
  c8vec_norm_l1_test ( );
  c8vec_norm_l2_test ( );
  c8vec_norm_li_test ( );
  c8vec_print_test ( );
  c8vec_sort_a_l1_test ( );
  c8vec_sort_a_l2_test ( );
  c8vec_sort_a_li_test ( );
  c8vec_spiral_test ( );
  c8vec_uniform_01_test ( );
  c8vec_unity_test ( );

  cartesian_to_c8_test ( );

  polar_to_c8_test ( );

  r8_atan_test ( );
  r8_sign_test ( );
  r8_uniform_01_test ( );

  r8poly2_root_test ( );
  r8poly3_root_test ( );
  r8poly4_root_test ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8LIB_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function [params,X,means,wX,wM,dwM] = pre_sphere_symm(params, X, dim)
% Spheres the data with covarianve matrix eigenvalue decomposition.
%   [params,X,means,wX,wM,dwM] = pre_sphere_symm(params, X, dim)
%     X    Data to be sphered
%     dim  Requested result data dimension
%     wX   Sphered data
%     wM   Sphering matrix
%     dwM  De-sphering matrix

% Copyright (C) 2004, 2005 DSS MATLAB package team (dss@cis.hut.fi).
% Distributed by Laboratory of Computer and Information Science,
% Helsinki University of Technology. http://www.cis.hut.fi/projects/dss/.
% $Id$

if nargin<2
    params.name = 'Symmetric sphering';
    params.description = 'Spheres the data with covarianve matrix eigenvalue decomposition.';
    return;
end

if nargin<3; dim=0; end

[xdim,tdim] = size(X);
% removing the mean
means = mean(X,2);
X = X - repmat(means,1,tdim);

[wM, dwM] = dss_sphere(X, dim, 1);

wX = wM * X;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function y = kron(X,Y)
%KRON (overloaded)

if isequal(X,1)
    y = Y;
    return;
end

if isequal(Y,1);
    y = X;
end

if (isa(X,'sdpvar') & isa(Y,'sdpvar'))
    y = [];
    s.type = '()';
    for i = 1:size(X,1)
        this_row = [];
        for j = 1:size(X,2);
            s.subs = {i,j};
            this_row = [this_row subsref(X,s)*Y];
        end
        y = [y;this_row];
    end
    return
end

if isa(X,'sdpvar')
    lmi_variables = getvariables(X);
    nv = length(lmi_variables);
    y  = X;
    sparse_Y = sparse(Y);
    % This one used also for checking size
    temp = kron(getbasematrix(X,0),sparse_Y);
    if size(Y,2)==1
        % Special case
        % [kron([N1(:) N2(:)],vec(M))]=[vec(kron(N1,M)) vec(kron(N2,M))]
        y.basis = kron(X.basis,sparse_Y);
    else
        y.basis=temp(:);
        X_base = X.basis;
        for i = 1:nv
           % temp = kron(getbasematrix(X,lmi_variables(i)),sparse_Y);
            temp = kron(reshape(X_base(:,i+1),X.dim),sparse_Y);
            y.basis(:,i+1) = temp(:);
        end;
    end
end

if isa(Y,'sdpvar')
    lmi_variables = getvariables(Y);
    nv = length(lmi_variables);
    y  = Y;
    sparse_X = sparse(X);
    % This one used also for checking size
    temp = kron(sparse_X,getbasematrix(Y,0));
    if size(X,1)==1
        % In this special case
        %[kron(vec(M),[N1(:) N2(:)])]==[vec(kron(M,N1)) vec(kron(M,N2))]
        y.basis = kron(sparse_X(:),Y.basis);
    else
        y.basis = temp(:);
        y.basis = spalloc(length(temp(:)),nv+1,0);
        y.basis(:,1) = temp(:);
        Y_base = Y.basis;
        for i = 1:nv
            % temp = kron(sparse_X,getbasematrix(Y,lmi_variables(i)));
            temp = kron(sparse_X,reshape(Y_base(:,i+1),Y.dim));
            y.basis(:,i+1) = temp(:);
        end
    end
end
y.dim(1) = size(temp,1);
y.dim(2) = size(temp,2);
y = clean(y);
% Reset info about conic terms
if isa(y,'sdpvar')
    y.conicinfo = [0 0];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%% Hough Circle Transform
% An example using the Hough circle detector.
%
% This program demonstrates circle finding with the Hough transform. We show
% how to use the OpenCV function |cv.HoughCircles| to detect circles in an
% image.
%
% Sources:
%
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/houghcircles.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/python/houghcircles.py>
% * <https://docs.opencv.org/3.2.0/d4/d70/tutorial_hough_circle.html>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/ImgTrans/HoughCircle_Demo.cpp>
% * <https://docs.opencv.org/3.2.0/da/d53/tutorial_py_houghcircles.html>
%

%% Theory
%
% The Hough Circle Transform works in a roughly analogous way to the
% Hough Line Transform.
%
% In the line detection case, a line was defined by two parameters
% $(r, \theta)$. In the circle case, we need three parameters to define a
% circle:
%
% $$C : ( x_{center}, y_{center}, r )$$
%
% where $(x_{center}, y_{center})$ define the center position (green point)
% and $r$ is the radius, which allows us to completely define a circle, as it
% can be seen below:
%
% <<https://docs.opencv.org/3.2.0/Hough_Circle_Tutorial_Theory_0.jpg>>
%
% For sake of efficiency, OpenCV implements a detection method slightly
% trickier than the standard Hough Transform: _The Hough gradient method_,
% which is made up of two main stages:
%
% * the first stage involves edge detection and finding the possible circle
%   centers
% * the second stage finds the best radius for each candidate center
%
% For more details, please check the book _Learning OpenCV_ or your favorite
% Computer Vision bibliography.
%

%% Code
%
% This program:
%
% * Loads an image and blur it to reduce the noise
% * Applies the _Hough Circle Transform_ to the blurred image .
% * Display the detected circle in a window.
%

%%
% Input image
if ~mexopencv.isOctave() && mexopencv.require('images')
    fname = which('coins.png');
else
    fname = fullfile(mexopencv.root(), 'test', 'detect_blob.png');
end
assert(~isempty(fname) && exist(fname,'file')==2, 'Image not found');
img = cv.imread(fname, 'Color',true);
figure, imshow(img), title('source')

%%
% convert it to grayscale, and reduce noise to avoid false circles detection
gray = cv.cvtColor(img, 'RGB2GRAY');
if true
    gray = cv.medianBlur(gray, 'KSize',5);
else
    gray = cv.GaussianBlur(gray, 'KSize',[7,7], 'SigmaX',0.9, 'SigmaY',0.9);
end

%%
% Hough Circle Transform
%
% * change |MinDist| to detect circles with different distances to each other
% * change |MinRadius| and |MaxRadius| to detect circles of different sizes
tic
circles = cv.HoughCircles(gray, 'Method','Gradient', 'DP',2, ...
    'MinDist',size(gray,1)/8, 'Param1',200, 'Param2',100, ...
    'MinRadius',20, 'MaxRadius',80);
toc

%%
% draw detected circles (outlines and centers), and display the result
circles = cat(1, circles{:});
center = round(circles(:,1:2));
radius = round(circles(:,3));
out = cv.circle(img, center, radius, 'Color',[0 0 255], ...
    'Thickness',2, 'LineType','AA');
out = cv.circle(out, center, 3, 'Color',[0 255 0], ...
    'Thickness','Filled', 'LineType','AA');
figure, imshow(out), title('detected circles')

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function ArgStruct=parseArgs(args,ArgStruct,varargin)
% Helper function for parsing varargin.
%
%
% ArgStruct=parseArgs(varargin,ArgStruct[,FlagtypeParams[,Aliases]])
%
% * ArgStruct is the structure full of named arguments with default values.
% * Flagtype params is params that don't require a value. (the value will be set to 1 if it is present)
% * Aliases can be used to map one argument-name to several argstruct fields
%
%
% example usage:
% --------------
% function parseargtest(varargin)
%
% %define the acceptable named arguments and assign default values
% Args=struct('Holdaxis',0, ...
%        'SpacingVertical',0.05,'SpacingHorizontal',0.05, ...
%        'PaddingLeft',0,'PaddingRight',0,'PaddingTop',0,'PaddingBottom',0, ...
%        'MarginLeft',.1,'MarginRight',.1,'MarginTop',.1,'MarginBottom',.1, ...
%        'rows',[],'cols',[]);
%
% %The capital letters define abrreviations.
% %  Eg. parseargtest('spacingvertical',0) is equivalent to  parseargtest('sv',0)
%
% Args=parseArgs(varargin,Args, ... % fill the arg-struct with values entered by the user
%           {'Holdaxis'}, ... %this argument has no value (flag-type)
%           {'Spacing' {'sh','sv'}; 'Padding' {'pl','pr','pt','pb'}; 'Margin' {'ml','mr','mt','mb'}});
%
% disp(Args)
%
%
%
%
% % Aslak Grinsted 2003

Aliases={};
FlagTypeParams='';

if (length(varargin)>0)
    FlagTypeParams=strvcat(varargin{1});
    if length(varargin)>1
        Aliases=varargin{2};
    end
end

%---------------Get "numeric" arguments
NumArgCount=1;
while (NumArgCount<=size(args,2))&(~ischar(args{NumArgCount}))
    NumArgCount=NumArgCount+1;
end
NumArgCount=NumArgCount-1;
if (NumArgCount>0)
    ArgStruct.NumericArguments={args{1:NumArgCount}};
else
    ArgStruct.NumericArguments={};
end

%--------------Make an accepted fieldname matrix (case insensitive)
Fnames=fieldnames(ArgStruct);
for i=1:length(Fnames)
    name=lower(Fnames{i,1});
    Fnames{i,2}=name; %col2=lower
    AbbrevIdx=find(Fnames{i,1}~=name);
    Fnames{i,3}=[name(AbbrevIdx) ' ']; %col3=abreviation letters (those that are uppercase in the ArgStruct) e.g. SpacingHoriz->sh
    %the space prevents strvcat from removing empty lines
    Fnames{i,4}=isempty(strmatch(Fnames{i,2},FlagTypeParams)); %Does this parameter have a value? (e.g. not flagtype)
end
FnamesFull=strvcat(Fnames{:,2});
FnamesAbbr=strvcat(Fnames{:,3});

if length(Aliases)>0
    for i=1:length(Aliases)
        name=lower(Aliases{i,1});
        FieldIdx=strmatch(name,FnamesAbbr,'exact'); %try abbreviations (must be exact)
        if isempty(FieldIdx)
            FieldIdx=strmatch(name,FnamesFull); %&??????? exact or not?
        end
        Aliases{i,2}=FieldIdx;
        AbbrevIdx=find(Aliases{i,1}~=name);
        Aliases{i,3}=[name(AbbrevIdx) ' ']; %the space prevents strvcat from removing empty lines
        Aliases{i,1}=name; %dont need the name in uppercase anymore for aliases
    end
    %Append aliases to the end of FnamesFull and FnamesAbbr
    FnamesFull=strvcat(FnamesFull,strvcat(Aliases{:,1}));
    FnamesAbbr=strvcat(FnamesAbbr,strvcat(Aliases{:,3}));
end

%--------------get parameters--------------------
l=NumArgCount+1;
while (l<=length(args))
    a=args{l};
    if ischar(a)
        paramHasValue=1; % assume that the parameter has is of type 'param',value
        a=lower(a);
        FieldIdx=strmatch(a,FnamesAbbr,'exact'); %try abbreviations (must be exact)
        if isempty(FieldIdx)
            FieldIdx=strmatch(a,FnamesFull);
        end
        if (length(FieldIdx)>1) %shortest fieldname should win
            [mx,mxi]=max(sum(FnamesFull(FieldIdx,:)==' ',2));
            FieldIdx=FieldIdx(mxi);
        end
        if FieldIdx>length(Fnames) %then it's an alias type.
            FieldIdx=Aliases{FieldIdx-length(Fnames),2};
        end

        if isempty(FieldIdx)
            error(['Unknown named parameter: ' a])
        end
        for curField=FieldIdx' %if it is an alias it could be more than one.
            if (Fnames{curField,4})
                val=args{l+1};
            else
                val=1; %parameter is of flag type and is set (1=true)....
            end
            ArgStruct.(Fnames{curField,1})=val;
        end
        l=l+1+Fnames{FieldIdx(1),4}; %if a wildcard matches more than one
    else
        error(['Expected a named parameter: ' num2str(a)])
    end
end

author: konstantinberlin

repo: malware-windows-audit-log-detection


text:
function [x,pca_x,pca_coeff] = preprocess(x)
    [pca_coeff,pca_x,latent] = pca(x);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [count, h, parent, post, L] = symbfact2 (A, mode, Lmode)	    %#ok
%SYMBFACT2  symbolic factorization
%
%   Analyzes the Cholesky factorization of A, A'*A, or A*A'.
%
%   Example:
%   count = symbfact2 (A)               returns row counts of R=chol(A)
%   count = symbfact2 (A,'col')         returns row counts of R=chol(A'*A)
%   count = symbfact2 (A,'sym')         same as symbfact2(A)
%   count = symbfact2 (A,'lo')          same as symbfact2(A'), uses tril(A)
%   count = symbfact2 (A,'row')         returns row counts of R=chol(A*A')
%
%   The flop count for a subsequent LL' factorization is sum(count.^2)
%
%   [count, h, parent, post, R] = symbfact2 (...) returns:
%
%       h: height of the elimination tree
%       parent: the elimination tree itself
%       post: postordering of the elimination tree
%       R: a 0-1 matrix whose structure is that of chol(A) for the symmetric
%           case, chol(A'*A) for the 'col' case, or chol(A*A') for the
%           'row' case.
%
%   symbfact2(A) and symbfact2(A,'sym') uses the upper triangular part of A
%   (triu(A)) and assumes the lower triangular part is the transpose of
%   the upper triangular part.  symbfact2(A,'lo') uses tril(A) instead.
%
%   With one to four output arguments, symbfact2 takes time almost proportional
%   to nnz(A)+n where n is the dimension of R, and memory proportional to
%   nnz(A).  Computing the 5th argument takes more time and memory, both
%   O(nnz(L)).  Internally, the pattern of L is computed and R=L' is returned.
%
%   The following forms return L = R' instead of R.  They are faster and take
%   less memory than the forms above.  They return the same count, h, parent,
%   and post outputs.
%
%   [count, h, parent, post, L] = symbfact2 (A,'col','L')
%   [count, h, parent, post, L] = symbfact2 (A,'sym','L')
%   [count, h, parent, post, L] = symbfact2 (A,'lo', 'L')
%   [count, h, parent, post, L] = symbfact2 (A,'row','L')
%
%   See also CHOL, ETREE, TREELAYOUT, SYMBFACT

%   Copyright 2006-2007, Timothy A. Davis
%   http://www.cise.ufl.edu/research/sparse

error ('symbfact2 mexFunction not found!') ;


text:
function [node,elem]=meshresample(v,f,keepratio)
%
% [node,elem]=meshresample(v,f,keepratio)
%
% resample mesh using CGAL mesh simplification utility
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2007/11/12
%
% input:
%    v: list of nodes
%    f: list of surface elements (each row for each triangle)
%    keepratio: decimation rate, a number less than 1, as the percentage
%               of the elements after the sampling
%
% output:
%    node: the node coordinates of the sampled surface mesh
%    elem: the element list of the sampled surface mesh
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

[node,elem]=domeshsimplify(v,f,keepratio);

if(length(node)==0)
    warning(['Your input mesh contains topological defects, and the ',...
           'mesh resampling utility aborted during processing. Now iso2mesh ',...
           'is trying to repair your mesh with meshcheckrepair. ',...
           'You can also call this manually before passing your mesh to meshresample.'] );
    [vnew,fnew]=meshcheckrepair(v,f);
    [node,elem]=domeshsimplify(vnew,fnew,keepratio);
end
[node,I,J]=unique(node,'rows');
elem=J(elem);
saveoff(node,elem,mwpath('post_remesh.off'));

end

% function to perform the actual resampling
function [node,elem]=domeshsimplify(v,f,keepratio)
  exesuff=getexeext;
  exesuff=fallbackexeext(exesuff,'cgalsimp2');

  saveoff(v,f,mwpath('pre_remesh.off'));
  deletemeshfile(mwpath('post_remesh.off'));
  system([' "' mcpath('cgalsimp2') exesuff '" "' mwpath('pre_remesh.off') '" ' num2str(keepratio) ' "' mwpath('post_remesh.off') '"']);
  [node,elem]=readoff(mwpath('post_remesh.off'));
end


text:
classdef IN1KMOP3 < EvoXBenchProblem
% <multi> <real> <large/none> <expensive/none>

%------------------------------- Reference --------------------------------
% Z. Lu, R. Cheng, Y. Jin, K. C. Tan, and K. Deb, Neural architecture
% search as multiobjective optimization benchmarks: Problem formulation and
% performance assessment, IEEE Transactions on Evolutionary Computation,
% 2023.
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            config.name = 'resnet';
            config.args.objs = 'err&params&flops';
            config.args.normalized_objectives = false;
            obj.Setting@EvoXBenchProblem(config);
        end
    end
end


text:
function plot3(catalogObject, varargin)
    %CATALOG.PLOT3 Plot hypocenters in 3-D
    %   catalogObject.plot3()
    %
    %   Optional name/value pairs:
    %     'nsigma' - controls how zoomed-in the axes are (default
    %     5)

    % Glenn Thompson 2014/06/01
    if all(isnan(catalogObject.lat))
        warning('No hypocenter data to plot');
        return
    end
    p = inputParser;
    p.addParamValue('nsigma', '5', @isstr);
    p.parse(varargin{:});
    nsigma = p.Results.nsigma;

    % change region
    region = get_region(catalogObject, nsigma);

    % Compute Marker Size
    symsize = get_symsize(catalogObject);

    % 3D plot
    figure
    set(gcf,'Color', [1 1 1]);
    scatter3(catalogObject.lon, catalogObject.lat, catalogObject.depth, symsize);
    set(gca, 'ZDir', 'reverse');
    grid on;
    set(gca, 'XLim', [region(1) region(2)]);
    set(gca, 'YLim', [region(3) region(4)]);
    xlabel('Longitude');
    ylabel('Latitude');
    zlabel('Depth (km)');

end

author: geoscience-community-codes

repo: GISMO


text:
function test_xformInplaneToVolume
%Validate transforming data from INPLANE view to VOLUME (gray) view
%
%  test_xformInplaneToVolume()
%
% Tests: initHiddenGray, loadMeanMap, loadCorAnal, ip2volParMap,
% ip2volCorAnal, ip2volTSeries
%
% INPUTS
%  No inputs
%
% RETURNS
%  No returns
%
% Example: test_getCurDataROI()
%
% See also MRVTEST
%
% Copyright Stanford team, mrVista, 2011

%% Initialize the key variables and data path
% Data directory (where the mrSession file is located)
dataDir = mrtInstallSampleData('functional','mrBOLD_01');

% This is the validation file
val = mrtGetValididationData('ipToVolumeData');

% These are the items we stored in the validation file
%
% val.codim      = size(co);
% val.comed      = nanmedian(co);
% val.cosample   = co(1000);
% val.mapdim     = size(map);
% val.mapmed     = nanmedian(map);
% val.mapsample  = map(1000);
% val.tSdim     = size(tSeries);
% val.tSsample  = tSeries(1000);
%
% save(vFile, '-struct', 'val')

%% Retain original directory, change to data directory
curDir = pwd;
cd(dataDir);

% There can be several data types - name the one you want to probe
dataType = 'Original';

% Which scan number from that data type?
scan = 1;

%% Get data structure:
ip  = initHiddenInplane(); % Foregoes interface - loads data silently
vol = initHiddenGray();

%% Set dataTYPE:
ip  = viewSet(ip, 'Current DataType', dataType); % Data type
vol = viewSet(vol, 'Current DataType', dataType); % Data type

%% Load coranal and mean map into INPLANE view
ip = loadCorAnal(ip);
ip = loadMeanMap(ip);

%% Transform data to VOLUME view
vol = ip2volParMap( ip, vol, scan, [], 'linear', 1);
vol = ip2volCorAnal(ip, vol, scan, -1);
vol = ip2volTSeries(ip, vol, scan);

co  = viewGet(vol, 'scan coherence', scan);
map = viewGet(vol, 'scan map', scan);
tSeries = loadtSeries(vol, scan, 1); %We can hardcode slice since gray view only has 1 slice

%% Go home
cd(curDir);

%% Validate..

% check the dimensions of coherence map and mean map
assertEqual(val.codim, size(co));
assertEqual(val.mapdim, size(map));
assertEqual(val.tSdim, size(tSeries));

% check the median value of the coherence map and mean map
assertElementsAlmostEqual(val.comed,nanmedian(co));
assertElementsAlmostEqual(val.mapmed,nanmedian(map));

% check the values of the mean map and coherence map for an arbitrary voxel
% (to make sure the sequence is correct)
assertElementsAlmostEqual(val.cosample, co(1000));
assertElementsAlmostEqual(val.mapsample, map(1000));

% check the values of the time series for the max and an arbitrary voxel
assertElementsAlmostEqual(val.tSsample, tSeries(1000));
assertElementsAlmostEqual(val.tSmax,nanmax(tSeries(:)));

%% Cleanup
% clean up vistadata repository because this test script wrote new data
% test_CleanUpSVN

mrvCleanWorkspace;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function model = spm_mvb(X,Y,X0,U,V,nG,sG)
% Bayesian optimisation of a multivariate linear model with a greedy search
% FORMAT model = spm_mvb(X,Y,X0,U,V,nG,sG)
%
% X      - contrast or target vector
% Y      - date feature matrix
% X0     - confounds
% U      - patterns
% V      - observation noise covariance
% nG     - number of Greedy iterations (nG = 1 => uniform hyperpriors)
%        - if not specified, the search will terminate when F falls
% sG     - size of successive subdivisions [default is 1/2)
%
% returns model:
%                F: log-evidence [F(0), F(1),...]
%                G: covariance partition indices
%                h: covariance hyperparameters
%                U: ordered patterns
%                M: MAP projector: qE = M*X
%               qE: conditional expectation of voxel weights
%               qC: conditional variance of voxel weights
%               Cp: empirical prior covariance (ordered  pattern space)
%               cp: empirical prior covariance (original pattern space)
%__________________________________________________________________________
%
% model: X = Y*P + X0*Q + R
%        P = U*E;
%   cov(E) = h1*diag(G(:,1)) + h2*diag(G(:,2)) + ...
%
% This routine uses a multivariate Bayesian (MVB) scheme to decode or
% recognise brain states from neuroimages. It resolves the ill-posed
% many-to-one mapping, from voxel values or data features to a target
% variable, using a parametric empirical or hierarchical Bayesian model.
% This model is inverted using standard variational techniques, in this
% case Variational Laplace, to furnish the model evidence and the
% conditional density of the model's parameters. This allows one to compare
% different models or hypotheses about the mapping from functional or
% structural anatomy to perceptual and behavioural consequences (or their
% deficits). The aim of MVB is not to predict (because the outcomes are
% known) but to enable inference on different models of structure-function
% mappings; such as distributed and sparse representations. This allows one
% to optimise the model itself and produce predictions that outperform
% standard pattern classification approaches, like support vector machines.
% Technically, the model inversion and inference uses the same empirical
% Bayesian procedures developed for ill-posed inverse problems (e.g.,
% source reconstruction in EEG).
%
% CAUTION: MVB should not be used to establish a significant mapping
% between brain states and some classification or contrast vector. Its use
% is limited to comparison of different models under the assumption
% (hyperprior) that this mapping exists. To ensure the mapping exists, use
% CVA or related approaches.
%
% See spm_mvb_ui and:
%
% Bayesian decoding of brain images.
% Friston K, Chu C, Mourão-Miranda J, Hulme O, Rees G, Penny W, Ashburner J.
% Neuroimage. 2008 Jan 1;39(1):181-205
%
% Multiple sparse priors for the M/EEG inverse problem.
% Friston K, Harrison L, Daunizeau J, Kiebel S, Phillips C, Trujillo-Barreto
% N, Henson R, Flandin G, Mattout J.
% Neuroimage. 2008 Feb 1;39(3):1104-20.
%
% Characterizing dynamic brain responses with fMRI: a multivariate approach.
% Friston KJ, Frith CD, Frackowiak RS, Turner R.
% Neuroimage. 1995 Jun;2(2):166-72.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_mvb.m 7679 2019-10-24 15:54:07Z spm $

% defaults (use splits +/- one standard deviation by default)
%--------------------------------------------------------------------------
try, V;          catch, V  = [];             end
try, nG; aG = 0; catch, nG = 8; aG = 1;      end
try, sG;         catch, sG = 1/2;            end

% get orders
%--------------------------------------------------------------------------
ns     = size(Y,1);                 % number of samples
nv     = size(Y,2);                 % number of voxels
np     = size(U,2);                 % number of patterns or parameters

% confounds
%--------------------------------------------------------------------------
if isempty(X0), X0 = zeros(ns,1);  end
if isempty(U),  U  = zeros(nv,0);  end
if isempty(V),  V  = speye(ns,ns); end

% number of error components
%--------------------------------------------------------------------------
if iscell(V)
    nh = length(V);
else
    nh = 1;
end

% null model
%--------------------------------------------------------------------------
model  = spm_mvb_G(X,zeros(ns,0),X0,[],V);
if ~np, return, end

% initialise G and greedy search
%==========================================================================
F  = model.F;
L  = Y*U;
G  = ones(np,1);
for  i = 1:nG

    % invert
    %----------------------------------------------------------------------
    M  = spm_mvb_G(X,L,X0,G,V);

    % record conditional estimates (and terminate automatically if aG)
    %----------------------------------------------------------------------
    if i == 1
        save_model = 1;
    else
        save_model = M.F > max(F);
    end
    if  save_model
        model      = M;
    elseif aG
        break
    end

    % record free-energy
    %----------------------------------------------------------------------
    F(i + 1)     = M.F;
    lnh          = log(M.h');

    disp('log evidence & hyperparameters:')
    fprintf('% 8.2f',F-F(1)),fprintf('\n')
    fprintf('% 8.2f',full(lnh)),fprintf('\n\n')

    % eliminate redundant components
    %----------------------------------------------------------------------
    lnh          = lnh((nh + 1):end) - lnh(end);
    j            = find(lnh < -16);
    G(:,j)       = [];

    % create new spatial support
    %----------------------------------------------------------------------
    g            = find(G(:,end));
    ng           = ceil(length(g)*sG);
    [q,j]        = sort(-sum(M.qE(g,:).^2,2));
    q            = g(j(1:ng));
    G(q,end + 1) = 1;

    % break if cluster is one
    %----------------------------------------------------------------------
    if ng < 1/(1 - sG), break, end

end

% project pattern weights to feature (voxel) weights
%==========================================================================

% remove some patterns if there are too many
%--------------------------------------------------------------------------
clear M X Y
qE       = sum(model.qE.^2,2);
[i,j]    = sort(-qE);
try
    i    = j(1:2^11);
catch
    i    = j;
end
L        = L(:,i);
U        = U(:,i);
cp       = model.Cp;
Cp       = cp(i,i);
MAP      = U*model.MAP(i,:);

% try to save conditional expectations (if there is enough memory)
%--------------------------------------------------------------------------
try
    qE   = U*model.qE(i,:);
catch
    qE   = [];
end

% remove confounds from L = Y*U
%--------------------------------------------------------------------------
L        = L - X0*pinv(full(X0))*L;

% conditional covariance in voxel space: qC  = U*( Cp - Cp*L'*iC*L*Cp )*U';
%--------------------------------------------------------------------------
UCp      = U*Cp;
qC       = sum(UCp.*U,2) - sum((UCp*L').*MAP,2);

model.F  = F;
model.U  = U;
model.M  = MAP;                             % MAP projector
model.qE = qE;                              % conditional expectation
model.Cp = Cp;                              % prior covariance (ordered)
model.cp = cp;                              % prior covariance (original)
model.qC = max(qC,exp(-16));                % conditional variance

author: spm


text:
function [sFile, ChannelMat] = in_fopen_smr(DataFile)
% IN_FOPEN_SMR: Open a Cambridge Electronic Design Spike2 file (.smr/.son).

% @=============================================================================
% 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:  Malcolm Lidierth, 2006-2007, King's College London
%           Adapted by Francois Tadel for Brainstorm, 2017-2021

%% ===== READ HEADER =====
% Get file type
[fPath, fBase, fExt]=fileparts(DataFile);
if strcmpi(fExt,'.smr') || strcmpi(fExt,'.smrx')
    byteorder = 'l';  % Spike2 for Windows source file: little-endian
elseif strcmpi(fExt,'.son')
    byteorder = 'b';  % Spike2 for Mac file: Big-endian
else
    error('Not a Spike2 file.');
end
% Open file
fid = fopen(DataFile, 'r', byteorder);
if (fid == -1)
    error('Could not open file.');
end

% Get file header
hdr = SONFileHeader(fid);
% Get list of channels
iChan = 0;
iMarkerChan = [];
for i = 1:hdr.channels
    bst_progress('text', sprintf('Reading channel info... [%d%%]', round(i/hdr.channels*100)));
    % Read channel information
    c = SONChannelInfo(fid, i);
    % Check type of the channel
    switch c.kind
        % ADC channels: read as signals
        case {1,9}
            iChan = iChan + 1;
            hdr.chaninfo(iChan).number    = i;
            hdr.chaninfo(iChan).kind      = c.kind;
            hdr.chaninfo(iChan).title     = c.title;
            hdr.chaninfo(iChan).comment   = c.comment;
            hdr.chaninfo(iChan).phyChan   = c.phyChan;
            hdr.chaninfo(iChan).idealRate = c.idealRate;
            % Read blocks
            hdr.chaninfo(iChan).blocks = SONGetBlockHeaders(fid, i);
        % Markers and events
        case {2,3,4,5,6,7,8}
        	iMarkerChan(end+1) = i;
    end
end
hdr.num_channels = length(hdr.chaninfo);
% Get maximum sampling rate
[sfreq, chMax] = max([hdr.chaninfo.idealRate]);

%% ===== FILL STRUCTURE =====
% Initialize returned file structure
sFile = db_template('sfile');
% Add information read from header
sFile.byteorder    = byteorder;
sFile.filename     = DataFile;
sFile.format       = 'EEG-SMR';
sFile.prop.sfreq   = sfreq;
sFile.prop.times   = [0, sum(hdr.chaninfo(chMax).blocks(5,:)) - 2] ./ sFile.prop.sfreq;
sFile.prop.nAvg    = 1;
sFile.channelflag  = ones(hdr.num_channels,1);
sFile.device       = 'CED Spike2';
sFile.header       = hdr;
% Comment: short filename
[fPath, fBase, fExt] = bst_fileparts(DataFile);
sFile.comment = fBase;

%% ===== CREATE EMPTY CHANNEL FILE =====
ChannelMat = db_template('channelmat');
ChannelMat.Comment = [sFile.device ' channels'];
ChannelMat.Channel = repmat(db_template('channeldesc'), [1, hdr.num_channels]);
% For each channel
for iChan = 1:hdr.num_channels
    ChannelMat.Channel(iChan).Type = 'EEG';
    ChannelMat.Channel(iChan).Name = hdr.chaninfo(iChan).title;
end

%% ===== READ MARKER INFORMATION =====
iChanMissing = [];
for iChan = 1:length(iMarkerChan)
    % Read channel
    try
        [d,header] = SONGetChannel(fid, iMarkerChan(iChan));
    catch
        iChanMissing = [iChanMissing, iMarkerChan(iChan)];
        continue;
    end
    if isempty(d) || isempty(header)
        continue;
    end
    % Get timing
    switch (header.kind)
        case {2,3,4}
            timeEvt = d(:)';
        case {5,6,7,8}
            timeEvt = d.timings(:)';
    end
    % Create event structure
    iEvt = length(sFile.events) + 1;
    if ~isempty(header.title)
        sFile.events(iEvt).label = header.title;
    else
        sFile.events(iEvt).label = sprintf('unknown_%02d', iEvt);
    end
    sFile.events(iEvt).times    = round(double(timeEvt).* sFile.prop.sfreq) ./ sFile.prop.sfreq;
    sFile.events(iEvt).epochs   = ones(size(sFile.events(iEvt).times));
    sFile.events(iEvt).select   = 1;
    sFile.events(iEvt).channels = [];
    sFile.events(iEvt).notes    = [];
end
% Display missing channels
if ~isempty(iChanMissing)
    disp(['SON> Missing channels: ' sprintf('%d ', iChanMissing)]);
end

% Close file
fclose(fid);


text:
function [V_hat, lam_hat] = cpf_predictor(V, lam, z, step, pv, pq)
%CPF_PREDICTOR  Performs the predictor step for the continuation power flow
%   [V_HAT, LAM_HAT] = CPF_PREDICTOR(V, LAM, Z, STEP, PV, PQ)
%
%   Computes a prediction (approximation) to the next solution of the
%   continuation power flow using a normalized tangent predictor.
%
%   Inputs:
%       V : complex bus voltage vector at current solution
%       LAM : scalar lambda value at current solution
%       Z : normalized tangent prediction vector from previous step
%       STEP : continuation step length
%       PV : vector of indices of PV buses
%       PQ : vector of indices of PQ buses
%
%   Outputs:
%       V_HAT : predicted complex bus voltage vector
%       LAM_HAT : predicted lambda continuation parameter

%   MATPOWER
%   Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC)
%   by Shrirang Abhyankar, Argonne National Laboratory
%   and Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% sizes
nb = length(V);

Va = angle(V);
Vm = abs(V);
Va_hat = Va;
Vm_hat = Vm;

%% prediction for next step
Va_hat([pv; pq]) = Va([pv; pq]) + step * z([pv; pq]);
Vm_hat([pq])     = Vm([pq])     + step * z([nb+pq]);
lam_hat = lam + step * z(2*nb+1);
V_hat = Vm_hat .* exp(1j * Va_hat);

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function a = lesp ( m, n )

%*****************************************************************************80
%
%% LESP returns the LESP matrix.
%
%  Formula:
%
%    if ( I - J == 1 )
%      A(I,J) = 1 / I
%    else if ( I - J == 0 )
%      A(I,J) = - ( 2*I+3 )
%    else if ( I - J == 1 )
%      A(I,J) = J
%    else
%      A(I,J) = 0.0
%
%  Example:
%
%    M = 5, N = 5
%
%     -5    2    .    .     .
%     1/2  -7    3    .     .
%      .   1/3  -9    4     .
%      .    .   1/4 -11     5
%      .    .    .   1/5  -13
%
%
%  Properties:
%
%    The matrix is tridiagonal.
%
%    Because A is tridiagonal, it has property A (bipartite).
%
%    A is generally not symmetric: A' /= A.
%
%    The eigenvalues are real, and smoothly distributed in [-2*N-3.5, -4.5].
%
%    The eigenvalues are sensitive.
%
%    The matrix is similar to the symmetric tridiagonal matrix with
%    the same diagonal entries and with off-diagonal entries 1,
%    via a similarity transformation using the diagonal matrix
%    D = diagonal ( 1!, 2!, ..., N! ).
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Wim Lenferink, MN Spijker,
%    On the use of stability regions in the numerical analysis of initial
%    value problems,
%    Mathematics of Computation,
%    Volume 57, 1991, pages 221-237.
%
%    Lloyd Trefethen,
%    Pseudospectra of matrices,
%    in Numerical Analysis 1991,
%    Proceedings of the 14th Dundee Conference,
%    D F Griffiths and G A Watson, editors,
%    Pitman Research Notes in Mathematics, volume 260,
%    Longman Scientific and Technical, Essex, UK, 1992, pages 234-266.
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of A.
%
%    Output, real A(M,N), the matrix.
%
  a = zeros ( m, n );

  for i = 1 : m
    for j = 1 : n

      if ( i - j == 1 )
        a(i,j) = 1.0 / i;
      elseif ( i - j == 0 )
        a(i,j) = - ( 2 * i + 3 );
      elseif ( i - j == -1 )
        a(i,j) = j;
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end


text:
function [R] = averageRotationTransform(R1s,R2s)
%AVERAGEROTATION Summary of this function goes here
%   takes from R1s to R2s
[R1s,R2s] = removeNanRots(R1s,R2s);
R = eye(3);
N = length(R1s);
maxDiff = pi/3;

for i=1:N
    diffs{i} = R1s{i}'*R2s{i};
end
err = Inf;
iter = 0;

while(err > 1e-4 && iter < 50)
    iter = iter+1;
    delta = zeros(3,3);ct = 0;
    for i=1:N
        erri = norm(logm(R'*diffs{i}),'fro')/sqrt(2);
        if(iter < 5 || erri < maxDiff)
            delta = delta + logm(R'*diffs{i});
            ct = ct+1;
        end
    end
    delta = (delta/ct);
    err = (norm(delta,'fro')/sqrt(2));
    %disp(err);
    delta = expm(delta);
    R = R*delta;
    iter = iter + 1;
end
R = real(R);

end

function [nR1s,nR2s]=removeNanRots(R1s,R2s)
    n1 = cellfun(@(x)(any(isnan(x(:)))),R1s(:));
    n2 = cellfun(@(x)(any(isnan(x(:)))),R2s(:));
    n12 = n1 | n2;
    nR1s = R1s(~n12);
    nR2s = R2s(~n12);
end

repo: CategoryShapes


text:
function display(pot)

disp('utility potential object');
disp(struct(pot));

author: bayesnet

repo: bnt


text:
function value = rayleigh ( n, a, x )

%*****************************************************************************80
%
%% RAYLEIGH returns the Rayleigh quotient of the matrix A and the vector X.
%
%  Formula:
%
%    RAYLEIGH = X' * A * X / ( X' * X )
%
%  Properties:
%
%    If X is an eigenvector of A, then RAYLEIGH will equal the
%    corresponding eigenvalue.
%
%    The set of all Rayleigh quotients for a matrix is known
%    as its "field of values".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real A(N,N), the matrix.
%
%    Input, real X(N), the vector used in the Rayleigh quotient.
%
%    Output, real VALUE, the Rayleigh quotient of A and X.
%
  value = ( x(1:n) * a(1:n,1:n) * x(1:n))') /  ( x(1:n) * x(1:n)' );

  return
end


text:
%%*******************************************************************
%% HKMrhsfun: compute the right-hand side vector of the
%%            Schur complement equation for the  HKM direction.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*******************************************************************

function [rhs,EinvRc,hRd] = HKMrhsfun(blk,At,par,X,Z,rp,Rd,sigmu,hRd,dX,dZ)

m = length(rp);
if (nargin > 8)
    corrector = 1;
else
    corrector = 0;
    hRd = zeros(m,1);
end
hEinvRc = zeros(m,1);
EinvRc = cell(size(blk,1),1);
rhsfree = [];
%%
for p = 1:size(blk,1)
    pblk = blk(p,:);
    n = sum(pblk{2});
    if strcmp(pblk{1},'l')
        if iscell(sigmu)
            EinvRc{p} = sigmu{p}./Z{p} -X{p};
        else
            EinvRc{p} = sigmu./Z{p} -X{p};
        end
        Rq = sparse(n,1);
        if (corrector) && (norm(par.parbarrier{p})==0)
            Rq = dX{p}.*dZ{p}./Z{p};
        else
            tmp  = par.dd{p}.*Rd{p};
            tmp2 = mexMatvec(At{p},tmp,1);
            hRd = hRd + tmp2;
        end
        EinvRc{p} = EinvRc{p} - Rq;
        tmp2 = mexMatvec(At{p,1},EinvRc{p},1);
        hEinvRc = hEinvRc + tmp2;
    elseif strcmp(pblk{1},'q')
        if iscell(sigmu)
            EinvRc{p} = qops(pblk,sigmu{p},par.Zinv{p},3) -X{p};
        else
            EinvRc{p} = sigmu*par.Zinv{p} -X{p};
        end
        Rq = sparse(n,1);
        if (corrector) && (norm(par.parbarrier{p})==0)
            ff{p} = qops(pblk,1./par.gamz{p},Z{p},3); %#ok
            hdx = qops(pblk,par.gamz{p},ff{p},5,dX{p});
            hdz = qops(pblk,par.gamz{p},ff{p},6,dZ{p});
            hdxdz = Arrow(pblk,hdx,hdz);
            Rq = qops(pblk,par.gamz{p},ff{p},6,hdxdz);
        else
            tmp = par.dd{p}.*Rd{p} ...
                + qops(pblk,qops(pblk,Rd{p},par.Zinv{p},1),X{p},3) ...
                + qops(pblk,qops(pblk,Rd{p},X{p},1),par.Zinv{p},3);
            tmp2 = mexMatvec(At{p,1},tmp,1);
            hRd = hRd + tmp2;
        end
        EinvRc{p} = EinvRc{p} - Rq;
        tmp2 = mexMatvec(At{p,1},EinvRc{p},1);
        hEinvRc = hEinvRc + tmp2;
    elseif strcmp(pblk{1},'s')
        if iscell(sigmu)
            %%ss = [0,cumsum(pblk{2})];
            %%sigmuvec = zeros(n,1);
            %%for k = 1:length(pblk{2});
            %%   sigmuvec(ss(k)+1:ss(k+1)) = sigmu{p}(k)*ones(pblk{2}(k),1);
            %%end
            sigmuvec = mexexpand(pblk{2},sigmu{p});
            EinvRc{p} = par.Zinv{p}*spdiags(sigmuvec,0,n,n) -X{p};
        else
            EinvRc{p} = sigmu*par.Zinv{p} -X{p};
        end
        Rq = sparse(n,n);
        if (corrector) && (norm(par.parbarrier{p})==0)
            Rq = Prod3(pblk,dX{p},dZ{p},par.Zinv{p},0);
            Rq = 0.5*(Rq+Rq');
        else
            tmp = Prod3(pblk,X{p},Rd{p},par.Zinv{p},0,par.nzlistAy{p});
            tmp = 0.5*(tmp+tmp');
            tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),{tmp});
            hRd = hRd + tmp2;
        end
        EinvRc{p} = EinvRc{p} - Rq;
        tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),EinvRc(p));
        hEinvRc = hEinvRc + tmp2;
    elseif strcmp(pblk{1},'u')
        rhsfree = [rhsfree; Rd{p}]; %#ok
    end
end
%%
rhs = rp + hRd - hEinvRc;
rhs = full([rhs; rhsfree]);
%%*******************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [xRot,rotMat]=GCRS2MOD(xVec,TT1,TT2)
%%GCRS2MOD Rotate a vector from the geocentric celestial reference system
%          (GCRS) to the mean of date (MOD) coordinate system, which is
%          the coordinate system using the mean equinox and
%          ecliptic of date, IAU 2006/2000A model. The transformation is
%          performed by adding the precession and frame bias.
%
%INPUTS: xVec The 3XN matrix of N 3X1 Cartesian vectors that are to be
%             rotated from the GCRS coordinate system into the MOD
%             coordinate system.
% TT1, TT2 Jul1,Jul2 Two parts of a Julian date given in TT. The units
%             of the date are days. The full date is the sum of both
%             terms. The date is broken into two parts to provide more
%             bits of precision. It does not matter how the date is
%             split.
%
%OUTPUTS: xRot The 3XN matrix of the N 3X1 input vector rotated into the
%              MOD coordinate system.
%       rotMat The 3X3 rotation matrix such that
%              xRot(:,i)=rotMat*xVec(:,i).
%
%This uses functions in the International Astronomical Union's (IAU)
%Standard's of Fundamental Astronomy (SOFA) library to obtain the product
%of the precession rotation matrix and the frame rotation bias matrix. One
%goes from GCRS to mean of date by applying a frame bias and then
%precession. Thus this function removes those rotations. The rotations are
%discussed in the documentation for the SOFA library as well as in [1]
%among other sources.
%
%The algorithm can be compiled for use in Matlab  using the
%CompileCLibraries function.
%
%The algorithm is run in Matlab using the command format
%[xRot,rotMat]=GCRS2MOD(xVec,TT1,TT2);
%
%Different celestial coordinate systems are compared in [2].
%
%REFERENCES:
%[1] G. Petit and B. Luzum, IERS Conventions (2010), International Earth
%    Rotation and Reference Systems Service Std. 36, 2010.
%[2] D. F. Crouse, "An Overview of Major Terrestrial, Celestial, and
%    Temporal Coordinate Systems for Target Tracking," Formal Report,
%    Naval Research Laboratory, no. NRL/FR/5344--16-10,279, 10 Aug. 2016,
%    173 pages.
%
%March 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

error('This function is only implemented as a mexed C or C++ function. Please run CompileCLibraries.m to compile the function for use.')

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 c = mvnvbcost(x,Covx,mu,Covmu,pCholCov,pLogDetCov)
% c = mvnvbcost(x,Covx,mu,Covmu,pCholCov,pLogDetCov)
%
% Calculates c = -KL(q||p) = <log p(X)> - <log q(X)>
% where the expectation <.> is taken over q(X). This is used for
% calculating the lower bound of the marginal loglikelihood in VB
% models.
%
% p(X) = N(M,T)
% q(X) = N(x,Covx)
%
% Rest of the parameters are defined as:
% q(M) = N(mu,Covmu)
% <inv T> = inv(pCholCov * pCholCov')
% <log det T> = pLogDetCov
%
% If Covx==[], then the term <log q(X)> is not calculated.
% This is useful when X is, e.g., observations.

c = 0;
d = length(x);

% Cost from q-posterior
if ~isempty(Covx)
  % TODO:
  % if ~isempty(pCholCov), entropy = mvnentropy(pCholCov,true); else.... end
%  entropy = mvnentropy(pCholCov,true);
  entropy = mvnentropy(Covx);
  c = entropy;
else
  Covx = 0;
end

% Use Cholesky of the prior covariance
% $$$ opts.UT = true;
% $$$ opts.TRANSA = true;
z = solve_tril(pCholCov, x-mu);

% Cost from prior

% Below: (x-mu)' * Cov^(-1) * (x-mu) + trace(Cov^(-1)*(Covx+Covmu))
err2 = z'*z + trace(solve_triu(pCholCov',solve_tril(pCholCov,Covx+Covmu)));
%err2 = (x-mu)'*pInvCov*(x-mu) + (Covx(:)+Covmu(:))'*pInvCov(:);
%((x-mu).^2 + diag(Covx) + diag(Covmu));
c = c - 0.5*pLogDetCov - 0.5*err2 - 0.5*d*log(2*pi);


text:
function foot_clearance(nlp, bounds, frame)
    % constraints for swing foot clearance

    domain = nlp.Plant;
    x = domain.States.x;

    pos = getCartesianPosition(domain, frame);
    constraint_func = SymFunction(['foot_clearance_',domain.Name], pos(3), {x});

    % Foot Clearance Middle
    numNode = nlp.NumNode;
    addNodeConstraint(nlp, constraint_func, {'x'}, 'first', 0, 0, 'Nonlinear');
    addNodeConstraint(nlp, constraint_func, {'x'}, floor(numNode/2)+1, ...
        bounds.constrBounds.foot_clearance.lb, ...
        bounds.constrBounds.foot_clearance.ub,'Nonlinear');
    addNodeConstraint(nlp, constraint_func, {'x'}, 'last', 0, 0, 'Nonlinear');
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [kx, ky] = gradientKernels(sigma)
% Create kernels for computing gradient within 2D images.
%
%   KX = gradientKernels(SIGMA)
%   [KX, KY] = gradientKernels(SIGMA)
%   Generates the gradient kernel(s) used for computing gradients within 2D
%   images.
%   This is a low-level function that is called by the functions that
%   compute gradients.
%
%
%   Example
%   [kx, ky] = gradientKernels(2);
%
%   See also
%     imGradient, imGradientFilter, orientedGaussianKernel
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
% Created: 2022-06-01,    using Matlab 9.9.0.1570001 (R2020b) Update 4
% Copyright 2022 INRAE.

if nargin == 0 || sigma == 0
    % Default case: normalised sobel matrix
    kx = fspecial('sobel')' / 8;

else
    % compute size according to sigma
    Nx = ceil((3*sigma));
    lx = -Nx:Nx;
    ky = exp(-((lx / sigma) .^2) * .5);
    kx = -(lx / sigma) .* ky;
    kx = ky' * kx;
    kx = kx / sum(kx(kx > 0));
end

% optional kernel for gradient in Y-direction
if nargout > 1
    ky = kx';
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function vot_wrapper(tracker_name, runfile_name, do_cleanup)

if nargin < 3
    do_cleanup = true;
end

% *************************************************************
% VOT: Always call exit command at the end to terminate Matlab!
% *************************************************************
if do_cleanup
    cleanup = onCleanup(@() exit() );
else
    [pathstr, ~, ~] = fileparts(mfilename('fullpath'));
    cd_ind = strfind(pathstr, filesep());
    pathstr = pathstr(1:cd_ind(end)-1);
    cleanup = onCleanup(@() cd(pathstr));
end

try

% *************************************************************
% VOT: Set random seed to a different value every time.
% *************************************************************
RandStream.setGlobalStream(RandStream('mt19937ar', 'Seed', sum(clock)));

% **********************************
% VOT: Get initialization data
% **********************************
[images, region] = vot_initialize();

results = cell(length(images), 1);

bb_scale = 1;

% If the provided region is a polygon ...
if numel(region) > 4
    % Init with an axis aligned bounding box with correct area and center
    % coordinate
    cx = mean(region(1:2:end));
    cy = mean(region(2:2:end));
    x1 = min(region(1:2:end));
    x2 = max(region(1:2:end));
    y1 = min(region(2:2:end));
    y2 = max(region(2:2:end));
    A1 = norm(region(1:2) - region(3:4)) * norm(region(3:4) - region(5:6));
    A2 = (x2 - x1) * (y2 - y1);
    s = sqrt(A1/A2);
    w = s * (x2 - x1) + 1;
    h = s * (y2 - y1) + 1;
else
    cx = region(1) + (region(3) - 1)/2;
    cy = region(2) + (region(4) - 1)/2;
    w = region(3);
    h = region(4);
end

init_c = [cx cy];
init_sz = bb_scale * [w h];

im_size = size(imread(images{1}));
im_size = im_size([2 1]);

init_pos = min(max(round(init_c - (init_sz - 1)/2), [1 1]), im_size);
init_sz = min(max(round(init_sz), [1 1]), im_size - init_pos + 1);

seq.s_frames = images;
seq.init_rect = [init_pos, init_sz];

[file_path, file_name_start, file_ext] = fileparts(seq.s_frames{1});
[~, file_name_end, ~] = fileparts(seq.s_frames{end});
seq.path = file_path;
seq.name = 'vot_seq';
seq.ext = file_ext(2:end);
seq.len = length(seq.s_frames);
seq.nz = length(file_name_start);
seq.startFrame = str2num(file_name_start);
seq.endFrame = str2num(file_name_end);

% setup_tracker_paths(tracker_name);

otb_res = eval([runfile_name '(seq, [], []);']);
%convert the results to rectangle format
% otb_res = convert_to_rect(otb_res);

num_frames = numel(images);

for frame = 1:num_frames
    bb = otb_res.res(frame,:);
    sz = bb(3:4);
    c = bb(1:2) + (sz - 1)/2;
    new_sz = sz / bb_scale;
    new_tl = c - (new_sz - 1)/2;
    results{frame} = round([new_tl, new_sz]);
end

% **********************************
% VOT: Output the results
% **********************************
vot_quit(results);

catch err
    [wrapper_pathstr, ~, ~] = fileparts(mfilename('fullpath'));
    cd_ind = strfind(wrapper_pathstr, filesep());
    VOT_path = wrapper_pathstr(1:cd_ind(end));

    error_report_path = [VOT_path 'error_reports\'];
    if ~exist(error_report_path, 'dir')
        mkdir(error_report_path);
    end

    report_file_name = [error_report_path tracker_name '_' runfile_name datestr(now,'_yymmdd_HHMM') '.mat'];

    save(report_file_name, 'err')

    rethrow(err);
end

repo: STRCF