text:
% This script stores in a file the standard configuration of RP for
% sampling from the 4 segmentations discussed in the paper. The parameters
% were trained with the VOC 2007 dataset. Further details can be found in
% the paper:
%
%"S. Manen, M. Guillaumin, and L. Van Gool. Prime Object Proposals with
%Randomized Prim's Algorithm. In ICCV, 2013."

function generateRPConfig_4segs(rppath)

	%% Parameter specification:
	params.approxFinalNBoxes = 10000;                                    %Approximate number of proposals
	params.rSeedForRun = -1;                                            %Random seed to be used (-1 to generate it with a hashing function)
	params.q = 10;                                                      %Parameter to eliminate near duplicates (raise it to eliminate more duplicates)

	% LAB segmentation
	params.segmentations{1}.colorspace = 'LAB';                         %Colorspace: 'RGB', 'LAB', 'opponent', 'rg', 'HSV'
	% --> Segmentation parameters:
	params.segmentations{1}.superpixels.sigma = 0.8;
	params.segmentations{1}.superpixels.c = 100;
	params.segmentations{1}.superpixels.min_size = 100;
	% --> Parameters trained from VOC07:
	%   --> Feature weights:
	params.segmentations{1}.simWeights.wBias = 3.0017;
	params.segmentations{1}.simWeights.wCommonBorder = -1.0029;
	params.segmentations{1}.simWeights.wLABColorHist = -2.6864;
	params.segmentations{1}.simWeights.wSizePer = -2.3655;
	%   --> Size term alpha, as explained in the paper sec. 4.2.
	%       It is quantized to contain exactly 2^16 elements for speed
	%       purposes.
	params.segmentations{1}.alpha = dlmread(fullfile(rppath,'config','alpha/alpha_voc07.dat'));
	params.segmentations{1}.verbose = false;                            %Set to true to display more information during execution

	% Opponent
	params.segmentations{2} = params.segmentations{1};
	params.segmentations{2}.colorspace = 'opponent';                    %Colorspace: 'RGB', 'LAB', 'opponent', 'rg', 'HSV'

	% rg
	params.segmentations{3} = params.segmentations{1};
	params.segmentations{3}.colorspace = 'rg';                          %Colorspace: 'RGB', 'LAB', 'opponent', 'rg', 'HSV'

	% Hue
	params.segmentations{4} = params.segmentations{1};
	params.segmentations{4}.colorspace = 'HSV';                         %Colorspace: 'RGB', 'LAB', 'opponent', 'rg', 'HSV'

	%% Save parameters:
	save(fullfile(rppath, 'config', 'rp_4segs.mat'),'params');
end

author: Cloud-CV

repo: object-proposals


text:
function kern = rbfardKernParamInit(kern)

% RBFARDKERNPARAMINIT RBFARD kernel parameter initialisation.
% The automatic relevance determination version of the radial basis
% function kernel (RBFARD) is a very smooth non-linear kernel and is a
% popular choice for generic use.
%
% k(x_i, x_j) = sigma2 * exp(-gamma/2 *(x_i - x_j)'*A*(x_i - x_j))
%
% The parameters are sigma2, the process variance (kern.variance), the
% diagonal matrix of input scales (kern.inputScales, constrained to be
% between zero and one) and gamma, the inverse width
% (kern.inverseWidth). The inverse width controls how wide the basis
% functions are, the larger gamma, the smaller the basis functions
% are.
%
% SEEALSO : rbfKernParamInit
%
% FORMAT
% DESC initialises the automatic relevance determination radial basis function
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006

% KERN

% This parameter is restricted positive.
kern.inverseWidth = 2/kern.inputDimension;
kern.variance = 1;
% These parameters are restricted to lie between 0 and 1.
kern.inputScales = 0.999*ones(1, kern.inputDimension);
kern.nParams = 2 + kern.inputDimension;

kern.transforms(1).index = [1 2];
kern.transforms(1).type = optimiDefaultConstraint('positive');
kern.transforms(2).index = [3:kern.nParams];
kern.transforms(2).type = optimiDefaultConstraint('zeroone');

kern.isStationary = true;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ obj ] = estimate( obj, varargin )
% Estimate parameters of the HUGE model.
%
% INPUTS:
%   obj - A tapas_Huge object containing fMRI time series.
%
% OPTIONAL INPUTS:
%   This function accepts optional name-value pair arguments. For a list of
%   valid name-value pairs, see the user manual or type 'help
%   tapas_huge_property_names'.
%
% OUTPUTS:
%   obj - A tapas_Huge object containing the estimation result in the
%         'posterior' property.
%
% EXAMPLES:
%   [obj] = ESTIMATE(obj)    Invert the HUGE model stored in obj.
%
%   [obj] = ESTIMATE(obj, 'K', 2)    Set the number of clusters to 2 and
%       invert the HUGE model stored in obj.
%
%   [obj] = ESTIMATE(obj, 'Verbose', true)    Print progress of estimation
%       to command line.
%
%   [obj] = ESTIMATE(obj, 'Dcm', dcms, 'OmitFromClustering', 1)    Import
%       data stored in 'dcms' and omit self-connections from clustering.
%
% See also TAPAS_HUGE_DEMO
%

% Author: Yu Yao (yao@biomed.ee.ethz.ch)
% Copyright (C) 2019 Translational Neuromodeling Unit
%                    Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of TAPAS, which is released under the terms of the GNU
% General Public Licence (GPL), version 3. For further details, see
% <https://www.gnu.org/licenses/>.
%
% This software is provided "as is", without warranty of any kind, express
% or implied, including, but not limited to the warranties of
% merchantability, fitness for a particular purpose and non-infringement.
%
% This software is intended for research only. Do not use for clinical
% purpose. Please note that this toolbox is under active development.
% Considerable changes may occur in future releases. For support please
% refer to:
% https://github.com/translationalneuromodeling/tapas/issues
%

%% process and check input
obj = obj.optional_inputs(varargin{:});
assert(obj.N > 0, 'TAPAS:HUGE:NoData', 'Cannot estimate without data.')

tapas_huge_compile( );

% save random number generator seed
seed = rng();

% check if filename for storing result is valid
filename = obj.options.nvp.saveto;
if ~isempty(filename)
    [pathStr, fileStr, extStr] = fileparts(filename);
    if ~strcmpi(extStr,'.mat')
        filename = [filename '.mat'];
        [pathStr, fileStr, extStr] = fileparts(filename);
    end
    if isempty(fileStr)
        fileStr = ['huge' datestr(now,'-yyyymmdd-HHMMSS')];
    end
    filename = fullfile(pathStr, [fileStr, extStr]);
    save(filename, obj);
end

% decide how to treat confounds
switch lower(obj.options.nvp.confoundsvariant)
    case 'default'
        obj.options.confVar = double(obj.M > 0);
    case 'none'
        obj.options.confVar = 0;
    case 'global'
        obj.options.confVar = 1;
    case 'cluster'
        obj.options.confVar = 2;
end
if obj.options.confVar && ~obj.M
    obj.options.confVar = 0;
    warning('TAPAS:HUGE:MissingConfounds',...
        ['HUGE is configured to use group-level confounds, but no values ' ...
         'for confounds were supplied. Proceeding without confounds.']);
end

if obj.K > obj.N
    warning('TAPAS:HUGE:LargeK',...
        'Number of clusters K exceeds number of subjects N.');
end

% check priors
obj.prior = build_prior(obj);

% starting values
obj = prepare_starting_values(obj);

%% model inversion
obj.aux = struct();
fprintf('Inverting HUGE with %u cluster', obj.K);
if obj.K > 1
    fprintf('s');
end
fprintf(' using %s.\n', obj.options.nvp.method);

switch obj.options.nvp.method
    case 'VB'
        obj = obj.vb_invert( );
    case 'MH'
        obj = obj.mh_invert( );
end

%% post processing
obj.aux = [];
obj.options.start = rmfield(obj.options.start, {'subjects', 'clusters'});

% if working with simulated data, calculate balanced purity
obj.posterior.bPurity = [];
if ~isempty(obj.model)
    model = obj.model;
    [~, labels] = max(model.d, [], 2);
    obj.posterior.bPurity = tapas_huge_bpurity(labels, obj.posterior.q_nk);
end
% save additional information
obj.posterior.method = obj.options.nvp.method;
obj.posterior.version = obj.version;
obj.posterior.seed = seed;

% save result
if ~isempty(filename)
    save(filename, obj);
    fprintf('Saved result to "%s".\n', filename);
end

end

function [ prior ] = build_prior( obj )

prior = struct();

% alpha_0
prior.alpha_0 = obj.options.prior.alpha_0;
if isscalar(prior.alpha_0)
    prior.alpha_0 = repmat(prior.alpha_0, obj.K, 1);
end
assert(size(prior.alpha_0, 1) == obj.K && size(prior.alpha_0, 2) == 1 && ...
    all(prior.alpha_0 > 0), 'TAPAS:HUGE:Prior', ...
    'Prior alpha0 must be positive column vector of length %u.', obj.K);

% S_0
prior.S_0 = obj.options.nvp.priorclustervariance;
if isscalar(prior.S_0)
    prior.S_0 = eye(obj.idx.P_c)*prior.S_0;
end
assert(size(prior.S_0, 1) == obj.idx.P_c, 'TAPAS:HUGE:Prior',...
     'PriorClusterVariance must be of size %ux%u.', obj.idx.P_c, obj.idx.P_c);

% m_0
prior.m_0 = obj.options.nvp.priorclustermean;
if isscalar(prior.m_0)
    prior.m_0 = repmat(prior.m_0, 1, obj.idx.P_c);
end
assert(numel(prior.m_0) == obj.idx.P_c, 'TAPAS:HUGE:Prior', ...
       'PriorClusterMean must be row vector of length %u.', obj.idx.P_c);
prior.m_0 = prior.m_0(:)';

% nu_0
prior.nu_0 = obj.options.nvp.priordegree;
assert(all(prior.nu_0(:) > 0), 'TAPAS:HUGE:Prior', ...
    'PriorDegree must be positive scalar.');

% tau_0
prior.tau_0 = obj.options.nvp.priorvarianceratio;
assert(all(prior.tau_0(:) > 0), 'TAPAS:HUGE:Prior', ...
    'PriorCovarianceRatio must be positive scalar.');

% mu_h
prior.mu_h = obj.options.prior.mu_h;
if isscalar(prior.mu_h)
    prior.mu_h = repmat(prior.mu_h, 1, obj.idx.P_h);
end
assert(size(prior.mu_h, 2) == obj.idx.P_h && ...
       size(prior.mu_h, 1) == 1, 'TAPAS:HUGE:Prior', ...
    'Prior homogenous mean must be row vector of length %u.', obj.idx.P_h);

% Sigma_h
prior.Sigma_h = obj.options.prior.Sigma_h;
if isscalar(prior.Sigma_h)
    prior.Sigma_h = eye(obj.idx.P_h)*prior.Sigma_h;
elseif isvector(prior.Sigma_h)
    assert(length(prior.Sigma_h)>=4, 'TAPAS:HUGE:Prior', ...
        'Prior homogenous covariance must be a vector of length 4.')
    A = double(obj.dcm.a - diag(diag(obj.dcm.a)))*prior.Sigma_h(2) + ...
        eye(obj.R)*prior.Sigma_h(1);
    tmp = double([obj.dcm.b(:);obj.dcm.c(:);obj.dcm.d(:)])*prior.Sigma_h(3);
    tmp = [A(:); tmp; ones(2*obj.R+1,1)*prior.Sigma_h(4)];
    prior.Sigma_h = diag(tmp(obj.idx.homogenous));
end
try
    assert(size(prior.Sigma_h, 1) == obj.idx.P_h && ...
        issymmetric(prior.Sigma_h), '', '');
    chol(prior.Sigma_h);
catch
    error('TAPAS:HUGE:Prior', ['Prior hemodynamic covariance must ' ...
        'be symmetric, positive definite matrix of size %u.'], obj.idx.P_h);
end

% mean and variance of lambda
prior.mu_lambda = obj.options.prior.mu_lambda;
assert(prior.mu_lambda > 0, 'TAPAS:HUGE:Prior', ...
    'Prior mean of noise precision must be positive scalar.');
prior.s2_lambda = obj.options.prior.s2_lambda;
assert(prior.s2_lambda > 0, 'TAPAS:HUGE:Prior', ...
    'Prior variance of noise precision must be positive scalar.');

% m_beta_0
prior.m_beta_0 = obj.options.prior.m_beta_0;
if isscalar(prior.m_beta_0)
    prior.m_beta_0 = repmat(prior.m_beta_0, obj.M, 1);
end
assert(size(prior.m_beta_0, 1) == obj.M && ...
       size(prior.m_beta_0, 2) == 1, 'TAPAS:HUGE:Prior', ...
    'Prior confound mean must be column vector of length %u.', obj.M);

% S_beta_0
prior.S_beta_0 = obj.options.prior.S_beta_0;
if isscalar(prior.S_beta_0)
    prior.S_beta_0 = eye(obj.M)*prior.S_beta_0;
end
try
    assert(size(prior.S_beta_0, 1) == obj.M && ...
        issymmetric(prior.S_beta_0), '', '');
    chol(prior.S_beta_0);
catch
    error('TAPAS:HUGE:Prior', ['Prior confound covariance must ' ...
        'be symmetric, positive definite matrix of size %u.'], obj.M);
end

end

function [ obj ] = prepare_starting_values( obj )
% subject-level
% initialize starting values to prior mean
stVal1 = repmat([obj.prior.m_0, obj.prior.mu_h], obj.N, 1);

if ~isempty(obj.options.nvp.startingvaluedcm) && ...
        isnumeric(obj.options.nvp.startingvaluedcm) % custom starting values
    stVal1 = obj.options.nvp.startingvaluedcm;
    [stN, stP] = size(stVal1);
    assert(stN == obj.N && stP == obj.idx.P_c + obj.idx.P_h,...
        'TAPAS:HUGE:StartingValues',...
        'Size of subject-level starting values does not agree with model.')
elseif strcmpi(obj.options.nvp.startingvaluedcm, 'spm') % use SPM estimates
    for n = 1:obj.N
        try
            tmp = [obj.data(n).spm.A(:); obj.data(n).spm.B(:); ...
                obj.data(n).spm.C(:); obj.data(n).spm.D(:); ...
                obj.data(n).spm.transit(:); obj.data(n).spm.decay(:);...
                obj.data(n).spm.epsilon(:)];
            stVal1(n,:) = full(tmp([obj.idx.clustering; obj.idx.homogenous]));
        catch
            fprintf('Missing or invalid SPM estimate for subject %u.' , n);
            fprintf(' Using prior as starting value instead.\n');
        end
    end
end

% randomize starting values
if obj.options.nvp.randomize && ~(strcmpi(obj.options.nvp.method, 'MH') && ...
        ~obj.options.mh.nSteps.dcm) %%% TODO remove exception and introduce new method for clustering only
    stVal1 = stVal1 + randn(size(stVal1)).*obj.options.start.dcm;
end

obj.options.start.subjects = stVal1;

% cluster-level
if ~isempty(obj.options.nvp.startingvaluegmm) && ...
        isnumeric(obj.options.nvp.startingvaluegmm) % custom starting values
    stVal2 = obj.options.nvp.startingvaluegmm;
    [stN, stP] = size(stVal2);
    assert(stN == obj.K && stP == obj.idx.P_c,...
        'TAPAS:HUGE:StartingValues',...
        'Size of cluster-level starting values does not agree with model.')
else
    stVal2 = repmat(obj.prior.m_0, obj.K, 1);
end

% randomize starting values
if obj.options.nvp.randomize && ~(strcmpi(obj.options.nvp.method, 'MH') && ...
        ~obj.options.mh.nSteps.clusters)
    stVal2 = stVal2 + randn(size(stVal2)).*obj.options.start.gmm;
end

obj.options.start.clusters = stVal2;

end

author: translationalneuromodeling

repo: tapas


text:
function [L_t, S_t, iters, frob_err] = ncrpca(M, true_r, EPS, MAX_ITER, EPS_S, incoh, TOL)

% This matlab code implements Non-convex Robust PCA (NcRPCA)
% Input:
% M = given low rank+sparse matrix to be decomposed
% true_r = maximum rank of the low rank rank component
% EPS (optional) = convergence threshold for ||M-(L_t+S_t)||_F; default is 1e-3
% MAX_ITER (optional) = maximum iterations for NcRPCA; default is 51
% EPS_S (optional) = threshold for removing small entries in the sparse component; default is 1e-3
% incoh (optional) = incoherence of the low rank component; default is 1
% TOL (optional) = tolerance for relative error in ||M-(L_t+S_t)||_F in consecutive iterations; default is 1e-1
% Output:
% M_t = thresholded M at each iteration
% L_t = rank-k approximation of M_t
% S_t = sparse component, computed as M-M_t
% iters = number of iteration of NcRPCA
% frob_err = ||M-(L_t+S_t)||_F at each iteration

if nargin < 7, TOL = 1e-1; end
if nargin < 6, incoh = 1; end
if nargin < 5, EPS_S = 1e-3; end
if nargin < 4, MAX_ITER = 51; end
if nargin < 3, EPS = 1e-3; end

%addpath code_ncrpca/PROPACK;
%addpath PROPACK;
frob_err(1) = inf;
[~, n] = size(M);
t = 1;
idx = [];
thresh_const = incoh; % threshold constant: can be tuned depending on incoherence
thresh_red = 0.9; % parameter to reduce the threshold constant: can be tuned
r_hat = 1; % initial rank for stagewise algorithm
L_t = zeros(size(M));
Sig_t = lansvd(M,1,'L');
D_t = M-L_t;
thresh = thresh_const*Sig_t/sqrt(n);
idx = unique([find(abs(D_t) > thresh); idx]);
S_t = zeros(size(M));
S_t(idx) = D_t(idx); % initial thresholding
if max(idx(:))==0
    idx = [];
end
while frob_err(t)/norm(M, 'fro')>=EPS && t<MAX_ITER % convergence check
    if ~mod(t, 10)
        %fprintf('Iter no. %d\n', t);
    end
    t = t+1;
    [U_t, Sig_t, V_t] = lansvd(M-S_t, r_hat+1, 'L');
    %[U_t, Sig_t, V_t] = svds(M-S_t, r_hat+1); % use this if not using propack
    L_t= U_t(:,1:r_hat)*Sig_t(1:r_hat,1:r_hat)*V_t(:,1:r_hat)';
    D_t = M-L_t;
    thresh = (thresh_const/sqrt(n))*Sig_t(r_hat+1, r_hat+1); % use n instead of sqrt(n) for thresholding less aggresively
    idx = unique([find(abs(D_t) > thresh); idx]);
    S_t(idx) = D_t(idx);
    frob_err(t) = norm(M-(L_t+S_t), 'fro');
    if ((frob_err(t-1)-frob_err(t))/frob_err(t-1) <= TOL) && r_hat<true_r
%         r_hat = r_hat+1; % use this for incrementally updating rank by 1
        sig_t = lansvd(M-S_t, true_r, 'L'); % svd function from propack
        ratio_sig = sig_t(r_hat+1:end)./[sig_t(r_hat+2:end); sig_t(end)];
        [~, mx_idx] = max(ratio_sig);
        r_hat = r_hat+mx_idx; % update rank for the next stage
    elseif ((frob_err(t-1)-frob_err(t))/frob_err(t-1) <= TOL) && r_hat==true_r
        thresh_const = thresh_const*thresh_red; % tune threshold
    end
end
S_t(abs(S_t)<EPS_S) = 0; % threshold to remove small errors and obtain sparse component
iters = length(frob_err)-1; % no. of iters. of ncrpca = length(frob_err)-1
end

author: andrewssobral

repo: lrslibrary


text:
function [x,y,button, key] = PointInput(varargin)
%function [x,y,button, key] = PointInput(n,figh, axish)
% Input:
% n = number of clicks
% figh - which figure to analyze the clicks in
% axish - in which axes .. try to avoid this complications, use defaults
% Output:
% x,y - coordinates of the axes where the click happened,
% button - number of button presseed (1 -left,2 -middle, 3 -right
% if key not button was pressed then button=0, and key returns the key that
% was pressed
% also if it is not a character but a system key (like arrow)
% you can get double(key) and use my table from ~antsiro/matlab/keycodes
% file . enjoy :))
% Anton
[n, figh, axish] = DefaultArgs(varargin,{1,gcf,gca});
x=[];y=[];key={};
if n>0
    for i=1:n
        res = waitforbuttonpress;
        if res==0
            whatbutton = get(figh,'SelectionType');
            mousecoord = get(axish,'CurrentPoint');
            x(i)=mousecoord(1,1);
            y(i) = mousecoord(1,2);

            switch whatbutton
                case 'normal'
                    button(i)=1;
                case 'extend'
                    button(i)=2;
                case 'alt'
                    button(i)=3;
            end
            key{end+1}= '';
        else
            x=NaN; y=NaN; button=0;
            key{end+1} = get(gcf,'CurrentCharacter');
        end
    end
    if n==1
        key=key{1};
    end
else
    whatbutton = get(figh,'SelectionType');
    mousecoord = get(axish,'CurrentPoint');
    x= mousecoord(1,1);
    y = mousecoord(1,2);

    switch whatbutton
        case 'normal'
            button=1;
        case 'extend'
            button=2;
        case 'alt'
            button=3;
        otherwise
   %         fprintf('you pressed button %s\n',whatbutton);
            button = 0;
    end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function hu = imHuInvariants(img)
% Compute Hu's invariant for a 2D image.
%
%   HU = imHuInvariants(IMG)
%   HU is a row vector with 12 elements.
%
%   Example
%   imHuInvariants
%
%   See also
%
%   Current implementation is based on the paper:
%   "Noise tolerance of moment invariants in pattern recognition"
%   I. Baslev, Pattern Recognition Letters 19 (1998): 1183-1189
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2008-10-08,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

%% Compute image moments

% total mass of image
m00 = sum(img(:));

% center of mass (first order non-centered moment)
cx = imMoment(img, 1, 0)/m00;
cy = imMoment(img, 0, 1)/m00;

% second-order centered and scaled moments
m02 = imCSMoment(img, 0, 2, [cx cy], m00);
m11 = imCSMoment(img, 1, 1, [cx cy], m00);
m20 = imCSMoment(img, 2, 0, [cx cy], m00);

% third-order centered and scaled moments
m30 = imCSMoment(img, 3, 0, [cx cy], m00);
m21 = imCSMoment(img, 2, 1, [cx cy], m00);
m12 = imCSMoment(img, 1, 2, [cx cy], m00);
m03 = imCSMoment(img, 0, 3, [cx cy], m00);

% fourth-order centered and scaled moments
m40 = imCSMoment(img, 4, 0, [cx cy], m00);
m31 = imCSMoment(img, 3, 1, [cx cy], m00);
m22 = imCSMoment(img, 2, 2, [cx cy], m00);
m13 = imCSMoment(img, 1, 3, [cx cy], m00);
m04 = imCSMoment(img, 0, 4, [cx cy], m00);

%% computation shortcuts

% degree 2
a40 = m20 + m02;
a42 = m20 - m02;
b42 = 2*m11;

% degree 3
a51 = m30 + m12;
b51 = m21 + m03;
a53 = m30 - m12;
b53 = 3*m21 - m03;

% degree 4
a60 = m40 + 2*m22 + m04;
a62 = m40 - m04;
b62 = 2*(m31 + m13);
a64 = m40 - 6*m22 + m04;
b64 = 4*(m31 - m13);

%% Compute Hu's invariants

% init size
hu = zeros(1, 12);

% degree 2
hu(1)   = a40;
hu(2)   = a42^2 + b42^2;
hu(3)   = a42*(a51^2-b51^2) + 2*a51*b51*b42;

% degree 3
hu(4)   = a51^2 + b51^2;
hu(5)   = a53^2 + b53^2;
hu(6)   = a51*a53*(a51^2 - 3*b51^2) + b51*b53*(3*a51^2 - b51^2);
hu(7)   = a51*b53*(a51^2 - 3*b51^2) - b51*a53*(3*a51^2 - b51^2);

% degree 4
hu(8)   = a60;
hu(9)   = a62^2 + b62^2;
hu(10)  = a64^2 + b64^2;
hu(11)  = a64*(a62^2-b62^2) + 2*a62*b62*b64;
hu(12)  = b64*(a62^2-b62^2) - 2*a62*b62*b64;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
% HYPOTHESIS function.
% It predicts the output values y based on the input values X and model parameters.
function [predictions] = hypothesis(X, theta)
    % Input:
    % X - input features - (m x n) matrix.
    % theta - our model parameters - (n x 1) vector.
    %
    % Output:
    % predictions - output values that a calculated based on model parameters - (m x 1) vector.
    %
    % Where:
    % m - number of training examples,
    % n - number of features.

    predictions = sigmoid(X * theta);
end

author: trekhleb

repo: machine-learning-octave


text:
function [H,L,Zf1,Zf0] = nciab(x,BL)
% NCIAB(X) takes data in vector x, filters it with half-band allpass
% filter derived from 5th-order Butterworth filter (ref. Mitra P10.40)
% and returns sub-band coefficients for high-pass channel H and low-pass
% channel L.  BL is buffer length to frame data to transmit final all-
% pass filter conditions.  Zf1 and Zf0 are final condition vectors from
% outputs of analysis filters.
%
% Evan Ruzanski, 4/25/2003

% Set params
M = 2; % DS factor

% Prepare input data
while mod(length(x),2*BL) ~= 0 % Make input even length
    x = [x 0];
end

% Filters (5th order Butterworth per Problem 10.40 Mitra (2e))
numa0 = [0.1056 1];
dena0 = [1 0.1056];

numa1 = [0.5278 1];
dena1 = [1 0.5278];

% Create signals (ref. Fig. 10.55 Mitra (2e))
x0 = [x];
x1 = [x];

% Downsample x0
v0 = x0(1:M:length(x0));

% Downsample x1
v1 = x1(2:M:length(x1)); % Delay = 1

% Filter each channel
lenv = length(v0);
count = lenv/BL;
pointer = 1;
win = [pointer:pointer + BL - 1]; % Window frames of data

bcf = 1; % Init buffer switches
bcb = 4;

while count > 0
    buf0(bcf,1:BL) = v0(win); % Use recursive double buffers
    buf1(bcf,1:BL) = v1(win);

    [u01,Zf0(1,count)] = filter(numa0,dena0,buf0(bcf,1:BL));
    [u11,Zf1(1,count)] = filter(numa1,dena1,buf1(bcf,1:BL));

    buf0(bcb,1:BL) = u01;
    buf1(bcb,1:BL) = u11;

    u0p(win) = buf0(bcb,1:BL);
    u1p(win) = buf1(bcb,1:BL);

    bcf = mod(bcf+2,2)+1; % Switch input buffer
    bcb = mod(bcb+2,2)+3; % Switch output buffer
    pointer = pointer + BL;
    win = [pointer:pointer + BL - 1]; % Slide window
    count = count - 1;
end

% Create/tx band coefficients
L = u0p + u1p;
H = u0p - u1p;

author: Sable

repo: mcbench-benchmarks


text:
function [A,x,y,I]=geoimread_lcg_matlab2020a(filename,varargin)
%GEOIMREAD reads a sub region of a geotiff or geojp2 image.
%
%
%% Syntax
%
% A = geoimread(filename)
% A = geoimread(filename,xlim,ylim)
% A = geoimread(filename,latlim,lonlim)
% A = geoimread(...,buffer)
% [A,x,y,I] = geoimread(...)
% geoimread(...)
%
%
%% Description
%
% A = geoimread(filename) returns the full image given by a filename. This
% syntax is equivalent to A = geotiffread(filename).
%
% A = geoimread(filename,xlim,ylim) limits the region of the geotiff file to
% the limits given by xlim and ylim, which are map units (usually meters) relative
% to the data projection. For example, if the geotiff is projected in Texas Centric
% Mapping System/Lambert Conformal coordinates, xlim and ylim will have units of
% meters relative to the origin (100 W, 18 N). xlim and ylim can be multimensional,
% in which case the limits of the map will be taken as the limits of the limits of
% the distribution of all points in xlim, ylim.
%
% A = geoimread(filename,latlim,lonlim) if no values in xlim, ylim exceed
% normal values of latitudes and longitudes, geoimread assumes you've entered
% limits in geographic coordinates of degrees latitude and longitude. The first
% input is latitude, the second input is longitude.
%
% A = geoimread(...,buffer) adds a buffer in map units (usually meters or feet) to the
% limits of the region of interest.  This may be useful if you want to load an image
% surrounding scattered lat/lon data.  If you'd like an extra 2 kilometers of image around
% your data, enter 2000 as the buffer.  If buffer is a two-element vector, the first
% element is applied to the left and right extents of the image, and the second element
% is applied to the top and bottom extents of the image.
%
% [A,x,y,I] = geoimread(...) also returns pixel center coordinates (x,y) of the
% output image and a geotiff info structure I. I is a useful input for projfwd and projinv.
%
% geoimread(...) without any outputs shows the output image A without loading
% any data into the workspace.
%
%
%% Examples:
%
% % Show a whole geotiff:
% geoimread('boston.tif');
%
% % Compare results from above to a subset geotiff:
% mapx = [765884 766035 766963]; % units are feet
% mapy = [2959218 2957723 2958972];
% geoimread('boston.tif',mapx,mapy)
%
% % Or if you have coordinates in lat/lon and you want a 500 foot buffer:
% lat = [42.3675288 42.3634246 42.3668397];
% lon = [-71.0940009 -71.0934685 -71.0900125];
%
% geoimread('boston.tif',lat,lon,500);
%
%% Author Info:
%
% (c) Aslak Grinsted 2014- (http://www.glaciology.net/)
%     & Chad A. Greene (http://chadagreene.com/)
%
%%
%
% See also GEOTIFFREAD, GEOTIFFINFO, PIXCENTERS, and PROJFWD.

%%  CHAD'S CHANGES:
% The following changes were made by Chad A. Greene (http://chadagreene.com/)
% of the University of Texas Institute for Geophysics (UTIG) on Sept 17, 2014:
%
% * More input checking and error messages.
%
% * Clarified syntax and description in header.
%
% * The fileparts line now only writes the file extension because other outputs went unused.
%
% * If geographic coordinate limits are inferred, they are now ordered latlim,lonlim. <-- **FUNCTIONALITY CHANGE**
%
% * Limits xlim, ylim or latlim, lonlim can be scalar, vector, or matrix-- xlim is now taken
%   as xlim = [min(xlim(:)) max(xlim(:))]. This will save a small step if you have some data points
%   given by x,y.  Now you can simply enter your data coordinates and geoimread will
%   figure out the limits.
%
% * Output variable I now has correct corner coordinates for the subset image instead of
%   the previous version which returned corner coordinates of the full original image.
%
% * A buffer can be specified around input limits.
%
%
% Syntax before the changes:
% geoimread('myimage.tif',[min(easting)-buffer_m max(easting)+buffer_m],...
%     [min(northing)-buffer_m max(northing)+buffer_m]);
%
% Syntax after the changes:
% geoimread('myimage.tif',easting,northing,buffer_m)
%
%

%TODO: support downsampling (ReductionLevel parameter in imread)
%TODO: use map2pix and latlon2pix instead of projfwd and pixcenters. more robust if it is a rotational coordinate system.

%% Set defaults:

usegeocoordinates = false;
returnNONsubsetImage = true;
buffer_x = 0;
buffer_y = 0;

%% Input checks:

% Check for mapping toolbox:
% assert(license('test','map_toolbox')==1,'geoimread requires Matlab''s Mapping Toolbox.')

% Check file type:
assert(isnumeric(filename)==0,'Input filename must be a string.')
[~,~,ext] = fileparts(filename);
switch upper(ext)
    case {'.JP2' '.JPEG2000' '.GEOJP2'}
        I = jp2tiffinfo(filename);
    case {'.TIF' '.TIFF' '.GTIF' '.GTIFF'}
        I = robustgeotiffinfo(filename);
    otherwise
        error('Unrecognized image file type. Must be tif, tiff, gtif, gtiff, jp2, jpeg2000, or geojp2.')
end

% Parse optional inputs:
if nargin>1
    returnNONsubsetImage = false;
    assert(nargin>2,'If you specify an xlim or latlim, you must specify a corresponding ylim or lonlim.')

    % Parse limits:
    xlimOrLatlim = varargin{1}(:);
    ylimOrLonlim = varargin{2}(:);

    assert(isnumeric(xlimOrLatlim)==1,'Input xlim or latlim must be numeric.')
    assert(isnumeric(ylimOrLonlim)==1,'Input ylim or lonlim must be numeric.')

    % Assume geo coordinates if no input limits exceed normal lat/lon values:
    if max(abs(xlimOrLatlim))<=90 && max(abs(ylimOrLonlim))<=360
        usegeocoordinates = true;
    end

    % Parse buffer:
    if nargin>3
        buffer_m = varargin{3};
        assert(isnumeric(buffer_m)==1,'Buffer value must be either a scalar or a two-element array.')
        assert(numel(buffer_m)<3,'Buffer value must be either a scalar or a two-element array.')
        buffer_x = buffer_m(1);
        if isscalar(buffer_m)
            buffer_y = buffer_m(1);
        else
            buffer_y = buffer_m(2);
        end
    end

    if nargin>4
        error('Too many inputs in geoimread.')
    end

end

%% Begin work:

% Get pixel coordinates of full (non-subset) image:
[x,y]=robustpixcenters(I);

if nargout==0
    fprintf('No outputs specified.\n')
    fprintf('X-range: %.1f:%.1f:%.1f\n',x(1),x(2)-x(1),x(end))
    fprintf('Y-range: %.1f:%.1f:%.1f\n',y(1),y(2)-y(1),y(end))
    clear A x y I
    return
end

% Set xlim and ylim depending on user inputs:
if returnNONsubsetImage
    xlimOrLatlim = x(:);
    ylimOrLonlim = y(:);
end

if usegeocoordinates
    % lat/lon limits switch to x/y limits here:
    if ~strcmp(I.ModelType,'ModelTypeGeographic')
        assert(license('test','map_toolbox')==1,'Mapping toolbox needed to project between lat/lon limits and x,y limits. Specify limits in x,y coordinates.')
        [xlimOrLatlim,ylimOrLonlim]=projfwd(I,xlimOrLatlim,ylimOrLonlim);
    end
end

xlim = [min(xlimOrLatlim)-buffer_x max(xlimOrLatlim)+buffer_x];
ylim = [min(ylimOrLonlim)-buffer_y max(ylimOrLonlim)+buffer_y];

% Rows and columns of pixels to read:
rows=find((y>=ylim(1))&(y<=ylim(2)));
cols=find((x>=xlim(1))&(x<=xlim(2)));

%% Display messages if region of interest is partly or wholly outside the image region:

if xlim(1)<min(x)||xlim(2)>max(x)
    disp('geoimread limits extend beyond the available image output in the x direction.')
end

if ylim(1)<min(y)||ylim(2)>max(y)
    disp('geoimread limits extend beyond the available image output in the y direction.')
end

if isempty(rows)||isempty(cols)
    error('No image coordinates can be found inside the specified limits.')
end

%% Load region of interest:
reductionlevel=0;
if reductionlevel==0
    rows=sort(rows([1 end]));
    cols=sort(cols([1 end]));
else
    %% Load region of interest:
    dpx=2^reductionlevel;
    rows=round(rows/dpx);cols=round(cols/dpx);

    rows=sort(rows([1 end]));
    cols=sort(cols([1 end]));
end
x=x(cols(1):cols(end));
y=y(rows(1):rows(end));

A=imread(filename,'PixelRegion',{rows cols});

%% Update info structure to more accurately reflect the new image:

if nargout == 4
    I.FileSize = numel(A);
    I.Height = size(A,1);
    I.Width = size(A,2);
    try
        I.TiePoints.WorldPoints.X = x(1);
        I.TiePoints.WorldPoints.Y = y(1);
        I.SpatialRef.RasterSize = [size(A,1),size(A,2)];
        I.RefMatrix(3,1) = x(1);
        I.RefMatrix(3,2) = y(1);
        I.BoundingBox = [min(x) min(y); max(x) max(y)];
        I.CornerCoords.X = [min(x) max(x) max(x) min(x)];
        I.CornerCoords.Y = [max(y) max(y) min(y) min(y)];
        %TODO: check whether GTRasterTypeGeoKey is RasterPixelIsArea or RasterPixelIsPoint
        I.CornerCoords.Row = .5 + [0 0 size(A,1) size(A,1)]; %TODO: is this .5 always true?
        I.CornerCoords.Col = .5 + [0 size(A,2) size(A,2) 0];
        [I.CornerCoords.Lat,I.CornerCoords.Lon] = projinv(I,I.CornerCoords.X,I.CornerCoords.Y);
        I.GeoTIFFTags.ModelTiepointTag(4) = x(1);
        I.GeoTIFFTags.ModelTiepointTag(5) = y(1);
        I.SpatialRef.XLimWorld = [min(x),max(x)];
        I.SpatialRef.YLimWorld = [min(y),max(y)];
    catch,end
end

%% Clean up:

function I=jp2tiffinfo(fname)
%
% This is a function that attempts to extract the geo-info from a jp2 file.
% There are two "standards" for how this can be stored:
% * GeoJP2: an embedded geotiff is stored in a "uuid" box with a specific id.
% * GMLJP2: an embedded gml file is in an asoc box.
% (extremely annoying that there are multiple "standards".)
%

%Documents:
% JP2 standard: http://www.jpeg.org/public/15444-1annexi.pdf
% GEOJP2 standard: http://wiki.opf-labs.org/download/attachments/11337762/15444-1annexi.pdf
% GMLJP2 standard: http://www.opengeospatial.org/standards/gmljp2

fid=fopen(fname,'r','ieee-be');
lblcontents='';
while ~feof(fid)
    lbox=fread(fid,1,'uint32');
    type=fread(fid,[1 4],'uint8=>char');
    switch lbox
        case 1
            lbox=fread(fid,1,'uint64');lbox=lbox-16;
        case 2:7 %"reserved reserved for ISO use" (????)
            lbox=8;%No idea why ... there is no info on how to interpret values 2-7.
            %Setting this to 8 works for sentinel-2 images.
        otherwise
            lbox=lbox-8;
    end

    switch type
        case 'asoc'
            %the asoc box is a container with other boxes. Set lbox to zero
            %so that we also parse the containing boxes.
            lbox=0;
        case 'lbl '
            %in gmljp2 the xml box is preceded by a lbl box. Read the
            %lbl box contents for later...
            lblcontents=fread(fid,[1 lbox],'uint8=>char');
            lbox=0;
        case 'xml '
            if (strfind(lblcontents,'gml')>0) %gmljp2:
                gmlcontents=fread(fid,[1 lbox],'uint8=>char');
                fout=fopen('.temp.xml','w');
                fwrite(fout,gmlcontents);
                fclose(fout);
                X=xmlread('.temp.xml');
                delete('.temp.xml');

                Elements=X.getElementsByTagName('gml:origin');
                origin=str2num(Elements.item(0).getTextContent());
                Elements=X.getElementsByTagName('gml:offsetVector');
                R=[];
                for ii=1:Elements.getLength
                    R{ii}=str2num(Elements.item(ii-1).getTextContent());
                end
                R=cell2mat(R');
                %ModelPixelScaleTag: [15 15 0]
                %ModelTiepointTag: [0 0 0 4.695e+05 7240500 0]
                assert((R(1,2)==0)&(R(2,1)==0),'unexpected offsetvector in GML');
                I=imfinfo(fname);
                I.GML=gmlcontents;
                I.ModelPixelScaleTag=[R(1,1) -R(2,2) 0]; %order and sign..
                I.ModelTiepointTag=[0 0 0 origin 0];
                fclose(fid);
                return %we have what we need - return
            end
        case 'uuid' %
            uuid=fread(fid,[1 16],'uint8=>char');
            lbox=lbox-16;
            geo=[177 75 248 189 8 61 75 67 165 174 140 215 213 166 206 3];
            if all(uuid==geo)
                fout=fopen('.temp.tif','w');
                contents=fread(fid,lbox,'uint8');lbox=0;
                fwrite(fout,contents);
                fclose(fout);
                fclose(fid);
                I=geotiffinfo_lcog_matlab2020a('.temp.tif');
                m=imfinfo(fname); % a little silly to use imfinfo when I already have a tag reader
                I.Height=m.Height;
                I.Width=m.Width;
                delete('.temp.tif');

                fclose(fid);
                return %we have what we need - return
            end
    end

    fseek(fid,lbox,0);

end
fclose(fid);

%--- BELOW is to make geoimread work even if no mapping toolbox (but with reduced functionality) ---

function I = robustgeotiffinfo(fname)
if license('test','map_toolbox')
    I=geotiffinfo_lcog_matlab2020a(fname);
else
    I=imfinfo(fname);
    % I.ModelType='ModelTypeProjected';  %TODO: fix
    %     %TODO: generate home-made refmatrix(?)....
    %     if isfield(tags, 'ModelTransformationTag') && numel(tags.ModelTransformationTag) >= 8
    %         geoimread does not work for rotated systems
    %
    %     else %use ModelPixelScaleTag instead
    %         dx =  I.ModelPixelScaleTag(1); dy = -I.ModelPixelScaleTag(2);
    %         x0 = I.ModelTiepointTag(4) - dx * I.ModelTiepointTag(1);
    %         y0 = I.ModelTiepointTag(5) - dy * I.ModelTiepointTag(2);
    %         J = [dx 0; 0 dy];
    %     end
    %     I.RefMatrix=[flipud(J); x0-J(1,1)-J(1,2)];
end

function [x,y]=robustpixcenters(I)
if license('test','map_toolbox') && isfield(I,'GeoTIFFCodes')
    [x,y]=pixcenters(I);
else
    %I have not read documentation... but this only works for rectilinear systems.
    assert(I.ModelPixelScaleTag(3)==0,'unexpected ModelPixelScaleTag format.');
    assert(all(I.ModelTiepointTag(1:3)==0),'unexpected ModelTiepointTag format.');
    x=((0:I.Width-1)-I.ModelTiepointTag(1))*I.ModelPixelScaleTag(1)+I.ModelTiepointTag(4);
    y=((0:I.Height-1)-I.ModelTiepointTag(2))*-I.ModelPixelScaleTag(2)+I.ModelTiepointTag(5);
end

author: GERSL

repo: Fmask

save_path: github-repos/MATLAB/GERSL-Fmask


text:
function outdegs = findangle(boundvec, startidx)

% This function computes the bifurcation angle

npix = prod(size(boundvec));
if (sum(abs(boundvec))==0 | sum(abs(boundvec))== npix )
    outdegs =[];
    return;
end

% degree assignments
R  = npix/8;
dy = [0:1:R, R*ones(1, 2*R-1), R:-1:-R, -R* ones(1,2*R-1), -R:1:-1];
dx = [R*ones(1, R+1), R-1:-1:-R+1, -R*ones(1, 2*R+1), -R+1:1:R-1, R*ones(1, R)];
degs = atan2(dy, dx)*180/pi;
degs(degs<0) = degs(degs<0) + 360;

% find connected region
[labelmap, numlabel] = bwlabel(boundvec, 8);
if (boundvec(1)==1 & boundvec(end)==1)
    degs(labelmap == labelmap(end)) = degs(labelmap == labelmap(end))-360;
    labelmap(labelmap == labelmap(end)) = labelmap(1);
    numlabel= numlabel-1;
end

% angle assign to each region
for k=1:numlabel
    seeds = sort(degs(labelmap==k));
    m(k) = (seeds(1)+seeds(end))/2;
end
m(m<0) = m(m<0)+360;

for k = 1:numlabel-1
    outdegs(k)  = m(k+1)-m(k);
end
outdegs(numlabel)= m(1)-m(numlabel);
outdegs(outdegs<0) = outdegs(outdegs<0)+360;

% circulate the startidx to the first one
if (nargin == 2)
    while (boundvec(startidx)==0)
        startidx = startidx+1;
        if (startidx>npix)
            startidx = 1;
        end
    end
    outdegs = circshift(outdegs, [0, -labelmap(startidx)+1]);
end


text:
%%% Setup
% restoredefaultpath;
close, clear, clc;

%%% Path
addpath(fullfile(pwd,'utils'));
addpath(fullfile(pwd,'libs','PROPACK'));
addpath(genpath(fullfile(pwd,'libs','manopt','manopt')));
addpath(fullfile(pwd,'libs','tensor_toolbox'));
addpath(fullfile(pwd,'libs','lightspeed'));

%%% Params
params.debug = 1;
clc;

%%% Load video
load(fullfile('trafficdb','traffic_patches.mat'));
V = im2double(imgdb{100}); % show_3dvideo(V);

%%% Observed entries
%Omega = ones(size(V));
%Omega = randi([0 1],size(V));
[~,Omega] = subsampling(V, 0.75);
%[~,Omega] = subsampling(V, 0.15);
% show_3dvideo(Omega);

%% Matrix completion
clc;
params.algs_path = 'algs_mc';
params.algs_name = get_algs_name(params);
displog('--- Matrix Completion ---');
for alg = 1:size(params.algs_name,2)
  %alg = 1;
  displog(['Current algorithm: ' params.algs_name(alg).name]);
  params.current_algorithm = alg;

  %%% Load algorithm
  current_alg_name = params.algs_name(params.current_algorithm).name;
  current_alg_path = fullfile(params.algs_path,current_alg_name);
  displog(['Loading algorithm: ' current_alg_name]);
  addpath(genpath(current_alg_path));

  %%% Matrix completion
  displog('Performing matrix completion');
  M = convert_video3d_to_2d(V); % imagesc(M)
  Idx = convert_video3d_to_2d(Omega); % imagesc(Idx)
  M(M == 0) = 1e-3;
  M = M.*Idx; % imagesc(M)
  params_mc.M = M;
  params_mc.Idx = Idx;
  M_hat = run_mc(params_mc); % imagesc(M_hat)

  % Build background model
  M_bg = mean(M_hat,2);
  I_bg = reshape(M_bg,size(V,1),size(V,2));
  if(params.debug)
    clf,imshow(I_bg);
    title(current_alg_name);
    pause(1);
  end

  rmpath(genpath(current_alg_path));
  %break;
end
clear alg seq;

%% Tensor completion
clc;
params.algs_path = 'algs_tc';
params.algs_name = get_algs_name(params);
displog('--- Tensor Completion ---');
for alg = 1:size(params.algs_name,2)
  %alg = 1;
  displog(['Current algorithm: ' params.algs_name(alg).name]);
  params.current_algorithm = alg;

  %%% Load algorithm
  current_alg_name = params.algs_name(params.current_algorithm).name;
  current_alg_path = fullfile(params.algs_path,current_alg_name);
  displog(['Loading algorithm: ' current_alg_name]);
  addpath(genpath(current_alg_path));

  %%% Tensor completion
  displog('Performing tensor completion');
  V(V == 0) = 1e-3;
  V = V.*Omega;
  params_tc.T = V;
  params_tc.Idx = Omega;
  params_tc.blksize = [16 16];
  V_hat = run_tc(params_tc); % imagesc(M_hat)

  % Build background model
  I_bg = mean(V_hat,3);
  if(params.debug)
    clf,imshow(I_bg);
    title(current_alg_name);
    pause(1);
  end

  rmpath(genpath(current_alg_path));
  %break;
end
clear alg seq;


text:
function cx = cscal ( n, ca, cx, incx )

%*****************************************************************************80
%
%% CSCAL scales a complex vector by a constant.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, complex CA, the multiplier.
%
%    Input, complex CX(*), the vector to be scaled.
%
%    Input, integer INCX, the increment between successive entries of CX.
%
%    Output, complex CX(*), the scaled vector.
%
  cx(1:incx:1+(n-1)*incx) = ca * cx(1:incx:1+(n-1)*incx);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [trMatrix,elecStatus] = fns_read_transfer(datafile)
% Note that Data only stores the potential values for nodes, which  belong to
% the head

trMatrix      = hdf5read(datafile, '/forward/transfer_matrix');
elecStatus    = hdf5read(datafile, '/forward/status');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%kgextents 'Generate Geometry Representing Extents Around Data'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros gextents.pane file
%
% Parameters:
% InputFile: i 'Input Data Object', required: 'input data object'
% OutputFile: o 'Output Geometry Object', required: 'output geometry object'
% MultiChoice: color 'Colors', default: 2: 'If 1, the lines will be colored white, else they will use the red, green, cyan color scheme described above'
%    Choices are:
% MultiChoice: type 'Box Type', default: 2: 'Select between one of the two box types'
%    Choices are:
% String: name 'Name ', default: 'extents box': 'The name selected here will identify the extents object in the renderer'
%
% Example: o = kgextents(i, {'i','';'o','';'color',2;'type',2;'name','extents box'})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% gextents - Generate Geometry Representing Extents Around Data
%
%  DESCRIPTION
% The gextents module is used to create a geometry object depicting the
% bounds of the space occupied in R3 of the input data.
%
% If present, the location segment of the input data must be in R3
% (i.e., have X, Y and Z components), while the value portion of the
% input data may be of an arbitrary dimension.  In the instance when
% input data is provided which has no location segment, the coordinates
% used are the indeces into the (w,h,d) portion of the value segment.
% In other words, if the input data set has a value segment of size
% width=32, height=48, depth=24, elements=3, and time=N, and no location
% data is present, the resulting box would have it's two corners
% positioned at (0,0,0) and (31,47,23).
%
% Two different box types may be selected by the user.  The first type
% always results in a rectangular shaped box, regardless of the location
% data.  This box type is called a "min/max" box.  The maximum extents
% of the location (either from the explicit location data or from the
% scheme outlined above) of the data are computed, and a box is
% constructed with two corners at those locations.  The other box type
% may result in a non-rectangular box, depending upon the behavior of
% the explicit location data.
%
% The user may elect to have a "colored" box or a white box.  If we
% consider the three dimensions of the data to correspond to the logical
% variables (u,v,w), then the box will be colored so that red lines
% indicate constant (v,w); green lines indicate constant (u,w) and cyan
% lines indicate constant (u,v).
%
% The user also has the option to show the "sides" of the box as grids,
% by setting one or more of the "axis Min" and "axis Max" toggles to
% True.  For example, setting Vaxis Min to True would display the side
% of the box where v = 0.  While viewing the sides of the box, it may be
% helpful to set the "Show Bounding Box" toggle to False, which will
% remove the overall bounding box and show only the geometry for the
% sides that are currently selected.
%
% Additionally, the user may provide an explicit name (text string) for
% the resulting geometry object.  This name may be available to assist
% in object identification and manipulation in a downstream renderer,
% depending upon which one is being used.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1996,1997 , The Regents of the University of California.  All rights reserved.
%

function varargout = kgextents(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kgextents(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'o', '__output';'color', 2;'type', 2;'name', 'extents box'};
maxval={0,0,0,0,0};
minval={0,0,0,0,0};
istoggle=[0,0,0,0,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','MultiChoice','MultiChoice','String'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'gextents"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%-----------------------------------------------------------------------
% Job saved on 08-Feb-2017 12:38:53 by cfg_util (rev $Rev: 7703 $)
% spm SPM - SPM12 (12.3)
% cfg_basicio BasicIO - Unknown
%-----------------------------------------------------------------------
matlabbatch{1}.spm.meeg.source.headmodel.D = {''};
matlabbatch{1}.spm.meeg.source.headmodel.val = 1;
matlabbatch{1}.spm.meeg.source.headmodel.comment = '';
matlabbatch{1}.spm.meeg.source.headmodel.meshing.meshes.template = 1;
matlabbatch{1}.spm.meeg.source.headmodel.meshing.meshres = 2;
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(1).fidname = 'nas';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(1).specification.select = 'nas';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(2).fidname = 'lpa';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(2).specification.select = 'FIL_CTF_L';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(3).fidname = 'rpa';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.fiducial(3).specification.select = 'FIL_CTF_R';
matlabbatch{1}.spm.meeg.source.headmodel.coregistration.coregspecify.useheadshape = 0;
matlabbatch{1}.spm.meeg.source.headmodel.forward.eeg = 'EEG BEM';
matlabbatch{1}.spm.meeg.source.headmodel.forward.meg = 'Single Shell';
matlabbatch{2}.spm.tools.beamforming.data.dir = '<UNDEFINED>';
matlabbatch{2}.spm.tools.beamforming.data.D(1) = cfg_dep('Head model specification: M/EEG dataset(s) with a forward model', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','D'));
matlabbatch{2}.spm.tools.beamforming.data.val = 1;
matlabbatch{2}.spm.tools.beamforming.data.gradsource = 'inv';
matlabbatch{2}.spm.tools.beamforming.data.space = 'MNI-aligned';
matlabbatch{2}.spm.tools.beamforming.data.overwrite = 0;
matlabbatch{3}.spm.tools.beamforming.sources.BF(1) = cfg_dep('Prepare data: BF.mat file', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{3}.spm.tools.beamforming.sources.reduce_rank = [2 3];
matlabbatch{3}.spm.tools.beamforming.sources.keep3d = 1;
matlabbatch{3}.spm.tools.beamforming.sources.plugin.grid.resolution = 10;
matlabbatch{3}.spm.tools.beamforming.sources.plugin.grid.space = 'MNI template';
matlabbatch{3}.spm.tools.beamforming.sources.visualise = 1;
matlabbatch{4}.spm.tools.beamforming.features.BF(1) = cfg_dep('Define sources: BF.mat file', substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{4}.spm.tools.beamforming.features.whatconditions.all = 1;
matlabbatch{4}.spm.tools.beamforming.features.woi = [-Inf Inf];
matlabbatch{4}.spm.tools.beamforming.features.modality = {'MEG'};
matlabbatch{4}.spm.tools.beamforming.features.fuse = 'no';
matlabbatch{4}.spm.tools.beamforming.features.plugin.cov.foi = [0 100];
matlabbatch{4}.spm.tools.beamforming.features.plugin.cov.taper = 'none';
matlabbatch{4}.spm.tools.beamforming.features.regularisation.manual.lambda = 5;
matlabbatch{4}.spm.tools.beamforming.features.bootstrap = false;
matlabbatch{5}.spm.tools.beamforming.inverse.BF(1) = cfg_dep('Covariance features: BF.mat file', substruct('.','val', '{}',{4}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{5}.spm.tools.beamforming.inverse.plugin.lcmv.orient = true;
matlabbatch{5}.spm.tools.beamforming.inverse.plugin.lcmv.keeplf = false;
matlabbatch{6}.spm.tools.beamforming.output.BF(1) = cfg_dep('Inverse solution: BF.mat file', substruct('.','val', '{}',{5}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_mv.isdesign.custom.whatconditions.all = 1;
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_mv.isdesign.custom.contrast = [-1 1];
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_mv.isdesign.custom.woi = [-100 0
                                                                                100 200];
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_mv.datafeatures = 'sumpower';
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_mv.foi = [5 15
                                                                   15 30
                                                                   30 60
                                                                   60 90];
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_mv.result = 'BIC';
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_mv.sametrials = false;
matlabbatch{6}.spm.tools.beamforming.output.plugin.image_mv.modality = 'MEG';
matlabbatch{7}.spm.tools.beamforming.write.BF(1) = cfg_dep('Output: BF.mat file', substruct('.','val', '{}',{6}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','BF'));
matlabbatch{7}.spm.tools.beamforming.write.plugin.nifti.normalise = 'no';
matlabbatch{7}.spm.tools.beamforming.write.plugin.nifti.space = 'mni';

author: spm


text:
function features = values_to_features(values)
	% Convert a time-series obtained from a sample to a feature using
	% different kind of statistics over the values and the derivatives

	mean_val = nanmean(values, 2);
	variance = nanvar(values, 0, 2);
	% feature_skewness = skewness(values(~isnan(values)));
	% feature_kurtosis = kurtosis(values(~isnan(values)));

	abs_delta = abs(values(:, 2:end) - values(:, 1:end-1));
	mean_delta = nanmean(abs_delta, 2);
	var_delta = nanvar(abs_delta, 0, 2);

	maximum = max(values, [], 2);
	minimum = min(values, [], 2);

	features = [mean_val; variance; ...
	%           feature_skewness; feature_kurtosis; ...
                mean_delta; var_delta; maximum; minimum];
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function val = essential_costE2cost(X, costE)
% Cost evaluation at X given function handle in the Essential matrix E.
%
% function val = essential_costE2cost(X, costE)
%
% costE is the function handle for the cost function in E.
%
% See also: essential_egradE2egrad essential_ehessE2ehess

% This file is part of Manopt: www.manopt.org.
% Original author: Roberto Tron, Aug. 8, 2014
% Contributors: Bamdev Mishra, May 22, 2015.

    e3hat = [0 -1 0; 1 0 0; 0 0 0];

    RA = X(:,1:3,:);
    RB = X(:,4:6,:);
    E = multiprod(multiprod(multitransp(RA), e3hat), RB);

    val = costE(E);
end

author: NicolasBoumal

repo: manopt


text:
function SO3VF = abs(F)
% overloads the componentwise absolute value |abs(SO3VF)|
%
% Note that abs is not the length of a vector3d which is the evaluation of
% the SO3VectorField in one rotation.
%
% one
%
% Syntax
%   SO3VF = abs(SO3VF)
%
% Input
%  SO3VF - @SO3VectorField
%
% Output
%  SO3VF - @SO3VectorField
%

SO3VF = SO3VectorFieldHandle(@(rot) g(rot),F.SRight,F.SLeft);

function g = g(rot)
g = F.eval(rot);
g = vector3d(abs(g.xyz.'));
end

end

author: mtex-toolbox

repo: mtex

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


text:
function inspect_bug3157

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY

%%
ctf151      = ft_read_sens(dccnpath('/home/common/matlab/fieldtrip/data/test/original/meg/ctf151/Subject01.ds'));
neuromag122 = ft_read_sens(dccnpath('/home/common/matlab/fieldtrip/data/test/original/meg/neuromag122/jg_single_01raw.fif'));
neuromag306 = ft_read_sens(dccnpath('/home/common/matlab/fieldtrip/data/test/original/meg/neuromag306/raw.fif'), 'senstype', 'meg');

%%

rx = 30;
ry = 30;
rz = 30;

fake = [];
fake.unit = 'mm';
fake.coilpos = [
  10  0 0
  -10 -0 0
  ];
fake.coilpos = ft_warp_apply(rotate([rx ry rz]), fake.coilpos);
fake.coilori = [
  0 0 1
  0 0 1
  ];
fake.coilori = ft_warp_apply(rotate([rx ry rz]), fake.coilori);
fake.label = {
  '1'
  };
fake.tra = [
  1 -1
  ];

figure
ft_plot_sens(fake, 'coil', false, 'unit', 'mm', 'coilsize', 40, 'coilshape', 'square');
ft_plot_sens(fake, 'coil', true,  'unit', 'mm', 'coilsize', 0);
grid on

%%

figure
ft_plot_sens(ctf151, 'coil', false, 'unit', 'mm', 'coilsize', 15, 'coilshape', 'circle', 'chantype', 'meggrad');

%%

figure
ft_plot_sens(neuromag122, 'coil', false, 'unit', 'mm', 'coilsize', 35, 'coilshape', 'square', 'chantype', 'megplanar');

%%

figure
ft_plot_sens(neuromag306, 'coil', false, 'unit', 'mm', 'coilsize', 30, 'coilshape', 'square', 'chantype', 'megplanar');

%%

figure
ft_plot_sens(neuromag306, 'coil', false, 'unit', 'mm', 'coilsize', 30, 'coilshape', 'square', 'chantype', 'megplanar', 'facecolor', 'b', 'facealpha', 0.1);


text:
disp("Question 1")
A = [1 2; 3 4; 5 6]
B = [1 2 3; 4 5 6]

C = B' + A
#C = A' * B
#C = B + A
C = A * B

disp("Question 2")

A = [16,2,3,13;5,11,10,8;9,7,6,12;4,14,15,1]

#B = A(:, 0:2)
B = A(:, 1:2)
B = A(1:4, 1:2)
#B = A(0:4, 0:2)

disp("Question 3")
A = ones(10,10);
x = ones(10,1);

v = zeros(10, 1);

for i = 1:10
	 for j = 1:10
	   v(i) = v(i) + A(i, j) * x(j);
	 end
end

v

v = zeros(10, 1);
v = A * x

v = zeros(10, 1);
v = sum (A * x)

v = zeros(10, 1)
v = A .* x;

v = zeros(10, 1)
#v = Ax;

disp("Question 4")
v = ones(7,1);
w = ones(7,1);

z = 0;
for i = 1:7
  z = z + v(i) * w(i);
end

z

z = v' * w

#z = v * w

z = sum (v .* w)

z = v .* w

disp("Question 4")
X = ones(7,7)

for i = 1:7
  for j = 1:7
    	A(i, j) = log (X(i, j));
	    B(i, j) = X(i, j) ^ 2;
    	C(i, j) = X(i, j) + 1;
	    D(i, j) = X(i, j) / 4;
   end
end

A
B
C
D

B = X .^ 2

C = X + 1

D = X / 4

B = X ^ 2

author: vugsus

repo: coursera-machine-learning


text:
function [scoreMat_viewpoint] = generateViewPointScoreMatrix(viewPointQ, viewPointT)

    vQ = [cos(viewPointQ) sin(viewPointQ)];
    vT = [cos(viewPointT) sin(viewPointT)];

    scoreMat_viewpoint = vQ * (vT');

end

repo: MOTBeyondPixels


text:
function data = simulate_flight(setup)

%This function simulates the flight phase of the dynamics

%Set up the integration algorithm
Tspan = setup.Tspan;
IC = [  setup.IC.th;
    setup.IC.x;
    setup.IC.y;
    setup.IC.dth;
    setup.IC.dx;
    setup.IC.dy];
eventFunc = @(t,z)events_flight(t,z,setup.P);
options = odeset(...
    'RelTol',setup.tol,...
    'AbsTol',setup.tol,...
    'Vectorized','on',...
    'MaxStep',setup.odeMaxStep,...
    'Events',eventFunc);
userfun = @(t,z)dynamics_flight(t,z,setup.P);

%Run simulation
sol = feval(setup.solver,userfun,Tspan,IC,options);

%Format for post processing
tspan = [sol.x(1),sol.x(end)];
nTime = ceil(setup.dataFreq*diff(tspan));
t = linspace(tspan(1),tspan(2),nTime);
Z = deval(sol,t);
[~, E] = dynamics_flight(t,Z,setup.P);

%Store in a nice format for plotting
data.time = t;
data.state.th = Z(1,:);
data.state.x = Z(2,:);
data.state.y = Z(3,:);
data.state.dth = Z(4,:);
data.state.dx = Z(5,:);
data.state.dy = Z(6,:);
data.contact.h = zeros(1,nTime);
data.contact.v = zeros(1,nTime);
data.energy.potential = E(1,:);
data.energy.kinetic = E(2,:);
data.P = setup.P;

%Get transitions for finite state machine
data.phase = 'FLIGHT';
if isempty(sol.ie)
    data.exit = 'TIMEOUT';
else
    switch sol.ie(end)
        case 1   %close end of the stick hit the ground
            data.exit = 'STRIKE_0';
        case 2   %Far end of the stick hit the ground
            data.exit = 'STRIKE_2';
        otherwise
            error('Invalid exit condition in simulate flight')
    end
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function det = r8cb_det ( n, ml, mu, a_lu )

%*****************************************************************************80
%
%% R8CB_DET computes the determinant of a R8CB matrix factored by R8CB_NP_FA.
%
%  Discussion:
%
%    The R8CB storage format is appropriate for a compact banded matrix.
%    It is assumed that the matrix has lower and upper bandwidths ML and MU,
%    respectively.  The matrix is stored in a way similar to that used
%    by LINPACK and LAPACK for a general banded matrix, except that in
%    this mode, no extra rows are set aside for possible fillin during pivoting.
%    Thus, this storage mode is suitable if you do not intend to factor
%    the matrix, or if you can guarantee that the matrix can be factored
%    without pivoting.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative, and no greater than N-1.
%
%    Input, real A_LU(ML+MU+1,N), the LU factors from R8CB_NP_FA.
%
%    Output, real DET, the determinant of the matrix.
%
  det = prod ( a_lu(mu+1,1:n) );

  return
end


text:
function varargout=gabor3_fwb(aspect,theta,bw,psi,sigma,sz)
% Returns 3D gabor filter.
% gb=GABOR_FWB(aspect,theta,bw,psi,sigma,sz)
%
% [aspecta, aspectb]
%        = 2 element array giving aspect ratios for 2 minor axis
%           (eg: [0.5, 1], for major < minoraxis1, major = minoraxi2)
% [theta,phi]
%        = yaw and pitch of major axis (0-2*pi)
%           roll isn't implemented, sorry.
% bw     = spatial bandwidth in pixels (decreasing fine detail,), (eg: >=1)
%               scales the frequency of the cosine modulation
% psi,   = phase shift, [optional, default: 0]
% sigma  = scales the falloff of the gaussian, (must be >=2) [default: = bw]
%              + can set to 'auto' to maintain default functionality
% [x y z] = size of gabor kernel created  [optional, size set automatically
%           to 3 standard deviations of gaussian kernel]

% Frederick Bryan adapted from gabor_fn.m
% July 2013
% Vanderbilt Univ

% To Do
% implement roll:
% This isn't hard as far as the meshgrid goes. All that needs to be done for that is to
% uncomment the appropriate lines. I just didn't want to figure out the math
% for the auto-sizing option. (if nargin < 6)
%

% handle mandatory args
if numel(aspect) ~= 2;
    error('1st argument (aspect) must be 2 element array');
end
if numel(theta) ~=2;
    error('2nd argument ([theta,phi]) must be 2 element array');
end

% handle optional inputs
if nargin<4
    psi = [0 0];
end
if nargin<5
    sigma = 'auto';
end
% allow 'auto' sizing of guassian kernel
if strcmp(sigma,'auto');
    sigma = bw;
end
if nargin<6
    % figure out size
    phi = theta(2);
    theta = theta(1);
    len1 = 3*sigma; % length along theta direction
    len2 = len1/aspect(1);
    len3 = len1/aspect(2);
    sz(1) = (len1*cos(theta)+len2*sin(theta)); % column/x dimension
    sz(2) = (len1*sin(theta)+len2*cos(theta)); % row/y direction
    sz(1) = round(sz(1)*cos(phi)+len3*sin(phi));
    sz(2) = round(sz(2));
    sz(3) = round(sz(1)*sin(phi)+len3*cos(phi));
end

% handle incorrectly given optional args

if length(sz)<2; % allow just one number to be given for size
    sz(2) = sz(1);
    sz(3) = sz(1);
end
% keyboard;
sx = sz(1);
sy = sz(2);
sz = sz(3);

% figured out size above, now make matrix of points
[x y z]=meshgrid(-sx:sx, sy:-1:-sy, -sz:sz); % note that y goes backwards

% rotate reference frame to point in theta direction
% http://en.wikipedia.org/wiki/Rotation_matrix
% yaw - rotation about z
xp = x*cos(theta)+y*sin(theta);
yp = -x*sin(theta)+y*cos(theta);
zp = z;
% pitch - rotation about y
xp = xp*cos(phi)-zp*sin(phi);
yp = yp;
zp = xp*sin(phi)+zp*cos(phi);
% % roll - rotation about x - unimplemented, requires "roll" angle
% xp = xp;
% yp = yp*cos(roll) + zp*sin(roll);
% zp = -yp*sin(roll) + zp*cos(roll);

% create gaussian pointing in theta direction with size determined by
%   aspect ration
sigmajor = sigma;
sigminor1 = sigma/aspect(1);
sigminor2 = sigma/aspect(2);
h1 = exp(-(xp.^2/sigmajor^2+yp.^2/sigminor1^2+zp.^2/sigminor2^2));

% multiply by cosine with appropriate bw
F=1/bw;  % Frequency
h2 = cos(2*pi*(F*xp)+psi);

g = h1.*h2;

if nargout>0
    varargout{1} = g;
end


text:
function [MPa] = psi2MPa(psi)
% Convert units of pressure from pounds per square inch to megapascals.
% Chad Greene 2012
MPa = psi*0.00689476;


text:
function	[A,W,SW,S,Q,ix_act] = ...
				compo_update_act(XX,A,W,SW,S,Q,SY,i_up,jx,ix_act)
% Update component
% XX = X(ix_act,:)*X'  [M_act x M_all]

% New alpha
a0 = A(i_up);
s0 = S(i_up);
b0 = 1 - s0./a0;	% scale factor for Q & S
qq0  = SY*sum(Q(:,i_up).^2,1);
anew = s0.^2 ./(qq0 - s0.*b0);

%if (anew - a0) < eps, return; end;

A(i_up) = anew;

% Weight  & variance update
ss = 1./(SW(jx,jx) + 1./(anew - a0));
SX = SW(jx,:);
SW = SW - ss * (SX' * SX);
SW = 0.5 * (SW + SW');	% guarantee symmetric matrix
Wj = W(:,jx) * ss;
W  = W - Wj * SX;

% Residual correlation
% XX = X(ix_act,:)*X'  [M_act x M_all]
P  = SX * XX;
S  = abs(S + ss * P.^2);
Q  = Q + Wj * P;
return

author: KamitaniLab

repo: GenericObjectDecoding


text:
% Reference Parameters

% Finish Time and Pen Time Table
[time_finish, time_table_pen] = cal_time_data(cp_ptp, ts1, time_wait);
max_pen_idx = length(time_table_pen) - 1;

% Position Reference
cp_ptp = cell2mat(cp_ptp);
x_ref = cp_ptp(1, :);
y_ref = cp_ptp(2, :);

% Angle Reference
[theta1_ref, theta2_ref] = xy2theta(x_ref, y_ref, l1, l2);

% Theta and Omega Limitation Check
thetas = [theta1_ref; theta2_ref];
thetas_max = [theta1_max; theta2_max];
omega1m_max = (100 - pwm1_offset) / pwm1_gain;
omega2m_max = (100 - pwm2_offset) / pwm2_gain;
omegas_max = [omega1m_max / g1; omega2m_max / g2];
chk_limit(thetas, thetas_max, omegas_max, ts1)

clear cp_ptp
clear thetas thetas_max omegas_max omega1m_max omega2m_max


text:
function [nodes,weights]=gaussquadrule(n,class,alpha,beta)
% gaussquadrule: generate a gauss quadrature rule using associated orthogonal polynomials
% usage: [nodes,weights]=gaussquadrule(n,class,alpha,beta)
%
% arguments:
%        n - number of nodes in the guassian quadrature
%
%    class - (optional) - class of gaussian quadrature rule
%            any one of { 'Legendre' 'Laguerre' 'Hermite',
%            'Jacobi', '1Cheby', '2Cheby' }. '1cheby' refers
%            of course to a first kind Gauss-Chebychev rule.
%
%            Capitalization is ignored and class names
%            may be shortened as long as they are left
%            unambiguous. These are all minimally valid:
%            {'le', 'la', 'h' 'j' '1' '2'}
%
%            DEFAULT == 'Legendre'
%
%            The nominal integration intervals are:
%
%              'Legendre' --> [  -1,  1]
%              'Laguerre' --> [   0,inf]
%              'Hermite'  --> [-inf,inf]
%              'Jacobi'   --> [  -1,  1]
%              '1Cheby'   --> [  -1,  1]
%              '2Cheby'   --> [  -1,  1]
%
%            See Abramowitz & Stegun for in-depth information
%
% alpha, beta - (optional) parameters for the rules
%            DEFAULT == 0 for both
%
%            alpha and beta are ignored for the 'legendre',
%            'hermite', 'cheby1' and 'cheby2' rules. Only
%            alpha is used for the Laguerre quadrature rules.
%
% arguments (output)
%    nodes - x values for integration
%  weights - integration weights

% default parameters
if (nargin<2) || isempty(class)
  class='legendre';
end

class=lower(class);
valc={'legendre' 'laguerre' 'hermite' 'jacobi' 'cheby1' ...
      'cheby2' '1chebychev' '2chebychev'};
ind=strmatch(class,valc);
if isempty(ind)
  error 'Invalid quadrature class'
elseif length(ind)>1
  error 'Ambiguous quadrature class'
else
  class=valc{ind};
end

if (nargin<4) || isempty(beta)
  beta=0;
end
if (nargin<3) || isempty(alpha);
  alpha=0;
end

% generate the appropriate orthogonal sympolys.
% use the sympoly objects for this. orthpoly
% creates them.
pn = orthpoly(n,class,alpha,beta);
pnp1 = orthpoly(n+1,class,alpha,beta);

nodes=roots(pn);
nodes=sort(nodes');

pnprime=diff(pn);
weights=zeros(1,n);
for i=1:n
  weights(i)=-(pnp1.Coefficient(end))/(pn.Coefficient(end))./ ...
    double(subs(pnp1,'x',nodes(i)))./ ...
    double(subs(pnprime,'x',nodes(i)));
end

% norm
switch class
  case 'legendre'
    weights=weights*(2/(2*n+1));
  case 'hermite'
    weights=sqrt(pi)*(2^n)*(gamma(n+1))*weights;
  case 'laguerre'
    weights=weights*gamma(n+alpha+1)/gamma(n+1);
  case 'jacobi'
    weights=weights*(2/sum(weights));
    nodes=(nodes+1)/2;
  case {'cheby1' '1chebychev'}
    weights=weights*(2/sum(weights));
  case {'cheby2' '2chebychev'}
    weights=weights*(2/sum(weights));
end


text:
function subplots(data, varargin)
    figureTitle = '';
    numRows = [];
    numCols = [];
    xAxis = [];
    linkingAxes = 1;

    for i=1:length(varargin)
        if(ischar(varargin{i}) || isscalar(varargin{i}))
            switch(varargin{i})
                case 'Title'
                    figureTitle = varargin{i+1};
                case 'NumRows'
                    numRows = varargin{i+1};
                case 'NumCols'
                    numCols = varargin{i+1};
                case 'LinkingAxes'
                    linkingAxes = varargin{i+1};
                case {'XAxis', 'xAxis'}
                    xAxis = varargin{i+1};
            end
        end
    end

    if(isempty(numCols) && isempty(numRows))
        numRows = size(data,2);
        numCols = 1;
    elseif(~isempty(numCols) && isempty(numRows))
        numRows = ceil(size(data,2)/numCols);
    end

    figure('Name',figureTitle);
    %subplot(numRows, numCols, 1);
    if(isempty(xAxis))
        xAxis = (1:length(data));
    end
    ah = [];
    for i=1:size(data,2)
        ah(i) = subplot(numRows, numCols, i);
        plot(xAxis, data(:,i));
    end
    if(linkingAxes)
        linkaxes(ah);
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%

for i = 1:length(X)
    minDist = Inf;
    for j = 1:K
        dist = norm(X(i,:) - centroids(j,:)) ^ 2;

        if dist < minDist
            minDist = dist;
            idx(i) = j;
        end
    end
end

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

end

repo: MachineLearning


text:
function value = ellipse_circumference_2d ( r1, r2 )

%*****************************************************************************80
%
%% ELLIPSE_CIRCUMFERENCE_2D returns the circumference of an ellipse in 2D.
%
%  Discussion:
%
%    There is no closed formula for the circumference of an ellipse.
%
%    Defining the eccentricity by
%
%      E = sqrt ( 1 - ( r2 / r1 )**2 )
%
%    where R1 and R2 are the major and minor axes, then
%
%      circumference
%        = 4 * R1 * E(K,2*PI)
%        = R1 * Integral ( 0 <= T <= 2*PI ) sqrt ( 1 - E**2 * sin**2 ( T ) ) dT
%
%    This integral can be approximated by the Gauss-Kummer formula.
%
%  Integration region:
%
%    Points (X,Y) such that
%
%      ( ( X - XC ) / R1 )**2 + ( ( Y - YC ) / R2 )**2 <= 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    John Harris and Horst Stocker,
%    Handbook of Mathematics and Computational Science,
%    Springer, 1998.
%
%  Parameters:
%
%    Input, real R1, R2, the major and minor semi-axes.
%
%    Output, real VALUE, the circumference of the ellipse.
%
  if ( r1 == r2 ) then
    value = 2.0 * pi * r1;
    return
  end
%
%  Compute the eccentricity of the ellipse.
%
  e = sqrt ( 1.0 - ( min ( r1, r2 ) / max ( r1, r2 ) )^2 );

  value = 1.0;
  term = value;
  i = 0;

  while ( 1 )

    i = i + 1;
    term = term * ( 2 * i - 3 ) * ( 2 * i - 1 ) * e * e / ( 2 * 2 * i * i );

    if ( abs ( term ) <= eps * ( abs ( value ) + 1.0 ) )
      break
    end

    value = value + term;

  end

  value = 2.0 * pi * max ( r1, r2 ) * value;

  return
end


text:
%
% Usage:   A=mexOMPMask(X,D,B,param);
% or       [A path]=mexOMPMask(X,D,B,param);
%
% Name: mexOMPMask
%
% Description: mexOMPMask is a variant of mexOMP that allow using
%     a binary mask B
%
%     for all columns x of X, and columns beta of B, it computes a column
%         alpha of A by addressing
%         min_{alpha} ||alpha||_0  s.t  ||diag(beta)*(x-Dalpha)||_2^2
%                                                               <= eps*||beta||_0/m
%         or
%         min_{alpha} ||diag(beta)*(x-Dalpha)||_2^2  s.t. ||alpha||_0 <= L
%         or
%         min_{alpha} 0.5||diag(beta)*(x-Dalpha)||_2^2  + lambda||alpha||_0
%
%
% Inputs: X:  double m x n matrix   (input signals)
%            m is the signal size
%            n is the number of signals to decompose
%         D:  double m x p matrix   (dictionary)
%            p is the number of elements in the dictionary
%            All the columns of D should have unit-norm !
%         B:  boolean m x n matrix   (mask)
%               p is the number of elements in the dictionary
%         param: struct
%            param.L (optional, maximum number of elements in each decomposition,
%               min(m,p) by default)
%            param.eps (optional, threshold on the squared l2-norm of the residual,
%               0 by default
%            param.lambda (optional, penalty parameter, 0 by default
%            param.numThreads (optional, number of threads for exploiting
%            multi-core / multi-cpus. By default, it takes the value -1,
%            which automatically selects all the available CPUs/cores).
%
% Output: A: double sparse p x n matrix (output coefficients)
%         path (optional): double dense p x L matrix
%                                     (regularization path of the first signal)
%
% Note: this function admits a few experimental usages, which have not
%     been extensively tested:
%      - single precision setting (even though the output alpha is double
%        precision)
%      - Passing an int32 vector of length n to param.L provides
%        a different parameter L for each input signal x_i
%      - Passing a double vector of length n to param.eps and or param.lambda
%        provides a different parameter eps (or lambda) for each input signal x_i
%
% Author: Julien Mairal, 2010

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
%% tri6_subtri3
% Below is a demonstration of the features of the |tri6_subtri3| function

%%
clear; close all; clc;

%% Syntax
% |[Ft3]=tri6_subtri3(Ft6,Vt6);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function volr = volresize(vol, newsz, varargin)
    % volresize(vol, newsz, {interp_type}, {offset_type}
    interp_type = 1;
    if nargin >= 3
        interp_type = varargin{1};
    end

    use_offset = 'sample_std';
    if nargin >= 4
        use_offset = varargin{2};
    end

    tmp = zeros(max(newsz, size(vol)));
    tmp(1:size(vol,1), 1:size(vol,2), 1:size(vol,3)) = vol;
    [n1, n2, n3] = ndgrid(1:size(tmp,1), 1:size(tmp,2), 1:size(tmp,3));
    k = size(vol) ./ newsz;
%     k = (size(vol)-1) ./ (newsz - 1);

    %yeah
    k = size(vol) ./ newsz;
    if strcmp(use_offset, 'matlab')
        d = 0.5*(1-k);
    elseif strcmp(use_offset, 'for_upsampling')
        d = 0.5 * [1,1,1];
%         d = [0,0,0];
    elseif strcmp(use_offset, 'for_upsampling_2')
        k = (size(vol) - 1) ./ (newsz - 1);
        k = size(vol)./ newsz;
        d = 1 - k;
    elseif strcmp(use_offset, 'sample_std')
        k = (newsz - 2) ./ (size(vol) - 2);
        k = (size(vol) - 1) ./ (newsz - 1);
        d = (1 - k);
     end

    T = cat(4, k(1) * n1 - n1 + d(1), k(2) * n2 - n2 + d(2), k(3) * n3 - n3  + d(3));

%         if any(k < 1)
%             sgm = 0.00001 + max(k - 1, 0)*1.1;
%             tmp = imgaussfilt3(tmp, sgm);
%         end

%     T = cat(4, (k(1)-1) * n1, (k(2)-1) * n2, (k(3)-1) * n3);

    volr = imdeform3(tmp, T, interp_type);

    volr = volr(1:newsz(1), 1:newsz(2), 1:newsz(3));
end

repo: pTVreg


text:
function F = intkernel(X,vp,gp,avg_flag)
%INTKERNEL Expected GP kernel in scalar correlation

if nargin < 4 || isempty(avg_flag); avg_flag = false; end

K = vp.K;           % Number of components
[N,D] = size(X);
mu(:,:) = vp.mu;
sigma(1,:) = vp.sigma;
lambda(:,1) = vp.lambda(:);
w(1,:) = vp.w;

Ns = numel(gp.post);            % Hyperparameter samples

F = zeros(N,Ns);

if isfield(vp,'delta') && ~isempty(vp.delta)
    delta = vp.delta;
else
    delta = 0;
end

% Integrated mean function being used?
integrated_meanfun = isfield(gp,'intmeanfun') && gp.intmeanfun > 0;

if integrated_meanfun
    % Evaluate basis functions
    Hs = gplite_intmeanfun(X,gp.intmeanfun);
end

% Loop over hyperparameter samples
for s = 1:Ns
    hyp = gp.post(s).hyp;

    % Extract GP hyperparameters from HYP
    ell = exp(hyp(1:D));
    ln_sf2 = 2*hyp(D+1);
    sum_lnell = sum(hyp(1:D));

    if integrated_meanfun
        %betabar = gp.post(s).intmean.betabar';
        %KinvHtbetabar = gp.post(s).intmean.HKinv'*betabar;
        plus_idx = gp.intmeanfun_var > 0;
        HKinv = gp.post(s).intmean.HKinv(plus_idx,:);
        Tplusinv = gp.post(s).intmean.Tplusinv;
    end

    L = gp.post(s).L;
    Lchol = gp.post(s).Lchol;

    sn2_eff = 1/gp.post(s).sW(1)^2;

    ddl = sq_dist(bsxfun(@rdivide,X',ell),bsxfun(@rdivide,gp.X',ell));
    ll = exp(ln_sf2 -0.5*ddl);

    if Lchol
        zz = (L\(L'\ll'))/sn2_eff;
    else
        zz = -L*ll';
    end

    for k = 1:K
        tau_k = sqrt(sigma(k)^2*lambda.^2 + ell.^2 + delta.^2);
        lnnf_k = ln_sf2 + sum_lnell - sum(log(tau_k));  % Covariance normalization factor
        delta_k = bsxfun(@rdivide,bsxfun(@minus, mu(:,k), gp.X'), tau_k);
        z_k = exp(lnnf_k -0.5 * sum(delta_k.^2,1));

        dd_k = bsxfun(@rdivide,bsxfun(@minus, mu(:,k), X'), tau_k);
        zz_k = exp(lnnf_k -0.5 * sum(dd_k.^2,1));

        F(:,s) = F(:,s) + w(k)*(zz_k - z_k*zz)';

        % Contribution of integrated mean function
        if integrated_meanfun
            switch gp.intmeanfun
                case 1; u_k = 1;
                case 2; u_k = [1,mu(:,k)'];
                case 3; u_k = [1,mu(:,k)',(mu(:,k).^2 + sigma(k)^2*lambda.^2)'];
                case 4; u_k = [1,mu(:,k)',(mu(:,k).^2 + sigma(k)^2*lambda.^2)',mumu_mat(k,:)];
            end

            F(:,s) = F(:,s) + w(k)*((u_k(plus_idx)*(Tplusinv*Hs)) ...
                + ((z_k*HKinv')*(Tplusinv*(HKinv*ll'))) ...
                - (u_k(plus_idx)*(Tplusinv*(HKinv*ll'))) ...
                - ((z_k*HKinv')*(Tplusinv*Hs)))';
        end

    end
end

% Average multiple hyperparameter samples
if Ns > 1 && avg_flag
    F = mean(F,2);
end

end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc