text:
function [L,D] = spm_eeg_lgainmat(D,Is,channels)
% Load or compute if necessary a gain matrix
% FORMAT [L,D] = spm_eeg_lgainmat(D,Is,channels)
% D    - Data structure
% Is   - indices of vertices
%
% L    - Lead-field or gain matrix L(:,Is)
%__________________________________________________________________________
% Copyright (C) 2008-2017 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_eeg_lgainmat.m 7757 2019-12-16 15:36:06Z spm $

SVNrev = '$Rev: 7757 $';

%-Get gain or lead-field matrix
%--------------------------------------------------------------------------
val = D.val;

forward = D.inv{val}.forward;

for ind = 1:numel(forward)
    modality = forward(ind).modality;

    %-Channels
    %----------------------------------------------------------------------
    if isequal(modality, 'MEG')
        chanind = D.indchantype({'MEG', 'MEGPLANAR'}, 'GOOD');
    else
        chanind = D.indchantype(modality, 'GOOD');
    end

    if ~isempty(chanind)
        forward(ind).channels = D.chanlabels(chanind);
    else
        error(['No good ' modality ' channels were found.']);
    end
end

if nargin < 3
    channels = [forward(:).channels];
end

try
    fname = D.inv{val}.gainmat;
    G = load(fullfile(D.path, fname)); % Relative path

    label = G.label;
    G     = G.G;
    if numel(label) ~= size(G, 1) || ~all(ismember(channels, label))
        error('Gain matrix has an incorrect number of channels.');
    end
catch
    spm('sFnBanner', mfilename, SVNrev);
    spm('Pointer', 'Watch');

    G     = {};
    label = {};

    for ind = 1:numel(forward)
        %-Create a new lead-field matrix
        %==================================================================

        %-Head Geometry (create tesselation file)
        %------------------------------------------------------------------
        vert = forward(ind).mesh.vert;
        face = forward(ind).mesh.face;

        %-Normals
        %------------------------------------------------------------------
        norm = spm_mesh_normals(struct('faces',face,'vertices',vert),true);

        vol  = forward(ind).vol;

        if ischar(vol)
            vol = ft_read_headmodel(vol);
        end

        modality = forward(ind).modality;

        if isfield(forward, 'siunits') && forward(ind).siunits
            units = D.units(D.indchannel(forward(ind).channels));
            sens  = forward(ind).sensors;
            siunits = isempty(strmatch('unknown', units));
        else
            siunits = false;
            sens = D.inv{val}.datareg(ind).sensors;
        end

        %-Forward computation
        %------------------------------------------------------------------
        [vol, sens] = ft_prepare_vol_sens(vol, sens, 'channel', forward(ind).channels);
        nvert = size(vert, 1);

        spm_progress_bar('Init', nvert, ['Computing ' modality ' leadfields']);
        if nvert > 100, Ibar = floor(linspace(1, nvert,100));
        else Ibar = [1:nvert]; end

        if ~isequal(ft_headmodeltype(vol), 'interpolate')
            Gxyz = zeros(length(forward(ind).channels), 3*nvert);
            for i = 1:nvert
                if siunits
                    Gxyz(:, (3*i - 2):(3*i))  = ft_compute_leadfield(vert(i, :), sens, vol,...
                        'dipoleunit', 'nA*m', 'chanunit', units);
                else
                    Gxyz(:, (3*i - 2):(3*i))  = ft_compute_leadfield(vert(i, :), sens, vol);
                end

                if any(Ibar == i)
                    spm_progress_bar('Set', i);
                end
            end
        else
            if siunits
                Gxyz = ft_compute_leadfield(vert, sens, vol, 'dipoleunit', 'nA*m', 'chanunit', units);
            else
                Gxyz = ft_compute_leadfield(vert, sens, vol);
            end
        end

        spm_progress_bar('Clear');

        spm_progress_bar('Init', nvert, ['Orienting ' modality ' leadfields']);

        G{ind} = zeros(size(Gxyz, 1), size(Gxyz, 2)/3);
        for i = 1:nvert
            G{ind}(:, i) = Gxyz(:, (3*i- 2):(3*i))*norm(i, :)';
            if ismember(i,Ibar)
                spm_progress_bar('Set', i);
            end

        end

        spm_progress_bar('Clear');

        %-Condition the scaling of the lead-field
        %------------------------------------------------------------------
        [Gs, scale] = spm_cond_units(G{ind});

        if siunits && abs(log10(scale))>2
            warning(['Scaling expected to be 1 for SI units, actual scaling ' num2str(scale)]);
            G{ind} = Gs;
        else
            scale = 1;
        end

        label = [label; forward(ind).channels(:)];

        forward(ind).scale = scale;
    end

    if numel(G) > 1
        G = cat(1, G{:});
    else
        G = G{1};
    end

    %-Save
    %----------------------------------------------------------------------
    D.inv{val}.gainmat = ['SPMgainmatrix_' spm_file(D.fname, 'basename') '_' num2str(val) '.mat'];
    save(fullfile(D.path, D.inv{val}.gainmat), 'G', 'label', spm_get_defaults('mat.format'));
    D.inv{val}.forward = forward;
    save(D);

    spm('Pointer', 'Arrow');
end

[sel1, sel2] = spm_match_str(channels, label);

if length(sel2) ~= numel(channels)
    error('Did not find a match for all the requested channels.');
end

L = sparse(G(sel2, :));

%-Retain selected sources if necessary
%--------------------------------------------------------------------------
if nargin > 1 && ~isempty(Is)
    L = L(:,Is);
end

D.inv{val}.forward = forward;

author: spm


text:
% Copyright (C) 2007 Peter Carbonetto. All Rights Reserved.
% This code is published under the Eclipse Public License.
%
% Author: Peter Carbonetto
%         Dept. of Computer Science
%         University of British Columbia
%         May 19, 2007

function A = spzeros (m, n)
  A = sparse([],[],[],m,n,0);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function Z = uplus(Y)
%UPLUS (overloaded)

disp('Most likely you meant to write a + b, but you wrote a +b')
disp('This can easily lead to bugs, as [a +b] is a vector with two elements')
disp('If you really want to use unitary plus, you will have to edit sdpvar/uplus')
disp('and delete this message')
if isa(Y,'blkvar')
    Y = sdpvar(Y);
end

switch Y.typeflag
case {0,5}
	Z = Y;
case {1,3,4}
	error('Relational objects canot be manipulated')
otherwise
	error('Please report internal bug in sdpvar/uplus')
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function options = symbset(varargin)
%SYMBSET  Create or alter the options for Optimization with SymBuilder
%
% options = symbset('param1',value1,'param2',value2,...) creates an
% options structure with the parameters 'param' set to their corresponding
% values in 'value'. Parameters not specified will be set to the SymBuilder
% default.
%
% options = symbset(oldopts,'param1',value1,...) creates a copy of the old
% options 'oldopts' and then fills in (or writes over) the parameters
% specified by 'param' and 'value'.
%
% options = symbset() creates an options structure with all fields set to
% symbset defaults.
%
% symbset() prints a list of all possible fields and their function.

%   Copyright (C) 2013 Jonathan Currie (I2C2)

% Print out possible values of properties.
numberargs = nargin;
if (nargin == 0) && (nargout == 0)
    printfields();
    return
end
%Names and Defaults
Names = {'solver','display','solverOpts','use1stDerivs','use2ndDerivs','preallocate'}';
Defaults = {'auto','iter',[],'yes','yes','yes'}';

%Collect Sizes and lowercase matches
m = size(Names,1);
%Create structure with all names and default values
st = [Names,Defaults]'; options = struct(st{:});

% Check we have char or structure input. If structure, insert arguments
i = 1;
while i <= numberargs
    arg = varargin{i};
    if ischar(arg)
        break;
    end
    if ~isempty(arg)
        if ~isa(arg,'struct')
            error('An argument was supplied that wasn''t a string or struct!');
        end
        for j = 1:m
            if any(strcmp(fieldnames(arg),Names{j,:}))
                val = arg.(Names{j,:});
            else
                val = [];
            end
            if ~isempty(val)
                checkfield(Names{j,:},val);
                options.(Names{j,:}) = val;
            end
        end
    end
    i = i + 1;
end

%Check we have even number of args left
if rem(numberargs-i+1,2) ~= 0
    error('You do not have a value specified for each field!');
end

%Go through each argument pair and assign to correct field
expectval = 0; %first arg is a name
while i <= numberargs
    arg = varargin{i};

    switch(expectval)
        case 0 %field
            if ~ischar(arg)
                error('Expected field name to be a string! - Argument %d',i);
            end
            j = find(strcmp(arg,Names) == 1);
            if isempty(j)  % if no matches
                error('Unrecognised parameter %s',arg);
            elseif(length(j) > 1)
                error('Ambiguous parameter %s',arg);
            end
            expectval = 1; %next arg is a value
        case 1
            checkfield(Names{j,:},lower(arg));
            options.(Names{j,:}) = lower(arg);
            expectval = 0;
    end
    i = i + 1;
end

if expectval %fallen off end somehow
    error('Missing value for %s',arg);
end

%If not using 1st derivatives, turn off 2nd too
if(strcmpi(options.use1stDerivs,'no')), options.use2ndDerivs = 'no'; end

function checkfield(field,value)
%Check a field contains correct data type

if isempty(value)
    return % empty matrix is always valid
end

switch lower(field)
    %String
    case {'use1stderivs','use2ndderivs','preallocate'}
        if(ischar(value) && (strcmpi(value,'yes') || strcmpi(value,'no')))
            valid = true;
        else
            valid = false; errmsg = sprintf('Parameter %s should be a string with ''yes'' or ''no''',field);
        end

    case 'solver'
        if(ischar(value) && checkSolver(value))
            valid = true;
        else
            valid = false; errmsg = sprintf('Parameter %s should be a valid solver',field);
        end
    case 'display'
        if(ischar(value) && any(strcmpi(value,{'off','iter','final'})))
            valid = true;
        else
            valid = false; errmsg = sprintf('Parameter %s should be ''off'', ''iter'' or ''final''',field);
        end
    case 'solveropts'
        if(isstruct(value))
            valid = true;
        else
            valid = false; errmsg = sprintf('Parameter %s should be a structure',field);
        end

    otherwise
        valid = false;
        errmsg = sprintf('Unrecognized parameter name ''%s''.', field);
end

if valid
    return;
else %error
    ME = MException('symbset:FieldError',errmsg);
    throwAsCaller(ME);
end

function printfields()
%Print out fields with defaults
fprintf('        solver: [ Optimization Solver {''auto''} ] \n');
fprintf('       display: [ Solver Display Level ''off'', {''iter''}, ''final'' ]\n');
fprintf('    solverOpts: [ Solver specific options {[]} ] \n');
fprintf('  use1stDerivs: [ Use first derivatives {''yes''} (Turns off use2ndDerivs too if no) ] \n');
fprintf('  use2ndDerivs: [ Use second derivatives {''yes''} ] \n');
fprintf('   preallocate: [ Preallocate return matrices/vectors {''yes''} ] \n');
fprintf('\n');


text:
function out = cov(f, g, varargin)
%COV   Covariance of a CHEBFUN.
%   COV(F) is the same as VAR(F) if F is a scalar-valued CHEBFUN.
%   COV(F) returns the covariance of the array-valued CHEBFUN F.
%   COV(F, G) returns the covariance matrix of the columns of F and G.
%
% See also VAR, MEAN.

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

% Empty case:
if ( isempty(f) )
    out = NaN;
    return
end

% Error checking:
if ( (nargin == 2) && ~all(size(f) == size(g)) )
    error('CHEBFUN:CHEBFUN:cov:size',' CHEBFUN dimensions do not agree.');
end
if ( nargin == 3 )
    error('CHEBFUN:CHEBFUN:cov:nargin', ...
        'CHEBFUN/COV does not support normalization.');
end

% Deal with row CHEBFUN objects:
if ( f(1).isTransposed )
    if ( nargin == 1 )
        out = transpose(cov(transpose(f)));
    else
        out = transpose(cov(transpose(f), transpose(g)));
    end
    return
end

% Conditional on COV(f) and COV(f, g).
if ( nargin == 1 ) % COV(f)

    if ( numColumns(f) == 1 )
        % The covariance of a scalar-valued CHEBFUN is the same as the variance:
        out = var(f);
        return

    else
        % Array-valued CHEBFUN or quasimatrix.

        Y = f - mean(f);
        out = diag(mean(Y.*conj(Y)));
        % Convert Y to a cell array of scalar-valued CHEBFUN objects.
        Y = mat2cell(Y);
        % Loop over each of the columns:
        for j = 1:numel(Y)
            for k = j+1:numel(Y)
                % Compute the scaled inner product of the jth and kth columns:
                out(j,k) = mean(Y{j}.*conj(Y{k}));
                % Use symmetry:
                out(k,j) = conj(out(j,k));
            end
        end

    end

else               % COV(f, g)

    % Convert to cell arrays of scalar-valued CHEBFUN objects:
    Y = cheb2cell(f - mean(f));
    Z = cheb2cell(g - mean(g));
    % Initialise output matrix:
    out = zeros(numel(f));
    % Loop over each of the columns:
    for j = 1:numel(Y)
        for k = 1:numel(Y)
            out(j,k) = mean(Y{j}.*conj(Z{k}));
        end
    end

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% DEPRECATED...
%
% Author: Javier Lopez-Calderon
% Center for Mind and Brain
% University of California, Davis,
% Davis, CA
% 2009

function EEG = checkzerolat(EEG)
nepoch = EEG.trials;
if isempty(EEG.epoch)
      fprintf('\nWARNING: checkzerolat() only works for epoched datasets. This call was ignored.\n')
      return
end
for i=1:nepoch
      latcell  =  cell2mat(EEG.epoch(i).eventlatency);
      [v, indx] = min(abs(latcell));
      EEG.epoch(i).eventlatency = num2cell(latcell - latcell(indx));
end

auxtimes  = EEG.times;
[v, indx]  = min(abs(auxtimes));
EEG.times = auxtimes - auxtimes(indx);
EEG.xmin  = min(EEG.times)/1000;
EEG.xmax  = max(EEG.times)/1000;
EEG.srate = round(EEG.srate);
EEG = eeg_checkset( EEG );

if EEG.times(1)~=auxtimes(1)
      msg = ['\nWarning: zero time-locked stimulus latency values were not found.\n'...
      'Therefore, ERPLAB adjusted latency values at EEG.epoch.eventlatency, EEG.times, EEG.xmin,and EEG.xmax.\n\n'];
      fprintf(msg);
      fprintf('Time range is now [%.3f  %.3f] sec.\n', EEG.min, EEG.max )
else
      fprintf('Zero latencies OK.\n')
end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function a = rutis3_eigen_right ( )

%*****************************************************************************80
%
%% RUTIS3_EIGEN_RIGHT returns the right eigenvectors of the RUTIS3 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, complex A(4,4), the right eigenvector matrix.
%
  i = sqrt ( -1.0 );

  a(1:4,1:4) = [
     1.0,  - 1.0,    1.0,   1.0; ...
     1.0,  - i     - i,    -1.0; ...
     1.0,    i,      i,    -1.0; ...
     11.0,   1.0,  - 1.0,   1.0 ]';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%% simple example of use of renorm_sibling_3d
clear; close all;
x = mandrill;
%% compute roto-translation scattering of an image
options.Q = 2;
Wop = wavelet_factory_3d_pyramid(options, options, options);
Sx = scat(x, Wop);

%% L1 renormalization
op = @(x)(sum(x, 3));
Sx_renorm = renorm_sibling_3d(Sx, op);

%% L1 + smoothing renormalization
op = renorm_factory_L1_smoothing(1);
[Sx_renorm, siblings] = renorm_sibling_3d(Sx, op);

%%
image_scat(Sx, 0, 0);
%%
image_scat(Sx_renorm, 0, 0);
%%
close all;
imagesc(log(image_scat_layer(Sx{3},0,1)+0.001));

%%
image_scat(Sx_renorm, 1, 1);
%%
image_scat(Sx, 1, 1);

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function calibSize = calibrationSize(obj,mask,outlierRate,initValues)
% return the calibration size of the input k-space region
%
% (c) Thomas Kuestner
% ---------------------------------------------------------------------

sx = 2;
sy = 2;
sz = 2;

if(nargin < 3)
    outlierRate = [0 0 0];
end

if(nargin == 4)
    sx = initValues(1);
    sy = initValues(2);
    sz = initValues(3);
end

if(sx >= size(mask,1))
    sx = size(mask,1); % for guaranteed correct cropping
    xflag = true;
else
    xflag = false;
end
if(sy >= size(mask,2))
    sy = size(mask,2);
    yflag = true;
else
    yflag = false;
end
if(sz >= size(mask,3))
    sz = size(mask,3);
    zflag = true;
else
    zflag = false;
end

while(~(xflag && yflag && zflag))

    if(~xflag)
        mask_helper = crop(mask,[sx+1,sy,sz]);
        if(all(mask_helper))
            sx = sx + 1;
        else
            if(outlierRate(1) > 0)
                if(nnz(prod(prod(mask_helper,2),3) == 0) < outlierRate(1) * size(mask_helper,1))
                    sx = sx + 1;
                else
                    xflag = true;
                end
            else
                xflag = true;
            end
        end
    end

    if(~yflag)
        mask_helper = crop(mask,[sx,sy+1,sz]);
        if(all(mask_helper))
            sy = sy + 1;
        else
            if(outlierRate(2) > 0)
                if(nnz(prod(prod(mask_helper,1),3) == 0) < outlierRate(2) * size(mask_helper,2))
                    sy = sy + 1;
                else
                    yflag = true;
                end
            else
                yflag = true;
            end
        end
    end

    if(~zflag)
        mask_helper = crop(mask,[sx,sy,sz+1]);
        if(all(mask_helper))
            sz = sz + 1;
        else
            if(outlierRate(3) > 0)
                if(nnz(prod(prod(mask_helper,1),2) == 0) < outlierRate(3) * size(mask_helper,3))
                    sz = sz + 1;
                else
                    zflag = true;
                end
            else
                zflag = true;
            end
        end
    end

    if(sx == size(mask,1))
        xflag = true;
    end
    if(sy == size(mask,2))
        yflag = true;
    end
    if(sz == size(mask,3))
        zflag = true;
    end
end

if(size(mask,3) > 1)
    calibSize = [sx,sy,sz];
else
    calibSize = [sx,sy];
end
if(any(outlierRate < 0))
    maskSize = size(mask);
    calibSize(outlierRate < 0) = maskSize(outlierRate < 0);
end

end

author: thomaskuestner


text:
function K = covPERard(cov, hyp, x, z, i)

% Stationary periodic covariance function for a stationary covariance function
% k0 such as covMaternard, covPPard, covRQard and covSEard.
% Stationary means that the covariance function k0(x,z) depends on the
% data points x,z only through the squared distance
% dxz = (x-z)'*inv(P)*(x-z) where the P matrix is diagonal with ARD parameters
% ell_1^2,...,ell_D^2, where D is the dimension of the input space.
% The covariance function is parameterized as:
%
% k(x,z) = k0(u(x),u(z)), u(x) = [sin(pi*x/p); cos(pi*x/p)]
%
% where the period p belongs to covPERiso and hyp0 belong to k0:
%
% hyp = [ log(p_1)
%         log(p_2)
%          .
%         log(p_D)
%         hyp0 ]
%
% The first D hyperparameters of k0 are the log lengthscales such that
% hyp0(i) = log(ell(i)) for i=1..D.
% Note that for k0 = covSEard and D = 1, a faster alternative is covPeriodic.
%
% Copyright (c) by Hannes Nickisch, 2013-10-16.
%
% See also COVFUNCTIONS.M.

nocov = false;                   % default case when no cov argument is provided
if nargin==0, cov = {@covSEard}; nocov = true; end                % default case
if isnumeric(cov)       % detect old version where the cov parameter was missing
  % i <- z, z <- x, x <- hyp, hyp <- cov
  if nargin>3, i = z; end
  if nargin>2, z = x; end
  if nargin>1, x = hyp; end
  hyp = cov; cov = {@covSEard}; nocov = true;
end

if nocov && nargin<2 || ~nocov && nargin<3         % report number of parameters
  K = ['(D+',feval(cov{:}),')']; return
end
if nocov && nargin<3 || ~nocov && nargin<4, z = []; end    % make sure, z exists
xeqz = isempty(z); dg = strcmp(z,'diag');                       % determine mode

[n,D] = size(x);
p = exp(hyp(1:D));
lell = hyp(D+(1:D));
hyp0 =  [lell; lell; hyp(2*D+1:end)];

if nocov && nargin<4 || ~nocov && nargin<5
  [x,z] = u(x,z,p,dg);                      % apply the embedding u:IR^D->IR^2*D
  K = feval(cov{:},hyp0,x,z);
else
  if i<=D
    if dg                                                   % compute distance d
      di = zeros([n,1]);
    else
      if xeqz                                             % symmetric matrix Kxx
        di = repmat(reshape(x(:,i),n, 1),[1, n])...
            -repmat(reshape(x(:,i),1, n),[n, 1]);
      else                                               % cross covariances Kxz
        nz = size(z,1);
        di = repmat(reshape(x(:,i),n, 1),[1,nz])...
            -repmat(reshape(z(:,i),1,nz),[n, 1]);
      end
    end
    di = 2*pi*di/p(i); dD2_dlpi = -2*sin(di).*di;   % derivative dD2i/dlog(p(i))
    [x,z] = u(x,z,p,dg);                    % apply the embedding u:IR^D->IR^2*D
    if dg                                            % compute squared distances
      D2 = zeros(n,1);
    else
      if xeqz
        D2 = sq_dist(x(:,[i,i+D])');
      else
        D2 = sq_dist(x(:,[i,i+D])',z(:,[i,i+D])');
      end
    end
    % reconstruct derivative w.r.t. D2i from derivative w.r.t. log(ell(i))
    dK_dD2 = feval(cov{:},hyp0,x,z,i) + feval(cov{:},hyp0,x,z,i+D);
    dK_dD2 = dK_dD2./(-2*D2); dK_dD2(D2<1e-12) = 0;
    K = dK_dD2.*dD2_dlpi;                                     % apply chain rule
  else
    [x,z] = u(x,z,p,dg);                    % apply the embedding u:IR^D->IR^2*D
    if i<=2*D
      K = feval(cov{:},hyp0,x,z,i-D) + feval(cov{:},hyp0,x,z,i);
    else
      K = feval(cov{:},hyp0,x,z,i);
    end
  end
end

function [x,z] = u(x,z,p,dg)                % apply the embedding u:IR^D->IR^2*D
  x = 2*pi*x*diag(1./p); x = [sin(x), cos(x)];
  if numel(z)>0 && ~dg, z = 2*pi*z*diag(1./p); z = [sin(z), cos(z)]; end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [cmds, errs, msgs] = gitSubmodulesClone(repo, preview, options)
cmds = {};

if ~exist('repo','var') || isempty(repo)
    repo = [pwd, '/'];
end
if ~exist('preview','var')
    preview = false;
end
if ~exist('options','var')
    options = 'init';
end

repoFull = filesepStandard_startup(repo,'full');
ii = 1;

submodules = parseGitSubmodulesFile(repoFull);
url = gitGetOrigin(repoFull);
urlroot = fileparts(url);
branch = gitGetBranch(repoFull);

if strcmp(options, 'update')
    try
        rmdir([repoFull, '/submodules'],'s')
    catch
    end
end
if ~ispathvalid([repoFull, '/submodules'])
    mkdir([repoFull, '/submodules'])
end

cmds{ii,1} = sprintf('cd %s', repoFull); ii = ii+1;
for jj = 1:size(submodules,1)
    [~, submodulename] = fileparts(submodules{jj,1});
    if ~ispathvalid([repoFull, 'submodules/', submodulename])
        cmds{ii,1} = sprintf('git clone --branch %s %s %s', branch, [urlroot, '/', submodulename], [repoFull, 'submodules/', submodulename]); ii = ii+1;
    end
end

[errs, msgs] = exeShellCmds(cmds, preview);

author: BUNPC


text:
classdef BT3 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP with bias feature

%------------------------------- Reference --------------------------------
% H. Li, Q. Zhang, and J. Deng, Biased multiobjective optimization and
% decomposition algorithm, IEEE Transactions on Cybernetics, 2017, 47(1):
% 52-66.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            [N,D]  = size(X);
            I1     = 2 : 2 : D;
            I2     = 3 : 2 : D;
            Y      = X - sin(repmat(1:D,N,1)*pi/2/D);
            X(:,1) = abs(X(:,1)).^0.02;
            PopObj(:,1) = X(:,1)         + sum(Y(:,I1).^2+(1-exp(-Y(:,I1).^2/1e-8))/5,2);
            PopObj(:,2) = 1-sqrt(X(:,1)) + sum(Y(:,I2).^2+(1-exp(-Y(:,I2).^2/1e-8))/5,2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^0.5;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function fvc = u3d_pre
% Syntax:
%     fvc = u3d_pre
%
% Description:
%    This function generates the input for the MESH_TO_LATEX function from
%    Alexandre Gramfort from your surface-graphs. The surface graphs 3d-model can be
%    displayed in u3d - format in pdf or if you don't delete the u3d-files
%    which are generated by MESH_TO_LATEX you can use deepview from
%    righthemisphere(http://www.righthemisphere.com/products/client-products/deep-view)
%    to embed the modell in Microsoft-Office products (Word, Excel, PowerPoint).
%
%    The output is a structure generated by the standard MATLAB function
%    surf2patch which can used in MESH_TO_LATEX:
%       fvc.vertices         -->    points
%       fvc.faces            -->    faces
%       fvc.facevertexcdata  -->    face_vertex_data
%
% Example:
%       [X,Y,Z] = peaks(30);
%       surf(X,Y,Z);
%       fvc = u3d_pre;
%       mesh_to_latex('tex/mesh',fvc.vertices,uint32(fvc.faces),fvc.facevertexcdata);
%
% Bugs and suggestions:
%    Please send to Sven Koerner: koerner(underline)sven(add)gmx.de
%
% See also:
%    MESH_TO_LATEX on the File Exchange
%    http://www.righthemisphere.com/products/client-products/deep-view
%
%
% License to use and modify this code is granted freely to all interested, as long as the original author is
% referenced and attributed as such. The original author maintains the right to be solely associated with this work.
%
% Programmed by Sven Koerner: koerner(underline)sven(add)gmx.de
% Date: 2010/04/14

% search for surface-graph
h = findobj('type','surface');

if ~isempty(h)
    % get defined data-points
    X = get(h(end,1), 'XData');
    Y = get(h(end,1), 'YData');
    Z = get(h(end,1), 'ZData');

    % scaled color to unscaled r
    cdata = get(h(end,1), 'CData');
    siz = size(cdata);
    cmap = colormap;
    nColors = size(cmap,1);
    cax = caxis;
    idx = ceil( (double(cdata) - cax(1)) / (cax(2)-cax(1)) * nColors);
    idx(idx<1) = 1;
    idx(idx>nColors) = nColors;
    %handle nans in idx
    nanmask = isnan(idx);
    idx(nanmask)=1; %temporarily replace w/ a valid colormap index
    realcolor = zeros(siz);
    for i = 1:3,
        c = cmap(idx,i);
        c = reshape(c,siz);
        realcolor(:,:,i) = c;
    end

    fvc = surf2patch(X,Y,Z,realcolor, 'triangles' );
else
    fvc.vertices            = [];
    fvc.faces               = [];
    fvc.facevertexcdata     = [];
    disp('No surface graph found.');
end;

author: Sable

repo: mcbench-benchmarks


text:
function ROT = createRotationVector3d(A,B)
%CREATEROTATIONVECTOR3D Calculates the rotation between two vectors.
%
%   ROT = createRotationVector3d(A, B) returns the 4x4 rotation matrix ROT
%   to transform vector A in the same direction as vector B.
%
%   Example
%     A=[ .1  .2  .3];
%     B=-1+2.*rand(1,3);
%     ROT = createRotationVector3d(A,B);
%     C = transformVector3d(A,ROT);
%     figure('color','w'); hold on; view(3)
%     O=[0 0 0];
%     drawVector3d(O, A,'r');
%     drawVector3d(O, B,'g');
%     drawVector3d(O, C,'r');
%
%   See also
%   transformPoint3d, createRotationOx, createRotationOy, createRotationOz
%
%   Source
%     https://math.stackexchange.com/a/897677
%
% ---------
% Author: oqilipo
% Created: 2017-08-07
% Copyright 2017

if isParallel3d(A,B)
    if A*B'>0
        ROT = eye(4);
    else
        ROT = -1*eye(4); ROT(end)=1;
    end
else
    a=normalizeVector3d(A);
    b=normalizeVector3d(B);
    a=reshape(a,3,1);
    b=reshape(b,3,1);

    v = cross(a,b);
    ssc = [0 -v(3) v(2); v(3) 0 -v(1); -v(2) v(1) 0];
    ROT = eye(3) + ssc + ssc^2*(1-dot(a,b))/(norm(v))^2;

    ROT = [ROT, [0;0;0]; 0 0 0 1];
end

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


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

%%
clear; close all; clc;

%% Syntax
% |[a]=DCM2euler(Q);|

%% Description
% This function is the inverse of the euler2DCM function. The Euler angles
% |a| are derived based on the input rotatin tensor |Q|.

%% Examples
%

%%
% Plot settings
fontSize=25;

%% Retrieving the Euler andles from a rotation tensor

%%
% Get example patch data
[F,V]=parasaurolophus;

%%
% Defining sets of true Euler angles for X, Y and Z axis rotation
a_true=[0.25*pi 0.25*pi -0.25*pi]

%%
% Use |euler2DCM| function to define the rotation tensor
[Q]=euler2DCM(a_true);

%%
% Use |DCM2euler| to retrieve the Euler angles
a_fit=DCM2euler(Q)

%% Handling symbolic expressions

try
    syms a b c

    a_true=[a b c]
    Q=euler2DCM(a_true);

    a_fit=DCM2euler(Q)

catch
    warning('Symbolic toolbox likely missing')
end

%%
%
% <<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 [gKern, gVarmeans, gVarcovars, gInd] = rbfardjitVardistPsi2Gradient(rbfardKern, vardist, Z, covGrad, learnInducing)

% RBFARDJITVARDISTPSI2GRADIENT description.

% VARGPLVM

if nargin < 5
    learnInducing = 1;
end

[gKern, gVarmeans, gVarcovars, gInd] = rbfard2VardistPsi2Gradient(rbfardKern, vardist, Z, covGrad, learnInducing);

author: SheffieldML

repo: vargplvm


text:
function sorted_model_files = sort_nnet_by_dev(model_files)

for i=1:length(model_files)
    idx = regexp(model_files{i}, '\.CV');
    CV_score(i) = str2num( model_files{i}(idx(end)+3:end-4) );
end

[tmp, idx] = sort(CV_score);
sorted_model_files = model_files(idx);

end

author: singaxiong

repo: SignalGraph


text:
function coordl1bregdemo
%COORDL1BREGDEMO An example of using COORDL1BREG.

% Construct A as a rand matrix.
N = 512*2;
M = N/2;
A = randn(M,N);

% Construct u_exact as a sparse vector.
p = floor(0.05*N);
u_exact = zeros(N,1);
% u_exact has spikes at uniformly random locations with amplitudes
% distributed uniformly in [0,1]
a = randperm(N);
% u_exact(a(1:p)) = a(p+1:2*p)*0.08+5; % 0.08 + 5
u_exact(a(1:p)) = rand(p,1)*N;  %randn is fast.

% Construct f = A*u_exact.
f = A*u_exact;

% Precompute B = A'*A. This step is not necessary.
B = A'*A;

% Initialize lambda.
lambda = 10;

% Call the main function  COORDL1BREG.
[u,Energy] = coordl1breg(A,f,lambda,'B',B,'PlotFun',@myplotfun);

% Plot the Energy function.
figure(2);
semilogy(Energy,'.-');
xlabel('Iteration');
ylabel('Energy');

% COORDL1BREG will call this function every iteration to plot intermediate solutions.
% Show the comparison of the output u and u_exact.
    function myplotfun(u)

        figure(1);
        x = 1:length(u);
        plot(x,u,'.r',x,u_exact,'o');
        xlim([1,length(u)]);
    end

end


text:
function [state, xf] = ft_preproc_online_filter_apply(state, x)

% FT_PREPROC_ONLINE_FILTER_APPLY passes a signal through the online filter and
% returns the updated filter model (delay states) and the filtered signal.
%
% Use as
%   [state, dat] = ft_preproc_online_filter_apply(state, dat)
% where
%   dat   = Nchan x Ntime
%   state = filter state, see FT_PREPROC_ONLINE_FILTER_INIT
%
% See also PREPROC

% Copyright (C) 2010, Stefan Klanke
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

[dimX, numX] = size(x);

if dimX > numX
  % explicit algorithm - faster for many channels, 1 sample (=fMRI)
  xf = zeros(size(x));

  for k=1:numX;
    z_old = state.z;
    z0 = x(:,k) - z_old*state.A2;
    xf(:,k) = z_old * state.B2 + z0*state.B1;
    state.z(:,2:end) = z_old(:,1:end-1);
    state.z(:,1) = z0;
  end
else
  % use built-in MATLAB stuff - faster for many samples, few channels
  [xf, z] = filter(state.B, state.A, x, state.z',2);
  state.z = z';
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [RR_e_resamp,QRS_A_resamp] = clean_and_resample_intervals(onsets,QRS_A1,Fs)

% Use same pre-processing used for ECG in the paper ''
%
%

IBI = diff(onsets)./Fs;
t_IBI = onsets(2:end)./Fs;

RR_e0 = IBI;     % no ectopic rejection
QRS_A0 = QRS_A1;

% reject ectopic beats
% remove invalid R-R intervals if RR<0.3 or RR>2.0 & 20% outside of 41 points moving mean exclude

RR_mean_41 = meanfilt1(IBI,41);
RR = [];    % RR interval in seconds
RR_time=[]; % time of RR interval in seconds
v_RR = ones(1,length(IBI))*-1;
n_RR = 0;
for j=1:length(IBI)
    if abs(IBI(j) - RR_mean_41(j)) < 0.2*RR_mean_41(j) &&  IBI(j)>= 0.3 && IBI(j) <= 2.0
        n_RR = n_RR+1;
        RR(n_RR) = IBI(j);
        RR_time(n_RR) = onsets(j)/Fs;
        v_RR(j)=1;
    end
end
valid_RR=find(v_RR==1);
% remove invalid QRS peaks if 50% outside of 41 points moving mean
QRS_A_mean_41=meanfilt1(QRS_A1,41);
n_QRS=0;
QRS=[];
QRS_time=[];
v_QRS=ones(1,length(QRS_A1))*-1;
for k=1:length(QRS_A1)
    if abs(QRS_A1(k)-QRS_A_mean_41(k))<0.5*QRS_A_mean_41(k)
        n_QRS=n_QRS+1;
        QRS(n_QRS) = QRS_A1(k);
        QRS_time(n_QRS) = tIBI(j);
        v_QRS(k)=1;
    end
end
valid_QRS=find(v_QRS==1);
valid_RRQRS=intersect(valid_RR,valid_QRS);
valid_RR=valid_RRQRS;

onsets = onsets(valid_RR);
QRS_A1=QRS_A1(valid_RR);

RR_e1 = IBI(valid_RR);

% resample to 4Hz
onsets = onsets./Fs;
t2=round(onsets(end));
Fs_resamp=4;
t1=(1/Fs_resamp):(1/Fs_resamp):t2; % 4Hz resample

% linear interpolation
QRS_A_resamp = interp1(onsets,QRS_A1,t1);%,'spline');
RR_e_resamp = interp1(onsets,RR_e1,t1);%,'spline');
QRS_A_resamp(find(isnan(QRS_A_resamp)))=nanmean(QRS_A_resamp);
RR_e_resamp(find(isnan(RR_e_resamp)))=nanmean(RR_e_resamp);

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function V2=tform(T,V)

%%
nDim_T=size(T,2);
nDim_V=size(V,2);

if nDim_T==(nDim_V+1)
        VV=V;
        VV(:,end+1)=1;
        V2=(T*VV')';
        V2=V2(:,[1 2 3]);
elseif nDim_T==nDim_V
        V2=(T*V')';
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/>.


text:
function [out] = interception_5(p1,p2,In)
%interception_5

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Interception excess after a combined absolute amount and fraction are intercepted
% Constraints:  f >= 0
% @(Inputs):    p1   - fraction that is not throughfall [-]
%               p2   - constnat interception and evaporation [mm/d]
%               In   - incoming flux [mm/d]

out = max(p1.*In-p2,0);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function poly_coef_count_test ( )

%*****************************************************************************80
%
%% POLY_COEF_COUNT_TEST tests POLY_COEF_COUNT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLY_COEF_COUNT_TEST\n' );
  fprintf ( 1, '  POLY_COEF_COUNT counts the number of coefficients\n' );
  fprintf ( 1, '  in a polynomial of degree DEGREE and dimension DIM\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ' Dimension    Degree     Count\n' );

  for dim = 1 : 3 : 10
    fprintf ( 1, '\n' );
    for degree = 0 : 5
      fprintf ( 1, '  %8d  %8d  %8d\n', ...
        dim, degree, poly_coef_count ( dim, degree ) );
    end
  end

  return
end


text:
%  Figure 7.44      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
% script to generate Fig. 7.44
clf;
f=[-10 1 0;-16 0 1;0 0 0];
g=[0 0 10]';
h=[1 0 0];
j=0;
a=[-f' -h'*h;0*f f];
b=[0*g;g];
c=[g' 0*h];
d=0;
rlocus(a,b,c,d);
axis([-12, 12, -6, 6]);
grid;
title('Fig.7.44  Symmetric root locus for a DC Servo system')


text:
function [prob_lut, prob_map] = getForegroundBackgroundProbs(frame, obj_rect, num_bins, bin_mapping)
%GETFOREGROUNDBACKGROUNDPROBS Computes the probability lookup table for the
%object vs surrounding region model.
% Parameters:
%   frame       Input (color) image cropped to contain only the surrounding
%               region
%   obj_rect    Rectangular object region
%   num_bins    Number of bins per channel (scalar)
%   bin_mapping Maps intensity values to num_bins bins

[rows, cols, layers] = size(frame);
%%size of object
obj_row = round(obj_rect(2));
obj_col = round(obj_rect(1));
obj_width = round(obj_rect(3));
obj_height = round(obj_rect(4));

if obj_row + obj_height > rows, obj_height = rows - obj_row; end
if obj_col + obj_width > cols, obj_width = cols - obj_col; end

if nargout > 1
  prob_map = zeros(rows, cols);
end

if layers == 3
  % Color image
  obj_hist = zeros(num_bins, num_bins, num_bins, 'double');
  surr_hist = zeros(num_bins, num_bins, num_bins, 'double');

  % Histogram over full image
  [x,y] = meshgrid(1:cols, 1:rows);
  idx_map = sub2ind([rows, cols], y(:), x(:));
  idx_1 = sub2ind([rows, cols, layers], y(:), x(:), ones(numel(x), 1));
  idx_2 = sub2ind([rows, cols, layers], y(:), x(:), 2.*ones(numel(x), 1));
  idx_3 = sub2ind([rows, cols, layers], y(:), x(:), 3.*ones(numel(x), 1));

  bin_1 = bin_mapping(frame(idx_1)+1);
  bin_2 = bin_mapping(frame(idx_2)+1);
  bin_3 = bin_mapping(frame(idx_3)+1);

  idx_hist_full = sub2ind(size(surr_hist), bin_1, bin_2, bin_3);
  bins = unique(idx_hist_full);
  for b = bins
    surr_hist(b) = nnz(idx_hist_full == b);
  end

  % Histogram over object region
  [x,y] = meshgrid(max(1,obj_col):(obj_col+obj_width), max(1,obj_row):(obj_row+obj_height));
  idx_1 = sub2ind([rows, cols, layers], y(:), x(:), ones(numel(x), 1));
  idx_2 = sub2ind([rows, cols, layers], y(:), x(:), 2.*ones(numel(x), 1));
  idx_3 = sub2ind([rows, cols, layers], y(:), x(:), 3.*ones(numel(x), 1));

  bin_1o = bin_mapping(frame(idx_1)+1);
  bin_2o = bin_mapping(frame(idx_2)+1);
  bin_3o = bin_mapping(frame(idx_3)+1);

  idx_hist_obj = sub2ind(size(obj_hist), bin_1o, bin_2o, bin_3o);
  bins = unique(idx_hist_obj);
  for b = bins
    obj_hist(b) = nnz(idx_hist_obj == b);
  end

  prob_lut = (obj_hist + 1) ./ (surr_hist + 2);
  if nargout > 1
    prob_map(idx_map) = prob_lut(idx_hist_full);
  end
elseif layers == 10
  obj_hist = zeros(1,1,layers, 'double');
  surr_hist = zeros(1,1,layers, 'double');

  for r = 1:rows
    for c = 1:cols
      surr_hist = surr_hist + (frame(r,c,:)+1)/2;
      if r >= obj_row && r < obj_row+obj_height && c >= obj_col && c < obj_col+obj_width
        obj_hist = obj_hist + (frame(r,c,:)+1)/2;
      end
    end
  end
  prob_lut = (obj_hist + 1) ./ (surr_hist + 2);

elseif layers == 1
  % Color image
  obj_hist = zeros(num_bins, 1, 'double');
  surr_hist = zeros(num_bins, 1, 'double');

  % Histogram over full image
  [x,y] = meshgrid(1:cols, 1:rows);
  idx_map = sub2ind([rows, cols], y(:), x(:));

  bin = bin_mapping(frame(idx_map)+1);

  idx_hist_full = sub2ind(size(surr_hist), bin);
  bins = unique(idx_hist_full);
  for b = bins
    surr_hist(b) = nnz(idx_hist_full == b);
  end

  % Histogram over object region
  [x,y] = meshgrid(max(1,obj_col):(obj_col+obj_width), max(1,obj_row):(obj_row+obj_height));
  idx = sub2ind([rows, cols], y(:), x(:));

  bin_1o = bin_mapping(frame(idx)+1);

  idx_hist_obj = sub2ind(size(obj_hist), bin_1o);
  bins = unique(idx_hist_obj);
  for b = bins
    obj_hist(b) = nnz(idx_hist_obj == b);
  end

  prob_lut = (obj_hist + 1) ./ (surr_hist + 2);
  if nargout > 1
    prob_map(idx_map) = prob_lut(idx_hist_full);
  end
else
  error('Not supported\n');
end

end

author: flyers

repo: drone-tracking


text:
function [x, infos] = deep_bidirectional_nmf(X, rank_layers, in_options)
% Deep Bidir-Semi-NMF.
%
% The problem of interest is defined as
%
%           min || X - Z_1 * Z_2 * ... * Z_n * H_n ||_F^2,
%           where
%           {H_n} >= 0.
%
% Given a matrix X, factor matrices {Z_1, Z_2, ..., Z_n, H_n} are calculated.
%
%
% Inputs:
%       X           : (m x n) matrix to factorize
%       rank_layers : ranks in each layer
%       in_options
%
%
% Output:
%       x           : non-negative matrix solution, i.e., x.Z (cell), x.H (cell)
%       infos       : log information
%           epoch   : iteration nuber
%           cost    : objective function value
%           optgap  : optimality gap
%           time    : elapsed time
%           grad_calc_count : number of sampled data elements (gradient calculations)
%
% References
%           G. Trigeorgis, K. Bousmalis, S. Zafeiriou and B. Schuller,
%           "A deep semi-NMF model for learning hidden representations",
%           ICML2014, 2014.
%
%           G. Trigeorgis, K. Bousmalis, S. Zafeiriou and B. Schuller
%           "A deep matrix factorization method for learning attribute representations,"
%           IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), vol.39, no.3, pp.417-429, 2017
%
%
% This file is part of NMFLibrary
%
% Originally created by G.Trigeorgis.
%
% Change log:
%
%       Jul. 26, 2018 (Hiroyuki Kasai): Modified code structures.
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

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

    % set the number of rank_layers
    num_of_layers = numel(rank_layers);

    % set local options
    local_options.bUpdateZ  = true;
    local_options.bUpdateH  = true;
    local_options.bUpdateLastH = true;
    local_options.deepZ     = true;
    local_options.deepH     = true;
    local_options.updateH_alg = 'mu';  % 'mu' or 'apg'
    local_options.eval_clustering_acc = 0;
    local_options.classnum      = 0;

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

    % initialize
    method_name = 'Deep-Bidir-SemiNMF';
    epoch = 0;
    grad_calc_count = 0;

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

    if ~isempty(options.gnd) && options.classnum > 1
        options.eval_clustering_acc = 1;
    end

    % initialize Z and H
    Z = cell(1, num_of_layers);
    H = cell(1, num_of_layers);
    P = cell(1, num_of_layers);
    Q = cell(1, num_of_layers);
    %A = cell(1, num_of_layers+1);
    %B = cell(1, num_of_layers);
    if ~isfield(options, 'x_init')

        if options.deepH
            for i_layer = 1:num_of_layers

                if options.verbose > 1
                    fprintf('### Initializing by %s for layer %d ... ', method_name, i_layer);
                end

                if i_layer == 1
                    % For the first layer we go linear from X to Z*H, so we use id
                    V = X;
                else
                    V = H{i_layer-1};
                end

                % For the later rank_layers we use nonlinearities as we go from
                % g(H_{k-1}) to Z*H_k
                semi_nmf_options.max_iter  = options.max_epoch;
                semi_nmf_options.bUpdateH  = options.bUpdateH;
                semi_nmf_options.bUpdateZ  = options.bUpdateZ;
                semi_nmf_options.verbose   = 0;

                [semi_nmf_x, ~] = semi_mu_nmf(V, rank_layers(i_layer), semi_nmf_options);
                Z{i_layer} = semi_nmf_x.W;
                H{i_layer} = semi_nmf_x.H;

                %fprintf('V: %5.5f, Zi: %5.5f, Hi: %5.5f\n', norm(V), norm(Z{i_layer}), norm(H{i_layer}));
                if options.verbose > 1
                    fprintf('done\n');
                end
            end
        end

        if options.deepZ
            XT = X';
            for i_layer = 1:num_of_layers

                if options.verbose > 1
                    fprintf('### Initializing by %s for layer %d ... ', method_name, i_layer);
                end

                if i_layer == 1
                    % For the first layer we go linear from X to Z*H, so we use id
                    V = XT;
                else
                    V = Q{i_layer-1};
                end

                % For the later rank_layers we use nonlinearities as we go from
                % g(H_{k-1}) to Z*H_k
                semi_nmf_options.max_iter  = options.max_epoch;
                semi_nmf_options.bUpdateH  = options.bUpdateH;
                semi_nmf_options.bUpdateZ  = options.bUpdateZ;
                semi_nmf_options.verbose   = 0;

                [semi_nmf_x, ~] = semi_mu_nmf(V, rank_layers(i_layer), semi_nmf_options);
                P{i_layer} = semi_nmf_x.W;
                Q{i_layer} = semi_nmf_x.H;

                %fprintf('V: %5.5f, Zi: %5.5f, Hi: %5.5f\n', norm(V), norm(Z{i_layer}), norm(H{i_layer}));
                if options.verbose > 1
                    fprintf('done\n');
                end
            end
        end

    else
        Z = options.Z;
        H = options.H;
    end

    Hm = H{num_of_layers};
    Qm = Q{num_of_layers};

    % select disp_freq
    disp_freq = set_disp_frequency(options);

    % store initial info
    clear infos;
    Z_rec = reconst_Z(Z, num_of_layers);

   % Concatinated_H = B{1};
    %Concatinated_Q = Z{1}';

    if options.deepH
        B{num_of_layers} = Hm;
        for i_layer = num_of_layers-1:-1:1
            B{i_layer} = Z{i_layer+1} * B{i_layer+1};
        end
        H_rec = B{1};
        Z_rec = Z{1};
        Concatinated_Q = Z{1}';
    elseif options.deepZ
        B{num_of_layers} = Qm;
        for i_layer = num_of_layers-1:-1:1
            B{i_layer} = P{i_layer+1} * B{i_layer+1};
        end
        H_rec = P{1}';
        Z_rec = B{1}';
    elseif options.deepH && options.deepZ

    end

    %[infos, f_val, optgap] = store_nmf_info(X, Z_rec, Hm, [], options, [], epoch, grad_calc_count, 0);
    [infos, f_val, optgap] = store_nmf_info(X, Z_rec, H_rec, [], options, [], epoch, grad_calc_count, 0);

    % evaluate clustering accuracy
    if ~isempty(options.gnd) && options.classnum > 1
        [infos] = store_clustering_accuracy(Hm, options.gnd, options.classnum, infos, options.eval_clustering_num, 0);
    end

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

    % set start time
    start_time = tic();

    % main loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        % update Z and deep H
        if options.deepH
            Z{1} = Concatinated_Q';
            [Z, Hm, Concatinated_H] = calc_deep_matrices(X, Z, Hm, num_of_layers, options, 1);
        else
            Concatinated_H = P{1}';
        end

%         [infos, f_val, optgap] = store_nmf_info(X, Z{1}, Concatinated_H, [], options, infos, epoch, grad_calc_count, 0);
%         % display infos
%         if options.verbose > 1
%             if ~mod(epoch, disp_freq)
%                 fprintf('Deep-SemiNMF (B): Epoch = %04d, cost = %.16e, optgap = %.4e\n', epoch, f_val, optgap);
%             end
%         end
%
        if options.deepZ
            % update H and deep Z
            P{1} = Concatinated_H';
            [P, Qm, Concatinated_Q] = calc_deep_matrices(XT, P, Qm, num_of_layers, options, 0);
        else
            Concatinated_Q = Z{1}';
        end

        %fprintf('%e\n', norm(XT-P{1}*Concatinated_Q, 'fro')^2 / 2 );

        % measure elapsed time
        elapsed_time = toc(start_time);

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

        % update epoch
        epoch = epoch + 1;

        % store info
        %Z_rec = reconst_Z(Z, num_of_layers);
        %[infos, f_val, optgap] = store_nmf_info(X, Z_rec, Hm, [], options, infos, epoch, grad_calc_count, elapsed_time);
        [infos, f_val, optgap] = store_nmf_info(X, Concatinated_Q', P{1}', [], options, infos, epoch, grad_calc_count, elapsed_time);

%         % display infos
%         if options.verbose > 1
%             if ~mod(epoch, disp_freq)
%                 fprintf('Deep-SemiNMF: Epoch = %04d, cost = %.16e, optgap = %.4e\n', epoch, f_val, optgap);
%             end
%         end

        % evaluate clustering accuracy
        if options.eval_clustering_acc
            [infos] = store_clustering_accuracy(Hm, options.gnd, options.classnum, infos, options.eval_clustering_num, epoch);
        end

        % display infos
        if options.verbose > 1
            if ~mod(epoch, disp_freq)

                if ~options.eval_clustering_acc
                    fprintf('%s: Epoch = %04d, cost = %.16e, optgap = %.4e\n', method_name, epoch, f_val, optgap);
                else
                    fprintf('%s: Epoch = %04d, cost = %.16e, optgap = %.4e, acc = %.4f, nmi = %.4f, purity = %.4f, f = %.4f\n', ...
                        method_name, epoch, f_val, optgap, infos.clustering_acc(end).acc, infos.clustering_acc(end).nmi, infos.clustering_acc(end).purity, infos.clustering_acc(end).f_val);
                end

            end

        elseif options.verbose == 1
            textwaitbar(epoch, options.max_epoch, '  progress');
        end

    end

    H{num_of_layers} = Hm;
    x.Z = Z;
    x.H = H;

end

function [Z, Hm, Concatinated_H] = calc_deep_matrices(X, Z, Hm, num_of_layers, options, Hm_nonnegative)

    A = cell(1, num_of_layers+1);
    B = cell(1, num_of_layers);

    m = size(Z, 1);

    % B{1} = Z{2} * Z{3} * ... * Z{num_of_layers} * Hm
    % B{2} = Z{3} * ... * Z{num_of_layers} * Hm
    % B{3} = ....
    % ....
    % B{num_of_layers} = Hm
    B{num_of_layers} = Hm;
    for i_layer = num_of_layers-1:-1:1
        B{i_layer} = Z{i_layer+1} * B{i_layer+1};
    end

    %% update Z
    % where Z = Z{1} = X * (B{1})^{-1} due to Z{1} B{1} = X.
    Z{1} = X  * pinv(B{1});

    A{1} = eye(m);
    % update Z{2}, ... , Z{num_of_layers}
    for i = 2 : num_of_layers

        % A{1} =
        % A{2} = Z{1}
        % A{3} = Z{1} * Z{2}
        % A{4} = Z{1} * Z{2} * Z{3}
        % ....
        % A{num_of_layers} = Z{1} * Z{2} * ... * Z{num_of_layers-1}
        if i == 2
            A{i} = Z{i-1};
        else
            A{i} = A{i-1} * Z{i-1};
        end

        Z{i} = pinv(A{i}) * X * pinv(B{i});

    end

    % update Hm
    A{num_of_layers+1} = A{num_of_layers} * Z{num_of_layers};
    if Hm_nonnegative
        if strcmp(options.updateH_alg, 'mu')
            P = A{num_of_layers+1}' * X;
            Pp = (abs(P)+P)./2;
            Pn = (abs(P)-P)./2;

            Q = A{num_of_layers+1}' * A{num_of_layers+1};

            Qp = (abs(Q)+Q)./2;
            Qn = (abs(Q)-Q)./2;

            Hm = Hm .* sqrt((Pp + Qn * Hm) ./ max(Pn + Qp * Hm, 1e-16));
        else
            % min_H 1/2 | X - A{num_of_layers+1}*H |^2_F
            % --> min_H 1/2 | X' - H' * A{num_of_layers+1}' |^2_F
            % ----> min_A 1/2 | X - A * B' |^2_F in "nesterov_mnls(X, B, A,..)"
            [tmpH, ~, ~] = nesterov_mnls(X', A{num_of_layers+1}, Hm', 1, options.apg_maxiter, 'basic');
            Hm = tmpH';
        end
    else
        Hm = pinv(A{num_of_layers+1}) * X;
    end

    %
    B{num_of_layers} = Hm;
    for i_layer = num_of_layers-1:-1:1
        B{i_layer} = Z{i_layer+1} * B{i_layer+1};
    end
    Concatinated_H = B{1};

%     % B{1} = Z{2} * Z{3} * ... * Z{num_of_layers} * H{num_of_layers}
%     % B{2} = Z{3} * ... * Z{num_of_layers} * H{num_of_layers}
%     % B{3} = ....
%     % ....
%     % B{num_of_layers} = H{num_of_layers}
%     B{num_of_layers} = H{num_of_layers};
%     for i_layer = num_of_layers-1:-1:1
%         B{i_layer} = Z{i_layer+1} * B{i_layer+1};
%     end
%
%     %% update Z
%     % where Z = Z{1} = X * (B{1})^{-1} due to Z{1} B{1} = X.
%     Z{1} = X  * pinv(B{1});
%
%     A{1} = eye(m);
%     % update Z{2}, ... , Z{num_of_layers}
%     for i = 2 : num_of_layers
%
%         % A{1} =
%         % A{2} = Z{1}
%         % A{3} = Z{1} * Z{2}
%         % A{4} = Z{1} * Z{2} * Z{3}
%         % ....
%         % A{num_of_layers} = Z{1} * Z{2} * ... * Z{num_of_layers-1}
%         if i == 2
%             A{i} = Z{i-1};
%         else
%             A{i} = A{i-1} * Z{i-1};
%         end
%
%         Z{i} = pinv(A{i}) * X * pinv(B{i});
%
%     end
%
%     % update H
%     A{num_of_layers+1} = A{num_of_layers} * Z{num_of_layers};
%     if strcmp(options.updateH_alg, 'mu')
%         P = A{num_of_layers+1}' * X;
%         Pp = (abs(P)+P)./2;
%         Pn = (abs(P)-P)./2;
%
%         Q = A{num_of_layers+1}' * A{num_of_layers+1};
%
%         Qp = (abs(Q)+Q)./2;
%         Qn = (abs(Q)-Q)./2;
%
%
%         H{num_of_layers} = H{num_of_layers} .* sqrt((Pp + Qn * H{num_of_layers}) ./ max(Pn + Qp * H{num_of_layers}, 1e-10));
%     else
%         % min_H 1/2 | X - A{num_of_layers+1}*H |^2_F
%         % --> min_H 1/2 | X' - H' * A{num_of_layers+1}' |^2_F
%         % ----> min_A 1/2 | X - A * B' |^2_F in "nesterov_mnls(X, B, A,..)"
%         [tmpH, ~, ~] = nesterov_mnls(X', A{num_of_layers+1}, H{num_of_layers}', 1, 100, 'basic');
%         H{num_of_layers} = tmpH';
%     end

end

% calculate Z = Z_1 * Z_2 * ... * Z_n
function Z_rec = reconst_Z(Z, num_of_layers)

    Z_rec = Z{num_of_layers};

    for k = num_of_layers-1 : -1 : 1
        Z_rec =  Z{k} * Z_rec;
    end

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function a = asinh(a)
%ASINH        Gradient inverse hyperbolic sine asinh(a)
%

% written  10/16/98     S.M. Rump
% modified 10/14/00     S.M. Rump  use Tony's trick
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/08/08     S.M. Rump  improved sparse multiplication: not using intval data type
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  N = getappdata(0,'INTLAB_GRADIENT_NUMVAR');
  % use full(a.x(:)): cures Matlab V6.0 bug
  % a=7; i=[1 1]; x=a(i), b=sparse(a); y=b(i)  yields row vector x but column vector y
  % ax is full anyway
  ax = 1 ./ sqrt( 1 + sqr(full(a.x(:))) );
  a.x = asinh(a.x);
  if issparse(a.dx)
    sizeax = size(a.dx,1);
    [ia,ja,sa] = find(a.dx);
    if isa(a.x,'intval')
      adx = times(ax(ia),sa,0);
      if adx.complex
        a.dx = intval( sparse(ia,ja,adx.mid,sizeax,N) , sparse(ia,ja,adx.rad,sizeax,N) , 'midrad' );
      else
        a.dx = intval( sparse(ia,ja,adx.inf,sizeax,N) , sparse(ia,ja,adx.sup,sizeax,N) , 'infsup' );
      end
    else
      a.dx = sparse(ia,ja,ax(ia).*sa,sizeax,N);
    end
  else
    a.dx = a.dx .* ax(:,ones(1,N));
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function varargout=loadh5(filename, varargin)
%
%    [data, meta] = loadh5(filename)
%    [data, meta] = loadh5(root_id)
%    [data, meta] = loadh5(filename, rootpath)
%    [data, meta] = loadh5(filename, rootpath, options)
%    [data, meta] = loadh5(filename, options)
%    [data, meta] = loadh5(filename, 'Param1',value1, 'Param2',value2,...)
%
%    Load data in an HDF5 file to a MATLAB structure.
%
%    author: Qianqian Fang (q.fang <at> neu.edu)
%
%    input
%        filename
%            Name of the file to load data from
%        root_id: an HDF5 handle (of type 'H5ML.id' in MATLAB)
%        rootpath : (optional)
%            Root path to read part of the HDF5 file to load
%        options: (optional) a struct or Param/value pairs for user specified options
%            Order: 'creation' - creation order (default), or 'alphabet' - alphabetic
%            Regroup: [0|1]: if 1, call regrouph5() to combine indexed
%                  groups into a cell array
%            PackHex: [1|0]: convert invalid characters in the group/dataset
%                  names to 0x[hex code] by calling encodevarname.m;
%                  if set to 0, call getvarname
%            ComplexFormat: {'realKey','imagKey'}: use 'realKey' and 'imagKey'
%                  as possible keywords for the real and the imaginary part
%                  of a complex array, respectively (sparse arrays not supported);
%                  a common list of keypairs is used even without this option
%
%    output
%        data: a structure (array) or cell (array)
%        meta: optional output to store the attributes stored in the file
%
%    example:
%        a={rand(2), struct('va',1,'vb','string'), 1+2i};
%        saveh5(a,'test.h5');
%        a2=loadh5('test.h5')
%        a3=loadh5('test.h5','regroup',1)
%        isequaln(a,a3.a)
%        a4=loadh5('test.h5','/a1')
%
%    This function was adapted from h5load.m by Pauli Virtanen <pav at iki.fi>
%    This file is part of EasyH5 Toolbox: https://github.com/fangq/easyh5
%
%    License: GPLv3 or 3-clause BSD license, see https://github.com/fangq/easyh5 for details
%

path = '';
if(bitand(length(varargin),1)==0)
    opt=varargin2struct(varargin{:});
elseif(length(varargin)>=3)
    path=varargin{1};
    opt=varargin2struct(varargin{2:end});
elseif(length(varargin)==1)
    path=varargin{1};
end

if(isa(filename,'H5ML.id'))
    loc=filename;
else
    if(exist('h5read','file'))
        loc = H5F.open(filename);
    else
        error('HDF5 is not supported');
    end
end

opt.rootpath=path;

if(~(isfield(opt,'complexformat') && iscellstr(opt.complexformat) && numel(opt.complexformat)==2))
    opt.complexformat={'Real','Imag'};
end

opt.releaseid=0;
vers=ver('MATLAB');
if(~isempty(vers))
    opt.releaseid=datenum(vers(1).Date);
end

if((isfield(opt,'order') && strcmpi(opt.order,'alphabet'))  || opt.releaseid<datenum('1-Jan-2015') )
    opt.order='H5_INDEX_NAME';
else
    opt.order='H5_INDEX_CRT_ORDER';
end

try
  if(nargin>1 && ~isempty(path))
      try
          rootgid=H5G.open(loc,path);
          [varargout{1:nargout}]=load_one(rootgid, opt);
          H5G.close(rootgid);
      catch
          [gname,dname]=fileparts(path);
          rootgid=H5G.open(loc,gname);
          [status, res]=group_iterate(rootgid,dname,struct('data',struct,'meta',struct,'opt',opt));
          if(nargout>0)
              varargout{1}=res.data;
          elseif(nargout>1)
              varargout{2}=res.meta;
          end
          H5G.close(rootgid);
      end
  else
      [varargout{1:nargout}]=load_one(loc, opt);
  end
  H5F.close(loc);
catch ME
  H5F.close(loc);
  rethrow(ME);
end

if(jsonopt('Regroup',0,opt))
    if(nargout>=1)
        varargout{1}=regrouph5(varargout{1});
    elseif(nargout>=2)
        varargout{2}=regrouph5(varargout{2});
    end
end

if(isfield(opt,'jdata') && opt.jdata && nargout>=1)
    varargout{1}=jdatadecode(varargout{1},'Base64',0,opt);
end

%--------------------------------------------------------------------------
function [data, meta]=load_one(loc, opt)

data = struct();
meta = struct();
inputdata=struct('data',data,'meta',meta,'opt',opt);

% Load groups and datasets
try
    [status,count,inputdata] = H5L.iterate(loc,opt.order,'H5_ITER_INC',0,@group_iterate,inputdata);
catch ME
    if(strcmp(opt.order,'H5_INDEX_CRT_ORDER'))
        [status,count,inputdata] = H5L.iterate(loc,'H5_INDEX_NAME','H5_ITER_INC',0,@group_iterate,inputdata);
    else
        rethrow(ME);
    end
end

data=inputdata.data;
meta=inputdata.meta;

%--------------------------------------------------------------------------
function [status, res]=group_iterate(group_id,objname,inputdata)
status=0;
attr=struct();

encodename=jsonopt('PackHex',1,inputdata.opt);

try
  data=inputdata.data;
  meta=inputdata.meta;

  % objtype index
  info = H5G.get_objinfo(group_id,objname,0);
  objtype = info.type;
  objtype = objtype+1;

  if objtype == 1
    % Group
    name = regexprep(objname, '.*/', '');

	group_loc = H5G.open(group_id, name);
	try
	  [sub_data, sub_meta] = load_one(group_loc, inputdata.opt);
	  H5G.close(group_loc);
	catch ME
	  H5G.close(group_loc);
	  rethrow(ME);
	end
	if(encodename)
        name=encodevarname(name);
    else
        name=genvarname(name);
    end
    data.(name) = sub_data;
    meta.(name) = sub_meta;

  elseif objtype == 2
    % Dataset
    name = regexprep(objname, '.*/', '');

	dataset_loc = H5D.open(group_id, name);
	try
	  sub_data = H5D.read(dataset_loc, ...
	      'H5ML_DEFAULT', 'H5S_ALL','H5S_ALL','H5P_DEFAULT');
          [status, count, attr]=H5A.iterate(dataset_loc, 'H5_INDEX_NAME', 'H5_ITER_INC', 0, @getattribute, attr);
	  H5D.close(dataset_loc);
	catch exc
	  H5D.close(dataset_loc);
	  rethrow(exc);
	end

	sub_data = fix_data(sub_data, attr, inputdata.opt);
	if(encodename)
        name=encodevarname(name);
    else
        name=genvarname(name);
    end
    data.(name) = sub_data;
    meta.(name) = attr;
  end
catch ME
    rethrow(ME);
end
res=struct('data',data,'meta',meta,'opt',inputdata.opt);

%--------------------------------------------------------------------------
function data=fix_data(data, attr, opt)
% Fix some common types of data to more friendly form.

if isstruct(data)
  fields = fieldnames(data);

  if(length(intersect(fields,{'SparseIndex',opt.complexformat{1}}))==2)
    if isnumeric(data.SparseIndex) && isnumeric(data.(opt.complexformat{1}))
      if(nargin>1 && isstruct(attr))
          if(isfield(attr,'SparseArraySize'))
              spd=sparse(1,prod(attr.SparseArraySize));
              if(isfield(data,opt.complexformat{2}))
                  spd(data.SparseIndex)=complex(data.(opt.complexformat{1}),data.(opt.complexformat{2}));
              else
                  spd(data.SparseIndex)=data.(opt.complexformat{1});
              end
              data=reshape(spd,attr.SparseArraySize(:)');
              return;
          end
      end
    end
  else
      if(numel(opt.complexformat)==2 && length(intersect(fields,opt.complexformat))==2)
        if isnumeric(data.(opt.complexformat{1})) && isnumeric(data.(opt.complexformat{2}))
            data = data.(opt.complexformat{1}) + 1j*data.(opt.complexformat{2});
        end
      else
        % if complexformat is not specified or not found, try some common complex number storage formats
        if(length(intersect(fields,{'Real','Imag'}))==2)
          if isnumeric(data.Real) && isnumeric(data.Imag)
            data = data.Real + 1j*data.Imag;
          end
        elseif(length(intersect(fields,{'real','imag'}))==2)
          if isnumeric(data.real) && isnumeric(data.imag)
            data = data.real + 1j*data.imag;
          end
        elseif(length(intersect(fields,{'Re','Im'}))==2)
          if isnumeric(data.Re) && isnumeric(data.Im)
            data = data.Re + 1j*data.Im;
          end
        elseif(length(intersect(fields,{'re','im'}))==2)
          if isnumeric(data.re) && isnumeric(data.im)
            data = data.re + 1j*data.im;
          end
        elseif(length(intersect(fields,{'r','i'}))==2)
          if isnumeric(data.r) && isnumeric(data.i)
            data = data.r + 1j*data.i;
          end
        end
      end
  end
end

if(isa(data,'uint8') || isa(data,'int8'))
  if(nargin>1 && isstruct(attr))
      if(isfield(attr,'MATLABObjectClass'))
         data=getArrayFromByteStream(data); % use undocumented function
      end
  end
end

%--------------------------------------------------------------------------
function [status, dataout]= getattribute(loc_id,attr_name,info,datain)
status=0;
attr_id = H5A.open(loc_id, attr_name, 'H5P_DEFAULT');
datain.(attr_name) = H5A.read(attr_id, 'H5ML_DEFAULT');
H5A.close(attr_id);
dataout=datain;

author: brainstorm-tools


text:
classdef PdeVariableToPrintGetter < handle

    properties (Access = private)
      physicalProblem
    end

    methods (Access = public)

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

        function v = compute(obj)
            p = obj.physicalProblem;
            if isempty(p.variables)
                p.computeChomog();
            end
            v.stress = p.variables.stress;
            v.strain = p.variables.strain;
%             v.u      = obj.splitDisplacement(p.variables.d_u,p.getDimensions().ndimf);
            v.u      = p.variables.d_u;
            v.quad   = p.getQuadrature();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.physicalProblem = cParams.physicalProblem;
        end

    end

    methods (Access = private, Static)

        function uM = splitDisplacement(u,nu)
            nnode = round(length(u)/nu);
            nodes = 1:nnode;
            uM = zeros(nnode,nu);
            for idim = 1:nu
                dofs = nu*(nodes-1)+idim;
                uM(:,idim) = u(dofs);
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [P,Uinit] = parafac_als(X,R,opts)
%PARAFAC_ALS Deprecated. Use CP_ALS instead.
%
%   See also CP_ALS.
%
%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

if (nargout == 2) & (nargin == 3)
    [P,Uinit] = cp_als(X,R,opts);
elseif (nargout == 2) & (nargin == 2)
    [P,Uinit] = cp_als(X,R);
elseif (nargout == 2) & (nargin == 1)
    [P,Uinit] = cp_als(X);
elseif (nargout == 1) & (nargin == 3)
    P = cp_als(X,R,opts);
elseif (nargout == 1) & (nargin == 2)
    P = cp_als(X,R);
elseif (nargout == 1) & (nargin == 1)
    P = cp_als(X);
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
classdef MimOrientationAxial < MimGuiPlugin
    % MimOrientationAxial. Gui Plugin for using a preset bone window/level
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimOrientationAxial is a Gui Plugin for the MIM Toolkit.
    %     The gui will create a button for the user to run this plugin.
    %     Running this plugin will change the window and level of the viewing
    %     panel to standard soft tissue values.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Axial'
        SelectedText = 'Axial'
        ToolTip = 'Changes the image orientation to a axial view'
        Category = 'View'
        Visibility = 'Dataset'
        Mode = 'Toolbar'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1
        Icon = 'axial_thumb.png'
        Location = 5
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.Orientation = GemImageOrientation.XY;
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded;
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.Orientation == GemImageOrientation.XY && ~strcmp(gui_app.ImagePanel.Mode, char(MimModes.View3DMode));
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% This function performs the binary decision tree testing.

% Copyright (C) 2012 Quan Wang <wangq10@rpi.edu>,
% Signal Analysis and Machine Perception Laboratory,
% Department of Electrical, Computer, and Systems Engineering,
% Rensselaer Polytechnic Institute, Troy, NY 12180, USA
%
% You are free to use this software for academic purposes if you cite our paper:
% Q. Wang, Y. Ou, A.A. Julius, K.L. Boyer, M.J. Kim,
% Tracking tetrahymena pyriformis cells using decision trees,
% in: 2012 International Conference on Pattern Recognition, Tsukuba Science City, Japan.
%
% For commercial use, please contact the authors.

function [y p]=decide01Tree(x,T)

% x: one data sample, row vector
% T: decision tree
% y: resulting label
% p: resulting probability of being label 1

k=1; % node index
d=1; % depth
c=1; % column

while d<T.depth
    if T.feature(k)==-1
        y=0;
        p=0;
        return;
    elseif x(T.feature(k))<=T.threshold(k)
        k=left_child(k);
        [d c]=index2depth(k);
    else
        k=right_child(k);
        [d c]=index2depth(k);
    end
end
p=T.leaf_prob(c);
if p>0.5
    y=1;
else
    y=0;
end


text:
function G = exchangenode(G,a,b)

%Exchange element at column a with element at column b;
buffer=G(:,a);
G(:,a)=G(:,b);
G(:,b)=buffer;

%Exchange element at row a with element at row b;
buffer=G(a,:);
G(a,:)=G(b,:);
G(b,:)=buffer;


text:
function s = hermite_integral ( p )

%*****************************************************************************80
%
%% HERMITE_INTEGRAL evaluates a monomial Hermite integral.
%
%  Discussion:
%
%    Integral ( -oo < x < +oo ) x^p exp(-x^2) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer P, the exponent.
%    0 <= P.
%
%    Output, real S, the value of the integral.
%
  if ( mod ( p, 2 ) == 0 )
    s = r8_factorial2 ( p - 1 ) * sqrt ( pi ) / 2.0 ^ ( p / 2 );
  else
    s = 0.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% absolute azimuthal angle to from fly to closest fly according to type
function [data,units] = compute_absanglefrom1to2(trx,n,type)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);

for i1 = 1:nflies,
  fly1 = flies(i1);
  data{i1} = abs(trx(fly1).(['anglefrom1to2_',type]));
end

units = parseunits('rad');

author: kristinbranson

repo: JAABA


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

%%
clear; close all; clc;

%% Syntax
% |[S]=cellPair2struct(fieldNameCell,fieldDataCell,convertOption,optionStruct);|

%% Description
% Converts a cell array pair to a structure. The cell |fieldNameCell|
% contains field names and the cell |fieldDataCell| contains the data. The
% data may be anything a structure can hold. The option parameter
% |convertOption| can be used to convert numeric data to text. The optional
% |optionStruct| structure can be used to control this conversion (see
% |mat2strIntDouble|).

%% Examples
%

fieldNameCell={'bob','mary','banana','qwerty_yuiop'};
fieldDataCell={5,'sure',[pi 1 3.5],[1 2 3]};

%%
% Convert all numerical data
convertOption=1;
[S]=cellPair2struct(fieldNameCell,fieldDataCell,convertOption)

%%
% No conversion
convertOption=0;
[S]=cellPair2struct(fieldNameCell,fieldDataCell,convertOption)

%%
% Convert some of the entries
convertOption=[1 0 0 1];
[S]=cellPair2struct(fieldNameCell,fieldDataCell,convertOption)

%%
%
% <<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/>.


text:
function out = SVT(X,tau)
    [U,S,V] = svd(X,'econ');
    out = U*shrink(S,tau)*V';
end

author: dynamicslab


text:
function out = DN_Quantile(y,p)
% DN_Quantile   Quantile of the data vector
%
% Calculates the quantile value at a specified proportion, p, using the
% Statistics Toolbox function, quantile.
%
%---INPUTS:
% y, the input data vector
% p, the quantile proportion

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

if nargin < 2
    fprintf(1,'Using quantile p = 0.5 (median) by default\n');
    p = 0.5;
end
if ~isnumeric(p) || (p < 0) || (p > 1)
    error('p must specify a proportion, in (0,1)');
end

out = quantile(y,p);

end

author: benfulcher

repo: hctsa


text:
function [projmat_art,D] = tomoproj_art_2(im,angles)
%TOMOPROJ_ART   [projmat_art,D] = tomoproj_art(im,angles)
%   unit of angles: Degree;
%   projection direction: When angle == 0, X-ray passes through up-down
%    direction;
%
%   Phymhan
%   02-Aug-2013 14:07:06

%Pad image
[im_pad,D] = impad(im);
%Calculate projection
num_proj = length(angles);
projmat_art = zeros(num_proj,2*D);
for k = 1:num_proj
    im_rot = imrotate(im_pad,-angles(k),'bilinear','crop');
    projmat_art(k,  1:  D) = sum(im_rot,1) ;
    projmat_art(k,D+1:2*D) = sum(im_rot,2)';
end
end

author: Sable

repo: mcbench-benchmarks


text:
function pde = elli3DcircIntf4(am,ap,bm,bp,r1,r2,n1,n2)
%% USAGE: polynomial solution for Poisson equation
%  Last Modified: 02/21/2020 by Xu Zhang

%% PDE Structure
pde = struct('intf',@intf,...
    'exactu1',@exactu1,'exactu2',@exactu2,'exactu3',@exactu3,...
    'um1',@um1,'um2',@um2,'um3',@um3,'up1',@up1,'up2',@up2,'up3',@up3,...
    'Dxu',@Dxu,'Dxum',@Dxum,'Dxup',@Dxup,'Dyu',@Dyu,...
    'Dyum',@Dyum,'Dyup',@Dyup,'Dzu',@Dzu,'Dzum',@Dzum,'Dzup',@Dzup,...
    'f1',@f1,'f2',@f2,'f3',@f3,...
    'fm1',@fm1,'fm2',@fm2,'fm3',@fm3,...
    'fp1',@fp1,'fp2',@fp2,'fp3',@fp3,...
    'A',@A,'Am',@Am,'Ap',@Ap,'one',@one,...
    'B',@B,'Bm',@Bm,'Bp',@Bp);

pde.am = am;
pde.ap = ap;
pde.bm = bm;
pde.bp = bp;
%% interface function
    function u = intf(x,y,z)
        u = (x.^2 + y.^2 + z.^2).^(1/2)/r1-1;
    end

%% exact solution
    function u = exactu1(x,y,z)
        u = um1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up1(x(id),y(id),z(id));
    end
    function u = exactu2(x,y,z)
        u = um2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up2(x(id),y(id),z(id));
    end
    function u = exactu3(x,y,z)
        u = um3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up3(x(id),y(id),z(id));
    end
    function u = um1(x,y,z)
        u = (-y + n1*uker1(x,y,z).*(y-z))/am;
    end
    function u = um2(x,y,z)
        u = (x + n1*uker1(x,y,z).*(z-x))/am;
    end
    function u = um3(x,y,z)
        u = (0 + n1*uker1(x,y,z).*(x-y))/am;
    end
    function u = up1(x,y,z)
        u = (-y + n2*uker1(x,y,z).*uker2(x,y,z).*(y-z))/ap;
    end
    function u = up2(x,y,z)
        u = (x + n2*uker1(x,y,z).*uker2(x,y,z).*(z-x))/ap;
    end
    function u = up3(x,y,z)
        u = (0 + n2*uker1(x,y,z).*uker2(x,y,z).*(x-y))/ap;
    end
% not Hcurl conforming so this one is not correct
    function u = uker1(x,y,z)
        u = r1^2 - (x.^2 + y.^2 + z.^2);
    end
    function u = uker2(x,y,z)
        u = r2^2 - (x.^2 + y.^2 + z.^2);
    end
%% Boundary Function
    function u = gD1(x,y,z)
        u = exactu1(x,y,z);
    end
    function u = gD2(x,y,z)
        u = exactu2(x,y,z);
    end
    function u = gD3(x,y,z)
        u = exactu3(x,y,z);
    end
%% Derivative of the exact solution
    function u = Dxu(x,y,z)
        u = Dxum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dxup(x(id),y(id),z(id));
    end
    function u = Dyu(x,y,z)
        u = Dyum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dyup(x(id),y(id),z(id));
    end
    function u = Dzu(x,y,z)
        u = Dzum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dzup(x(id),y(id),z(id));
    end
    function u = Dxum(x,y,z)
        u = -2*n1/am*(x.*y+x.*z-y.^2-z.^2 + uker1(x,y,z));
    end
    function u = Dyum(x,y,z)
        u = -2*n1/am*(y.*z+x.*y-z.^2-x.^2 + uker1(x,y,z));
    end
    function u = Dzum(x,y,z)
        u = -2*n1/am*(x.*z+y.*z-x.^2-y.^2 + uker1(x,y,z))+2/am;
    end
    function u = Dxup(x,y,z)
        u = -2*n2*uker2(x,y,z)./ap.*(x.*y+x.*z-y.^2-z.^2 + uker1(x,y,z)) - ...
            2*n2*uker1(x,y,z).*(x.*y+x.*z - y.^2-z.^2)/ap;
    end
    function u = Dyup(x,y,z)
        u = -2*n2*uker2(x,y,z)./ap.*(y.*z+x.*y-z.^2-x.^2 + uker1(x,y,z)) - ...
            2*n2*uker1(x,y,z).*(y.*z+x.*y - x.^2-z.^2)/ap;
    end
    function u = Dzup(x,y,z)
        u = -2*n2*uker2(x,y,z)./ap.*(x.*z+y.*z-x.^2-y.^2 + uker1(x,y,z)) - ...
            2*n2*uker1(x,y,z).*(x.*z+y.*z - x.^2-y.^2)/ap+2/ap;
    end

%% right hand side function
    function u = f1(x,y,z)
        u = fm1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp1(x(id),y(id),z(id));
    end
    function u = f2(x,y,z)
        u = fm2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp2(x(id),y(id),z(id));
    end
    function u = f3(x,y,z)
        u = fm3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp3(x(id),y(id),z(id));
    end

    function u = fm1(x,y,z)
        u = -10*n1*(z-y) + bm*um1(x,y,z);
    end
    function u = fm2(x,y,z)
        u = -10*n1*(x-z) + bm*um2(x,y,z);
    end
    function u = fm3(x,y,z)
        u = -10*n1*(y-x) + bm*um3(x,y,z);
    end
    function u = fp1(x,y,z)
        u = -10*n2*uker2(x,y,z).*(z-y) +8*n2*(z.^3-y.^3+z.*y.^2+z.*x.^2-y.*x.^2-y.*z.^2) -...
             10*n2*uker1(x,y,z).*(z-y) + bp*up1(x,y,z);
    end
    function u = fp2(x,y,z)
        u = -10*n2*uker2(x,y,z).*(x-z) +8*n2*(x.^3-z.^3+x.*z.^2+x.*y.^2-z.*y.^2-z.*x.^2) -...
             10*n2*uker1(x,y,z).*(x-z) + bp*up2(x,y,z);
    end
    function u = fp3(x,y,z)
        u = -10*n2*uker2(x,y,z).*(y-x) +8*n2*(y.^3-x.^3+y.*x.^2+y.*z.^2-x.*z.^2-x.*y.^2) -...
             10*n2*uker1(x,y,z).*(y-x) + bp*up3(x,y,z);
    end

%% Diffusion coefficient function
    function u = A(x,y,z)
        u = Am(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Ap(x(id),y(id),z(id));
    end
    function u = Am(x,y,z)
        u = am*ones(size(x));
    end
    function u = Ap(x,y,z)
        u = ap*ones(size(x));
    end

%% Mass coefficient function
    function u = B(x,y,z)
        u = Bm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Bp(x(id),y(id),z(id));
    end
    function u = Bm(x,y,z)
        u = bm*ones(size(x));
    end
    function u = Bp(x,y,z)
        u = bp*ones(size(x));
    end

%% Other function
    function u = one(x,y,z)
        u = ones(size(x));
    end
end

repo: ifem


text:
function out = isreal(f)
%ISREAL   True for real CLASSICFUN.
%   ISREAL(F) returns TRUE if F does not have an imaginary part and FALSE
%   otherwise.
%
% See also REAL, IMAG.

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

% Check if the ONEFUN of f is real:
out = isreal(f.onefun);

end

author: chebfun

repo: chebfun


text:
function test_bug1162

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_postamble ft_postamble_history ft_freqgrandaverage

timelock1 = [];
timelock1.label = {'1' '2'};
timelock1.time  = 1:5;
timelock1.dimord = 'rpt_chan_time';
timelock1.avg = randn(1,2,5);
timelock1.cfg = 'this is number 1';

timelock2 = [];
timelock2.label = {'1' '2'};
timelock2.time  = 1:5;
timelock2.dimord = 'rpt_chan_time';
timelock2.avg = randn(1,2,5);
timelock2.cfg = 'this is number 2';

cfg = [];
grandavg = ft_timelockgrandaverage(cfg, timelock1, timelock2);

% grandavg.cfg.previous{1} should be 'this is number 1'
% grandavg.cfg.previous{2} should be 'this is number 1'
disp(grandavg.cfg.previous);

assert(iscell(grandavg.cfg.previous));
assert(~isempty(grandavg.cfg.previous));
assert(~isempty(grandavg.cfg.previous{1}));
assert(~isempty(grandavg.cfg.previous{2}));


text:
function [amppts1,timepts,okflag] = ISM_RIRpow_approx(aa,room,cc,timepts,rt_type,rt_val)
%ISM_RIRpow_approx  Approximation of ISM RIR power (Lehmann & Johansson's method)
%
% [P_VEC,T_VEC,OK_FLAG] = ISM_RIRpow_approx(ALPHA,ROOM,C,T_VEC,RT_TYPE,RT_VAL)
%
% This function returns the predicted values of RIR power in P_VEC (as
% would result from ISM simulations) estimated by means of the EDC
% approximation method described in: "Prediction of energy decay in room
% impulse responses simulated with an image-source model", J. Acoust. Soc.
% Am., vol. 124(1), pp. 269-277, July 2008. The values of P_VEC are
% computed for the time points given as input in T_VEC (in sec), which is
% assumed to contain increasing values of time. The vector T_VEC (and
% corresponding vector P_VEC) will be cropped if the numerical computation
% limits are reached for the higher time values in T_VEC (for which NaNs
% are generated in P_VEC), in which case the output parameter OK_FLAG will
% be set to 0 (1 otherwise).
%
% The environmental setting is defined via the following input parameters:
%
%    ALPHA: 1-by-6 vector, corresponding to each wall's absorption
%           coefficient: [x1 x2 y1 y2 z1 z2]. Index 1 indicates wall closest
%           to the origin. E.g.: [0.5 0.5 0.45 0.87 0.84 0.32].
%  RT_TYPE: character string, measure of reverberation time used for the
%           definition of the coefficients in ALPHA. Set to either 'T60' or
%           'T20'.
%   RT_VAL: scalar, value of the reverberation time (in seconds) defined by
%           RT_TYPE. E.g.: 0.25.
%     ROOM: 1-by-3 vector, indicating the rectangular room dimensions
%           (in m): [x_length y_length z_length]. E.g.: [4 4 3].
%        C: scalar (in m/s), propagation speed of sound waves. E.g.: 343.

% Release date: November 2009
% Author: Eric A. Lehmann, Perth, Australia (www.eric-lehmann.com)
%
% Copyright (C) 2009 Eric A. Lehmann
%
% 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/>.

numradpts = length(timepts);
radpts = cc * timepts;              % radius values corresponding to time points

bxx = ( sqrt(1-aa(1))*sqrt(1-aa(2)) )^(1/room(1));
byy = ( sqrt(1-aa(3))*sqrt(1-aa(4)) )^(1/room(2));
bzz = ( sqrt(1-aa(5))*sqrt(1-aa(6)) )^(1/room(3));

if bxx==byy && byy==bzz,
    intcase = 1;
elseif bxx==byy && bxx~=bzz,
    intcase = 2;
elseif byy==bzz && bzz~=bxx,
    if bzz<bxx,     % coordinate swap x<->z
        foo = bxx; bxx = bzz; bzz = foo;
        intcase = 2;
    else
        intcase = 3;
    end
elseif bxx==bzz && bzz~=byy,
    if bzz<byy,     % coordinate swap y<->z
        foo = byy; byy = bzz; bzz = foo;
        intcase = 2;
    else
        intcase = 4;
    end
else
    intcase = 5;
    if bxx>bzz && bxx>byy,      % coordinate swap x<->z
        foo = bxx; bxx = bzz; bzz = foo;
    elseif byy>bzz && byy>bxx,	% coordinate swap y<->z
        foo = byy; byy = bzz; bzz = foo;
    end
end

amppts1 = zeros(1,numradpts);
for ss=1:numradpts,    % compute amplitude/energy estimates
    Bx = bxx^(radpts(ss)); Bx(Bx==0) = eps;
    By = byy^(radpts(ss)); By(By==0) = eps;
    Bz = bzz^(radpts(ss)); Bz(Bz==0) = eps;
    switch intcase
        case 1
            int2 = Bx;
        case 2
            int2 = (Bx-Bz) / log(Bx/Bz);
        case 3
            n1 = log(Bz/Bx);
            int2 = Bz*( expint(n1) + log(n1) + 0.5772156649 ) / n1;
        case 4
            n1 = log(Bz/By);
            int2 = Bz*( expint(n1) + log(n1) + 0.5772156649 ) / n1;
        otherwise
            n1 = log(Bz/By);
            n2 = log(Bz/Bx);
            int2 = Bz*(log(n1/n2) + expint(n1) - expint(n2)) / log(Bx/By);
    end
    amppts1(ss) = int2/radpts(ss);      % 'propto' really...
end

okflag = 1;
foo = find(isnan(amppts1),1);
if ~isempty(foo),
    amppts1 = amppts1(1:foo-1);
    timepts = timepts(1:foo-1);
    okflag = 0;
end

if nargin==6,
    switch lower(rt_type)  % offset correction
        case 't60', sl = exp(3.05*exp(-1.85*rt_val));
        case 't20', sl = exp(3.52*exp(-7.49*rt_val));
    end
    amppts1 = amppts1 ./ exp(sl*(timepts-timepts(1)));
end


text:
function [id,d] = df_query( DF, p)
% DF_QUERY Use the spatial index build with df_build to provide approximate NN
% and distance from the NN.
%
% [id,d] = df_query(DF,p)
%
% Input:
%   DF  output of df_build
%   p  query point
% Output:
%   id  index of the point closer to p (approximate)
%   d  distance from the point closer to p (approximate)

p = p-DF.MIN;
p = p./DF.S;
p = round(p)+1;

p(p<1) = 1;
p(p>DF.C') = DF.C(p>DF.C');

id = DF.N(p(1),p(2),p(3));
d  = DF.D(p(1),p(2),p(3));

end

author: alecjacobson

repo: gptoolbox


text:
 function [dim, varargout] = size(ob, varargin)
%function dim = size(ob, varargin)
% or [dim1 dim2] = size(ob)
% "size" method for this class

dim = ob.dim;

if length(varargin) == 1
	dim = dim(varargin{1});
	return
elseif length(varargin) ~= 0
	error 'bug'
end

if nargout > 1
	varargout{1} = dim(2);
	dim = dim(1);
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function n_im_array = applyOptimizedTransforms( im_array, n_transforms )
%%applyOptimizedTransforms Summary
%  Applies optimized camera path n_transforms to im_array

num_frames = size(im_array, 1);
crop_ratio = 0.8;
[height, width, ~] = size(im_array{1});
center_x = width/2; center_y = height/2;
x_length = crop_ratio*width; y_length = crop_ratio*height;

n_im_array = cell(num_frames, 1);

for k=2:num_frames
    % Transform crop window
    p1 = [center_x - x_length/2 center_y - y_length/2 1];
    p2 = [center_x - x_length/2 center_y + y_length/2 1];
    p3 = [center_x + x_length/2 center_y + y_length/2 1];
    p4 = [center_x + x_length/2 center_y - y_length/2 1];
    p1 = p1 * n_transforms{k - 1}; p2 = p2 * n_transforms{k - 1};
    p3 = p3 * n_transforms{k - 1}; p4 = p4 * n_transforms{k - 1};

    % Find homography and crop
    Xin = [p1(1) p2(1) p3(1) p4(1)];
    Yin = [p1(2) p2(2) p3(2) p4(2)];
    Xout = [1 1 width width];
    Yout = [1 height height 1];

    in = [Xin; Yin]; out = [Xout; Yout];
    tform = estimateGeometricTransform(in',out','projective');
    R = imref2d([360 640 3]);
    n_im_array{k} = imwarp(im_array{k},tform,'cubic','OutputView',R);

end

% Apply same transformation to 1st frame as in the 2nd frame
n_im_array{1} = n_im_array{2};

end

author: ishit