text:
function [C,h,Ph,F,Fa,Fc] = spm_sp_reml(YY,X,Q,N,hE)
% ReML estimation of covariance components from y*y' (for sparse patterns)
% FORMAT [C,h,Ph,F,Fa,Fc] = spm_sp_reml(YY,X,Q,N);
%
% YY  - (m x m) sample covariance matrix Y*Y'  {Y = (m x N) data matrix}
% X   - (m x p) design matrix
% Q   - {1 x q} components Q.q = eigenvectors; Q.v = eigenvalues
%               or (m x n) matrix of n basis functions
% N   - number of samples
%
% C   - (m x m) estimated errors = h(1)*Q{1} + h(2)*Q{2} + ...
% h   - (q x 1) ReML hyperparameters h
% Ph  - (q x q) conditional precision of log(h)
%
% F   - [-ve] free energy F = log evidence = p(Y|X,Q) = ReML objective
%
% Fa  - accuracy
% Fc  - complexity (F = Fa - Fc)
%
% Performs a Fisher-Scoring ascent on F to find ReML variance parameter
% estimates, using uninformative hyperpriors (this is effectively an ARD
% scheme).  The specification of components differs from spm_reml and
% spm_reml_sc.
%
%__________________________________________________________________________
%
% SPM ReML routines:
%
%      spm_reml:    no positivity constraints on covariance parameters
%      spm_reml_sc: positivity constraints on covariance parameters
%      spm_sp_reml: for sparse patterns (c.f., ARD)
%
%__________________________________________________________________________
% Copyright (C) 2006-2017 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_sp_reml.m 7192 2017-10-18 14:59:01Z guillaume $

% assume a single sample if not specified
%--------------------------------------------------------------------------
try, N; catch, N = 1;   end

% default number of iterations
%--------------------------------------------------------------------------
try, K; catch, K = 128; end

% number of hyperparameters
%--------------------------------------------------------------------------
n     = length(YY);
m     = length(Q);
h     = zeros(m,1);
dh    = zeros(m,1);
dFdh  = zeros(m,1);
dFdhh = zeros(m,m);

% uninformative hyperpriors
%--------------------------------------------------------------------------
hE    = sparse(m,1) - 32;
hC    = speye(m,m)*256;
hP    = inv(hC);

% ortho-normalise X
%--------------------------------------------------------------------------
if isempty(X)
    X = sparse(n,0);
    R = speye(n,n);
else
    X = orth(full(X));
    R = speye(n,n) - X*X';
end

% convert bssis to struct if necessary
%--------------------------------------------------------------------------
if ~iscell(Q)
    q     = cell(1,m);
    for i = 1:m
        q{i}.q = spm_en(Q(:,i));
    end
    Q     = q;
end

% find bases of Q if necessary
%--------------------------------------------------------------------------
for i = 1:m
    if isstruct(Q{i})
        try
            if ~isfield(Q{i},'v');
                Q{i}.v = ones(size(Q{i}.q,2),1);
            end
        catch
            warndlg('Please specify components with Q.q (basis) and Q.v');
            return
        end
    else
        [q,v] = spm_svd(Q{i});
        C.q   = q;
        C.v   = diag(v);
        Q{i}  = C;
    end
end

% scale YY
%--------------------------------------------------------------------------
sY    = n*trace(YY)/N;
YY    = YY/sY;

% scale Q
%--------------------------------------------------------------------------
for i = 1:m
    sh(i,1) = n*trace(Q{i}.q'*R*Q{i}.q*diag(Q{i}.v));
    Q{i}.q  = Q{i}.q/sqrt(sh(i));
end

% compute basis and dsdh
%--------------------------------------------------------------------------
q     = cell(1,m);
v     = cell(1,m);
for i = 1:m
    q{i} = Q{i}.q;
    v{i} = Q{i}.v(:);
end
q     = spm_cat(q);
dedh  = spm_cat(spm_diag(v));

% pre-compute bases
%--------------------------------------------------------------------------
[n,s] = size(q);
qq    = cell(s,1);
for i = 1:s
    qq{i} = q(:,i)*q(:,i)';
end

% ReML (EM/VB)
%--------------------------------------------------------------------------
dF    = Inf;
as    = 1:m;
for k = 1:K

    % compute current estimate of covariance
    %----------------------------------------------------------------------
    C     = sparse(n,n);
    e     = dedh*exp(h);
    for i = 1:s
        C = C + qq{i}*e(i);
    end
    iC    = inv(C + speye(n,n)/exp(32));

    % E-step: conditional covariance cov(B|y) {Cq}
    %======================================================================
    iCX    = iC*X;
    if ~isempty(X)
        Cq = inv(X'*iCX);
    else
        Cq = sparse(0);
    end

    % M-step: ReML estimate of hyperparameters
    %======================================================================

    % select relevant patterns
    %----------------------------------------------------------------------
    rel   = any(dedh(:,as),2);

    % Gradient dF/dh (first derivatives)
    %----------------------------------------------------------------------
    U     = iC - iCX*Cq*iCX';
    W     = U*(YY/N - C)*U';
    qr    = q(:,rel);
    P     = qr'*U*qr;

    % dF/dh = -trace(dF/diC*iC*Q{i}*iC)
    %----------------------------------------------------------------------
    dFde  =  N/2*sum(qr.*(W*qr))';
    dFdee = -N/2*P.*P';

    % dF/dhh = -trace{P*Q{i}*P*Q{j}}
    %----------------------------------------------------------------------
    dhdh  = dedh(rel,as)*diag(exp(h(as)));
    dFdh  = dhdh'*dFde;
    dFdhh = dhdh'*dFdee*dhdh;

    % add hyperpriors
    %----------------------------------------------------------------------
    e     = h     - hE;
    dFdh  = dFdh  - hP(as,as)*e(as);
    dFdhh = dFdhh - hP(as,as);

    % Fisher scoring: update dh = -inv(ddF/dhh)*dF/dh
    %----------------------------------------------------------------------
    dh    = spm_dx(dFdhh,dFdh)/log(k + 2);
    h(as) = h(as) + dh;

    % Convergence (1% change in log-evidence)
    %======================================================================
    % bar(h),drawnow
    dF      = dFdh'*dh;
    fprintf('%-30s: %i %30s%e\n','  ReML Iteration',k,'...',full(dF));

    % and ARD
    %----------------------------------------------------------------------
    if dF < 1e-2 || k == K
        break
    else
        as             = h > -16;
        h(~as)         = hE(~as);
        h(find(h > 1)) = 1;
    end

end

% log evidence = ln p(y|X,Q) = ReML objective = F = trace(R'*iC*R*YY)/2 ...
%--------------------------------------------------------------------------
Ph  = hP;
Ph(as,as) = Ph(as,as) - dFdhh;
if nargout > 3

    % tr(hP*inv(Ph)) - nh (complexity KL cost of parameters = 0)
    %----------------------------------------------------------------------
    Ft = trace(hP*inv(Ph)) - length(Ph);

    % complexity - KL(Ph,hP)
    %----------------------------------------------------------------------
    Fc = Ft/2 + e'*hP*e/2 + spm_logdet(Ph*inv(hP))/2;

    % Accuracy - ln p(Y|h)
    %----------------------------------------------------------------------
    Fa = Ft/2 - trace(C*U*YY*U')/2 - N*n*log(2*pi)/2 - N*spm_logdet(C)/2;

    % Free-energy
    %----------------------------------------------------------------------
    F  = Fa - Fc - N*n*log(sY)/2;

end

% return exp(h) if log-normal hyperpriors
%--------------------------------------------------------------------------
h  = sY*exp(h)./sh;
C  = sY*C;

author: spm


text:
function [P1, P2, spectrum] = slcopca(X1, X2, d, varargin)
%SLCOPCA Performs Coupled PCA Learning
%
% $ Syntax $
%   - [P1, P2] = slcopca(X1, X2, sch, ...)
%   - [P1, P2, spectrum] = slcopca(...)
%
% $ Arguments $
%   - X1:       The samples in the first modality
%   - X2:       The samples in the second modality
%   - d:        The target space dimension
%   - P1:       The projection matrix (bases) of the first modality space
%   - P2:       The projection matrix (bases) of the second modality space
%   - spectrum: The covariance energy along dimensions of target space
%
% $ Description $
%   - [P1, P2] = slcopca(X1, X2, sch, ...) performs coupled PCA learning
%     for two correlated sample spaces. The learning objective is to
%     pursue two subspaces such that they are maximally correlated. The
%     objective function is formulated as
%
%       maximize trace( P1^T * C_12 * P2 * P2^T * C_21 * P1 ) / n
%           s.t. P1^T * P1 = I, and  P2^T * P2 = I
%       where C_12 is the covariance between X1 and X2,
%             C_21 is the covariance between X2 and X1
%
%     Suppose the dimensions for the two spaces are d1 and d2 respectively,
%     and n pairs of corresponding samples are given in X1 and X2. Then X1
%     and X2 should be d1 x n and d2 x n matrices respectively.
%     You can further specify the following properties to control the
%     learning procedure:
%       - 'weights':    The weights of the samples, default = []
%       - 'mean1':      The pre-computed mean vector for X1, default = []
%                       if mean1 is set as 0, then it means that X1 has
%                       been centralized.
%       - 'mean2':      The pre-computed mean vector for X2, default = []
%
%   - [P1, P2, spectrum] = slcopca(...) also outputs the spectrum. You
%     can specify the following properties to control the type of the
%     output spectrum:
%       - 'spectype':   The type of output spectrum
%                       - 'normal': The average energies along target
%                                   dimensions.
%                       - 'ratio':  The ratio of preserved energy along
%                                   target dimensions to the total
%                                   energy.
%                       default = 'normal'.
%
% $ History $
%   - Created by Dahua Lin, on Sep 16, 2006
%

%% parse and verify input arguments

if nargin < 3
    raise_lackinput('slcopca', 3);
end

if ~isnumeric(X1) || ~isnumeric(X2) || ndims(X1) ~= 2 || ndims(X2) ~= 2
    error('sltoolbox:invalidarg', ...
        'X1 and X2 should be 2D numeric matrices');
end

[d1, n] = size(X1);
[d2, n2] = size(X2);
if n ~= n2
    error('sltoolbox:sizmismatch', ...
        'The numbers of samples in X1 and X2 do not match');
end

dmax = min(d1, d2);
if d > dmax
    error('sltoolbox:invalidarg', ...
        'The target dimension d should not exceed d1 or d2');
end

opts.weights = [];
opts.mean1 = [];
opts.mean2 = [];
opts.spectype = 'normal';
opts = slparseprops(opts, varargin{:});

w = opts.weights;
if ~isempty(w)
    if ~isequal(size(w), [1, n])
        error('sltoolbox:sizmismatch', ...
            'w should be a 1 x n row vector');
    end
end

vmean1 = opts.mean1;
vmean2 = opts.mean2;
if ~isempty(vmean1) && ~isequal(vmean1, 0) && ~isequal(size(vmean1), [d1, 1])
    error('sltoolbox:sizmismatch', ...
        'The size of mean1 is illegal');
end
if ~isempty(vmean2) && ~isequal(vmean2, 0) && ~isequal(size(vmean2), [d2, 1])
    error('sltoolbox:sizmismatch', ...
        'The size of mean1 is illegal');
end

%% main skeleton

% preprocess sample matrices

X1 = preprocess_samples(X1, vmean1, w);
X2 = preprocess_samples(X2, vmean2, w);

% construct problem
S = X1 * X2';

switch opts.spectype
    case 'normal'
        if isempty(w)
            tw = n;
        else
            tw = sum(w);
        end
    case 'ratio'
        tw = trace(S * S');
    otherwise
        error('sltoolbox:invalidarg', ...
            'Invalid spectrum type: %s', opts.spectype);
end

if d > dmax / 3;
    [P1, D, P2] = svd(S, 'econ');
    spectrum = diag(D);
    spectrum = spectrum(1:d);
    P1 = P1(:, 1:d);
    P2 = P2(:, 1:d);
else
    [P1, D, P2] = svds(S, d);
    spectrum = diag(D);
end

% post-process spectrum
if nargout >= 3
    spectrum = spectrum .* spectrum / tw;
end

%% Auxiliary functions

function Xp = preprocess_samples(X, vmean, w)

if ~isequal(vmean, 0)
    if isempty(vmean)
        vmean = slmean(X, w, true);
    end
    Xp = sladdvec(X, -vmean, 1);
else
    Xp = X;
end

if ~isempty(w)
    Xp = slmulvec(Xp, w, 2);
end

author: lmthang


text:
function [annolist_flat,single_person_flat] = flatten_annolist(annolist,single_person)

annolist_flat = struct('image',[],'annorect',[]);

n = 0;
single_person_flat = [];
for imgidx = 1:length(annolist)
    rect_gt = annolist(imgidx).annorect;
    for ridx = 1:length(rect_gt)
        if (isfield(rect_gt(ridx),'objpos') && ~isempty(rect_gt(ridx).objpos))
            n = n + 1;
            annolist_flat(n).image.name = annolist(imgidx).image.name;
            annolist_flat(n).annorect = rect_gt(ridx);
            single_person_flat(n) = ismember(ridx,single_person{imgidx});
        end
    end
end

end

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
% Author: Adapted by Jai Juneja from a function created by Ajmal Saeed
% Mian. Modified to include the following functionality:
%   - Calculation of estimated variance
%   - Initial guess input to reduce number of iterations for convergence
%   and increase likelihood of convergence
% Date Modified: 16/03/2013
%
% Original code includes the following copyright information:
% This code is written by Ajmal Saeed Mian {ajmal@csse.uwa.edu.au}
% Computer Science, The University of Western Australia. The code
% may be used, modified and distributed for research purposes with
% acknowledgement of the author and inclusion of this copyright information.
function [R, t, corr, icp_var, data2] = doICP(data1, data2, res, u, tri)
% [R, t, corr, error, data2] = icp2(data1, data2, res, u, tri)
%
% This is an implementation of the Iterative Closest Point (ICP) algorithm.
% The function takes two data sets and registers data2 with data1. It is
% assumed that data1 and data2 are in approximation registration. The code
% iterates till no more correspondences can be found.
%
% This is a modified version (12 April, 2005). It is more accurate and has
% less chances of getting stuck in a local minimum as opposed to my earlier
% version icp.m
%
% Arguments: data1 - 3 x n matrix of the x, y and z coordinates of data set 1
%            data2 - 3 x m matrix of the x, y and z coordinates of data set 2
%            res   - the tolerance distance for establishing closest point
%                     correspondences. Normally set equal to the resolution
%                     of data1
%            u     - 2 x 1 matrix in local polar co-ordinates with initial
%                    guess of rototranslation
%            tri   - optional argument. obtained by tri = delaunayn(data1');
%
% Returns: R - 3 x 3 accumulative rotation matrix used to register data2
%          t - 3 x 1 accumulative translation vector used to register data2
%          corr - p x 3 matrix of the index no.s of the corresponding points of
%                 data1 and data2 and their corresponding Euclidean distance
%          error - the mean error between the corresponding points of data1
%                  and data2 (normalized with res)
%          data2 - 3 x m matrix of the registered data2
%
%

maxIter = 30;
c1 = 0;
c2 = 1;
R = eye(2);
t = zeros(2,1);

% Determine R1 and t1 from odometry input u:

R1 = [cos(u(2)) -sin(u(2)); sin(u(2)) cos(u(2))];
t1 = [u(1); 0];

R_init = R1;
t_init = t1;

if nargin < 5
    tri = delaunayn(data1');
end

n = 0;
while c2 ~= c1
    c1 = c2;
    [corr, D] = dsearchn(data1', tri, data2');
    corr(:,2:3) = [[1 : length(corr)]' D];
    corr(find(D>2*res),:) = [];

    corr = -sortrows(-corr,3);
    corr = sortrows(corr,1);
    [B, Bi, Bj] = unique(corr(:,1));
    corr = corr(Bi,:);

    data2 = R1*data2;
    data2 = [data2(1,:)+t1(1); data2(2,:)+t1(2)];
    R = R1*R;
    t = R1*t + t1;
    c2 = length(corr);
    n = n + 1;
    if n > maxIter
        break;
    end
    % If there is no correlation revert to odometry data
    if isempty(corr)
        R = R_init;
        t = t_init;
        break;
    end
    [R1, t1] = reg(data1, data2, corr);
end

e1 = 1000001;
e2 = 1000000;
n = 0;
noChangeCount = 0;
while noChangeCount < 10
    e1 = e2;
    [corr, D] = dsearchn(data1', tri, data2');
    corr(:,2:3) = [[1 : length(corr)]' D];
    corr(D>2*res,:) = [];

    corr = -sortrows(-corr,3);
    corr = sortrows(corr,1);
    [B, Bi, Bj] = unique(corr(:,1));
    corr = corr(Bi,:);

    % If there is no correlation revert to odometry data
    if isempty(corr)
        R = R_init;
        t = t_init;
        break;
    end
    [R1 t1] = reg(data1, data2, corr);
    data2 = R1*data2;
    data2 = [data2(1,:)+t1(1); data2(2,:)+t1(2)];
    R = R1*R;
    t = R1*t + t1;
    e2 = sum(corr(:,3))/(length(corr)*res);

    n = n + 1;
    if n > maxIter
        break;
    end
    if abs(e2-e1)<res/10000
        noChangeCount = noChangeCount + 1;
    end
end
error = min(e1,e2);     % Minimum error
min_error2 = sum(corr(:,3))^2;  % Least squares minimum error
icp_var = min_error2/(2*length(corr)-3);    % Estimated variance

%-----------------------------------------------------------------
function [R1, t1] = reg(data1, data2, corr)
n = length(corr);
M = data1(:,corr(:,1));
mm = mean(M,2);
S = data2(:,corr(:,2));
ms = mean(S,2);
Sshifted = [S(1,:)-ms(1); S(2,:)-ms(2)];
Mshifted = [M(1,:)-mm(1); M(2,:)-mm(2)];
K = Sshifted*Mshifted';
K = K/n;
[U A V] = svd(K);
R1 = V*U';
if det(R1)<0
    B = eye(2);
    B(2,2) = det(V*U');
    R1 = V*B*U';
end
t1 = mm - R1*ms;

author: jaijuneja

repo: ekf-slam-matlab


text:
function headmodel = ft_headmodel_fns(seg, varargin)

% FT_HEADMODEL_FNS creates the volume conduction structure to be used
% in the FNS forward solver.
%
% Use as
%   headmodel = ft_headmodel_fns(seg, ...)
%
% Optional input arguments should be specified in key-value pairs and
% can include
%   tissuecond       = matrix C [9XN tissue types]; where N is the number of
%                      tissues and a 3x3 tensor conductivity matrix is stored
%                      in each column.
%   tissue           = see fns_contable_write
%   tissueval        = match tissues of segmentation input
%   transform        = 4x4 transformation matrix (default eye(4))
%   sens             = sensor information (for which ft_datatype(sens,'sens')==1)
%   deepelec         = used in the case of deep voxel solution
%   tolerance        = scalar (default 1e-8)
%
% Standard default values for conductivity matrix C are derived from
% Saleheen HI, Ng KT. New finite difference formulations for general
% inhomogeneous anisotropic bioelectric problems. IEEE Trans Biomed Eng.
% 1997
%
% Additional documentation available at:
% http://hunghienvn.nmsu.edu/wiki/index.php/FNS
%
% See also FT_PREPARE_VOL_SENS, FT_COMPUTE_LEADFIELD

% Copyright (C) 2011, Cristiano Micheli and Hung Dang
%
% 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$

ft_hastoolbox('fns', 1);

% get the optional arguments
tissue       = ft_getopt(varargin, 'tissue', []);
tissueval    = ft_getopt(varargin, 'tissueval', []);
tissuecond   = ft_getopt(varargin, 'tissuecond', []);
transform    = ft_getopt(varargin, 'transform', eye(4));
unit         = ft_getopt(varargin, 'unit', 'mm');
sens         = ft_getopt(varargin, 'sens', []);
deepelec     = ft_getopt(varargin, 'deepelec', []); % used in the case of deep voxel solution
tolerance    = ft_getopt(varargin, 'tolerance', 1e-8);

if isempty(sens)
  ft_error('A set of sensors is required')
end

if ispc
  ft_error('FNS only works on Linux and OS X')
end

% check the consistency between tissue values and the segmentation
vecval = ismember(tissueval,unique(seg(:)));
if any(vecval)==0
  ft_warning('Some of the tissue values are not in the segmentation')
end

% create the files to be written
try
  tmpfolder = pwd;

  cd(tempdir)
  [tmp,tname] = fileparts(tempname);
  segfile   = [tname];
  [tmp,tname] = fileparts(tempname);
  confile   = [tname '.csv'];
  [tmp,tname] = fileparts(tempname);
  elecfile = [tname '.h5'];
  [tmp,tname] = fileparts(tempname);
  exefile   = [tname '.sh'];
  [tmp,tname] = fileparts(tempname);
  datafile  = [tname '.h5'];

  % this requires the fieldtrip/fileio toolbox
  ft_hastoolbox('fileio', 1);

  % create a fake mri structure and write the segmentation on disk
  disp('writing the segmentation file...')
  mri = [];
  mri.dim = size(seg);
  mri.transform = eye(4);
  mri.seg = uint8(seg);

  cfg = [];
  cfg.datatype = 'uint8';
  cfg.coordsys  = 'ctf';
  cfg.parameter = 'seg';
  cfg.filename  = segfile;
  cfg.filetype  = 'analyze';
  ft_volumewrite(cfg, mri);

  % write the cond matrix on disk, load the default cond matrix in case not specified
  disp('writing the conductivity file...')
  condmatrix = fns_contable_write('tissue',tissue,'tissueval',tissueval,'tissuecond',tissuecond);
  csvwrite(confile,condmatrix);

  % write the positions of the electrodes on disk
  disp('writing the electrodes file...')
  pos = ft_warp_apply(inv(transform),sens.elecpos); % in voxel coordinates!

  % convert pos into int32 datatype.
  hdf5write(elecfile, '/electrodes/gridlocs', int32(pos));

  % Exe file
  efid = fopen(exefile, 'w');
  if ~ispc
    fprintf(efid,'#!/usr/bin/env bash\n');
    fprintf(efid,['elecsfwd1 -img ' segfile ' -electrodes ./' elecfile ' -data ./', ...
      datafile ' -contable ./' confile ' -TOL ' num2str(tolerance) ' \n']); %2>&1 > /dev/null
  end
  fclose(efid);

  % run the shell instructions
  dos(sprintf('chmod +x %s', exefile));
  dos(['./' exefile]);

  % FIXME: find a cleverer way to store the huge transfer matrix (vista?)
  [transfer,status] = fns_read_transfer(datafile);

  cleaner(segfile,confile,elecfile,exefile,datafile)

catch ME
  disp('The transfer matrix was not written')
  cleaner(segfile,confile,elecfile,exefile,datafile)
  cd(tmpfolder)
  rethrow(ME)
end

% start with an empty volume conductor
headmodel = [];
headmodel.tissue     = tissue;
headmodel.tissueval  = tissueval;
headmodel.transform  = transform;
headmodel.unit       = unit;
headmodel.segdim     = size(seg);
headmodel.type       = 'fns';
headmodel.transfer   = transfer;

if ~isempty(deepelec)
  headmodel.deepelec  = deepelec;
end

function cleaner(segfile,confile,elecfile,exefile,datafile)
delete([segfile '.hdr']);
delete([segfile '.img']);
delete(confile);
delete(elecfile);
delete(exefile);
delete(datafile);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% POP_EPOCH - Convert a continuous EEG dataset to epoched data by extracting
%               data epochs time locked to specified event types or event indices.
%               May also sub-epoch an already epoched dataset (if sub-epochs are
%               same size or smaller). This pop_function calls EPOCH.
% Usage:
%   >> OUTEEG = pop_epoch( EEG); % pop-up a data entry window
%   >> OUTEEG = pop_epoch( EEG, types, timelimits);
%   >> [OUTEEG, indices] = pop_epoch( EEG, types, timelimits,'key1', value1 ...);
%
% Graphic interface:
%   "Time-locking event type(s)" - [edit box] Select 'Edit > Event values'
%                to see a list of event.type values; else use the push button.
%                To use event types containing spaces, enter in single-quotes.
%                EPOCH function command line equivalent: 'typerange'
%   "..."      - [push button] scroll event types.
%   "Epoch limits" - [edit box] epoch latency range [start, end] in seconds relative
%                to the time-locking events. EPOCH function equivalent: 'timelim'
%   "Name for the new dataset" - [edit box]
%                EPOCH function equivalent: 'newname'
%   "Out-of-bounds EEG ..." - [edit box] Rejection limits ([min max], []=none).
%                EPOCH function equivalent: 'valuelim'
% Inputs:
%   EEG        - Input dataset. Data may already be epoched; in this case,
%                extract (shorter) subepochs time locked to epoch events.
%   types      - String (regular expression) or cell array of event types to time
%                lock to. The default is {} which means to extract time lock epochs
%                lock to every single event (Note: An event field called 'type' must
%                be defined in the 'EEG.event' structure). For regular
%                expression, see the REGEXP function.
%   timelim    - Epoch latency limits [start end] in seconds relative to
%                the time-locking event {default: [-1 2]}
%
% Optional inputs:
%   'eventindices'- [integer vector] Extract data epochs time locked to the
%                indexed event numbers. If type is defined as well (not
%                empty), take the intersection of the selected events.
%   'valuelim' - [min max] or [max]. Lower and upper bound latencies for
%                trial data. Else if one positive value is given, use its
%                negative as the lower bound. The given values are also
%                considered outliers (min max) {default: none}
%   'newname'  - [string] New dataset name {default: "[old_dataset] epochs"}
%   'epochinfo'- ['yes'|'no'] Propagate event information into the new
%                epoch structure {default: 'yes'}
%
% Outputs:
%   OUTEEG     - output dataset
%   indices    - indices of accepted events
%
% Authors: Arnaud Delorme and Hilit Serby, SCCN, INC, UCSD, 2001
%
% See also: eeglab, epoch

% deprecated
%   'timeunit' - Time unit ['seconds'|'points'] If 'seconds,' consider events
%                times to be in seconds. If 'points,' consider events as
%                indices into the data array. {default: 'points'}

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 01-25-02 reformated help & license -ad
% 02-13-02 introduction of 'key', val arguments -ad
% 02-13-02 rereferencing of events -ad
% 03-18-02 interface and debugging -ad
% 03-27-02 interface and debugging -ad & sm

function [EEG, indices, com] = pop_epoch( EEG, types, lim, varargin )

if nargin < 1
   help pop_epoch;
	return;
end
com = '';
indices = [];

if length(EEG) > 1 && isempty(EEG(1).event)
    error('Cannot process empty event structure for the first dataset');
end

if isempty(EEG(1).event)
    if EEG.trials > 1 && EEG.xmin <= 0 && EEG.xmax >=0
        disp('No EEG.event structure found: creating events of type ''TLE'' (Time-Locking Event) at time 0');
        EEG.event(EEG.trials).epoch = EEG.trials;
        for trial = 1:EEG.trials
            EEG.event(trial).epoch   = trial;
            EEG.event(trial).type    = 'TLE';
            EEG.event(trial).latency = -EEG.xmin*EEG.srate+1+(trial-1)*EEG.pnts;
        end
    else
        disp('Cannot epoch data with no events'); beep;
        return;
    end
end
if ~isfield(EEG(1).event, 'latency'),
    error( 'Absent latency field in event array/structure: must name one of the fields ''latency''');
end
if size(EEG(1).data,3) > 1
    epochlim = [num2str( round(EEG.xmin)) '  '  num2str(round(EEG.xmax))]; % Units in seconds as in GUI
%   epochlim = [num2str( round(EEG.xmin*1000)) '  '  num2str(round(EEG.xmax*1000))]; % Units in milliseconds
else
    epochlim = '-1 2';
end
OLDEEG = EEG;

if nargin < 3
	% popup window parameters
	% -----------------------
   promptstr    = { strvcat('Time-locking event type(s) ([]=all):', ...
                    'Select ''Edit > Event values'' to see type values.'), ...
                    'Epoch limits [start, end] in seconds:', ...
                    'Name for the new dataset:', ...
					'Out-of-bounds EEG rejection limits ([min max], []=none):'  };

   cbevent = ['tmpEEG = get(gcbf, ''userdata'');' ...
           'if ~isfield(tmpEEG(1).event, ''type'')' ...
				   '   errordlg2(''No type field'');' ...
				   'else' ...
                   '   tmpevent = tmpEEG(1).event;' ...
                   '   if isnumeric(tmpEEG(1).event(1).type),' ...
				   '        [tmps,tmpstr] = pop_chansel(unique([ tmpevent.type ]));' ...
				   '   else,' ...
                   '        [tmps,tmpstr] = pop_chansel(unique({ tmpevent.type }));' ...
                   '   end;' ...
				   '   if ~isempty(tmps)' ...
				   '       set(findobj(''parent'', gcbf, ''tag'', ''events''), ''string'', tmpstr);' ...
				   '   end;' ...
				   'end;' ...
				   'clear tmps tmpevent tmpEEG tmpv tmpstr tmpfieldnames;' ];

   geometry = { [2 1 0.5] [2 1 0.5] [2 1.5] [2 1 0.5] };
   uilist = { { 'style' 'text'       'string' 'Time-locking event type(s) ([]=all)' } ...
              { 'style' 'edit'       'string' '' 'tag' 'events' } ...
              { 'style' 'pushbutton' 'string' '...' 'callback' cbevent } ...
              { 'style' 'text'       'string' 'Epoch limits [start, end] in seconds' } ...
              { 'style' 'edit'       'string' epochlim } ...
              { } ...
              { 'style' 'text'       'string' 'Name for the new dataset' } ...
              { 'style' 'edit'       'string'  fastif(length(EEG) > 1 || isempty(EEG(1).setname), '', [ EEG(1).setname ' epochs' ]) 'enable' fastif(length(EEG) == 1, 'on', 'off') } ...
              { 'style' 'text'       'string' 'Out-of-bounds EEG limits if any [min max]' } ...
              { 'style' 'edit'       'string' '' } { } };

   result = inputgui( 'geometry', geometry, 'uilist', uilist, 'helpcom', 'pophelp(''pop_epoch'')', 'title', 'Extract data epochs - pop_epoch()', 'userdata', EEG);
   if length(result) == 0 return; end

   if strcmpi(result{1}, '[]'), result{1} = ''; end
   if ~isempty(result{1})
       if strcmpi(result{1}(1),'''')   % If event type appears to be in single-quotes, use comma
                                       % and single-quote as delimiter between event types. toby 2.24.2006
                                       % fixed Arnaud May 2006
            types = eval( [ '{' result{1} '}' ] );
       else types = parsetxt( result{1});
       end
   else types = {};
   end
   lim = eval( [ '[' result{2} ']' ] );

   args = {};
   if ~isempty( result{3} ),  args = { args{:}, 'newname', result{3} }; end
   if ~isempty( result{4} ),  args = { args{:}, 'valuelim', eval( [ '[' result{4} ']' ] ) }; end
   args = { args{:}, 'epochinfo', 'yes' };

else % no interactive inputs
    args = varargin;
end

% process multiple datasets
% -------------------------
if length(EEG) > 1
    if nargin < 2
        [ EEG, com ] = eeg_eval( 'pop_epoch', EEG, 'warning', 'on', 'params', { types lim args{:} } );
    else
        [ EEG, com ] = eeg_eval( 'pop_epoch', EEG, 'params', { types lim args{:} } );
    end
    return;
end

% create structure
% ----------------
if ~isempty(args)
   try, g = struct(args{:});
   catch, disp('pop_epoch(): wrong syntax in function arguments'); return; end
else
    g = [];
end

% test the presence of variables
% ------------------------------
try, g.epochfield; 	 	  catch, g.epochfield = 'type'; end; % obsolete
try, g.timeunit; 	 	  catch, g.timeunit = 'points'; end
try, g.verbose; 	      catch, g.verbose = 'on'; end % obsolete
try, g.newname; 	      catch, g.newname = fastif(isempty(EEG.setname), '', [EEG.setname ' epochs' ]); end
try, g.eventindices;      catch, g.eventindices = 1:length(EEG.event); end
try, g.epochinfo;         catch, g.epochinfo = 'yes'; end
try, if isempty(g.valuelim), g.valuelim = [-Inf Inf]; end; catch, g.valuelim = [-Inf Inf]; end

% transform string types into a int array of column indices
% ----------------------------------------------------------
eeglab_options;
tmpevent = EEG.event;
tmpeventlatency = [ tmpevent(:).latency ];
[tmpeventlatency Itmp] = sort(tmpeventlatency);
EEG.event = EEG.event(Itmp);  % sort by ascending time
Ievent = g.eventindices;
if ischar(EEG.data)
    EEG = eeg_checkset(EEG, 'loaddata'); % just in case
end

if ~isempty( types )
    % select the events for epoching
    % ------------------------------
    Ieventtmp = [];
    tmpevent = EEG.event;
    tmpeventtype = { tmpevent.type };
    if ischar(tmpeventtype{1}), tmpeventtype  = deblank(tmpeventtype); end
    if iscell(types)
		if ischar(EEG.event(1).type)
			for index2 = 1:length( types )
				tmpevent = types{index2};
				if ~ischar( tmpevent ), tmpevent = num2str( tmpevent ); end
                tmpEventList = strmatch(deblank(tmpevent), tmpeventtype, 'exact');
				Ieventtmp = [ Ieventtmp ; tmpEventList(:) ];
			end
		else
			for index2 = 1:length( types )
				tmpevent = types{index2};
				if ischar( tmpevent ),tmpevent = str2num( tmpevent ); end
				if isempty( tmpevent ), error('pop_epoch(): string entered in a numeric field'); end
				Ieventtmp = [ Ieventtmp find(tmpevent == [ tmpeventtype{:} ]) ];
			end
		end
    elseif ischar(types)
        disp('Using regular expression')
		if ~ischar(EEG.event(1).type)
            error('Event types must be string to use regular expression');
        end
        for iEvent = 1:length(EEG.event)
            if ~isempty(regexp(EEG.event(iEvent).type, types, 'once'))
                Ieventtmp = [ Ieventtmp iEvent ];
            end
        end
    else
        error('pop_epoch(): multiple event types must be entered as {''a'', ''cell'', ''array''} or a string');
    end
    Ievent = sort(intersect(Ievent, Ieventtmp));
end

% select event latencies for epoching
%------------------------------------
Ievent = sort(Ievent);
alllatencies = tmpeventlatency(Ievent);

if isempty(alllatencies)
   error('pop_epoch(): empty epoch range (no epochs were found).'); return;
end
fprintf('pop_epoch():%d epochs selected\n', length(alllatencies));

try
    % ----------------------------------------------------
    % For AMICA probabilities...Temporarily add model probabilities as channels
    %-----------------------------------------------------
    if isfield(EEG.etc, 'amica') && ~isempty(EEG.etc.amica) && isfield(EEG.etc.amica, 'v_smooth') && ~isempty(EEG.etc.amica.v_smooth) && ~isfield(EEG.etc.amica,'prob_added')
        if isfield(EEG.etc.amica, 'num_models') && ~isempty(EEG.etc.amica.num_models)
            if size(EEG.data,2) == size(EEG.etc.amica.v_smooth,2) && size(EEG.data,3) == size(EEG.etc.amica.v_smooth,3) && size(EEG.etc.amica.v_smooth,1) == EEG.etc.amica.num_models

                EEG = eeg_formatamica(EEG);
                %--------------------
                [EEG indices com] = pop_epoch(EEG,types,lim,args{:});
                %---------------------------------

                EEG = eeg_reformatamica(EEG);
                EEG = eeg_checkamica(EEG);
                return;
            else
                disp('AMICA probabilities not compatible with size of data, model probabilities cannot be epoched...')

            end
        end
    end
    % ----------------------------------------------------
catch
    warnmsg = strcat('your dataset contains amica information, but the amica plugin is not installed.  Continuing and ignoring amica information.');
    warning(warnmsg)
end

% change boundaries in rare cases when limits do not include time-locking events
% ------------------------------------------------------------------------------
tmpevents = EEG.event;
if lim(1) > 0
   % go through all onset latencies
   for Z1 = length(alllatencies):-1:1
       % if there is any event in between trigger and epoch onset which are boundary events
       selEvt = find([tmpevents.latency] > alllatencies(Z1) & [tmpevents.latency] < alllatencies(Z1) + lim(1) * EEG.srate);
       selEvt = selEvt( eeg_findboundaries(tmpevent(selEvt) ) ); % keep only boundary events

      if any(selEvt)
          if sum([tmpevents(selEvt).duration]) > lim(1) * EEG.srate
              alllatencies(Z1) = [];
          else
              % correct the latencies by the duration of the data that were cutout
              alllatencies(Z1) = alllatencies(Z1) - sum([tmpevents(selEvt).duration]);
          end
      end
   end
end
if lim(2) < 0
   % go through all onset latencies
   for Z1 = length(alllatencies):-1:1
      % if there is any event in between trigger and epoch onset which are boundary events
      selEvt = find([tmpevents.latency] < alllatencies(Z1) & [tmpevents.latency] > alllatencies(Z1) + lim(2) * EEG.srate);
      selEvt = selEvt( eeg_findboundaries(tmpevent(selEvt) ) ); % keep only boundary events
      if any(selEvt)
          if sum([tmpevents(selEvt).duration]) > -lim(2) * EEG.srate
              alllatencies(Z1) = [];
          else
              % correct the latencies by the duration of the data that were cutout
              alllatencies(Z1) = alllatencies(Z1) + sum([tmpevents(selEvt).duration]);
          end
      end
   end
end

% select event time format and epoch
% ----------------------------------
switch lower( g.timeunit )
 case 'points',	[EEG.data, tmptime, indices, epochevent]= epoch(EEG.data, alllatencies, [lim(1) lim(2)]*EEG.srate, ...
                                                    'valuelim', g.valuelim, 'allevents', tmpeventlatency);
  tmptime = tmptime/EEG.srate;
 case 'seconds',	[EEG.data, tmptime, indices, epochevent]= epoch(EEG.data, alllatencies, lim, 'valuelim', g.valuelim, ...
                                                    'srate', EEG.srate, 'allevents', tmpeventlatency);
	otherwise, disp('pop_epoch(): invalid event time format'); beep; return;
end
alllatencies = alllatencies(indices);
fprintf('pop_epoch():%d epochs generated\n', length(indices));

% update other fields
% -------------------
if lim(1) ~= tmptime(1) && lim(2)-1/EEG.srate ~= tmptime(2)
	fprintf('pop_epoch(): time limits have been adjusted to [%3.3f %3.3f] to fit data points limits\n', ...
		tmptime(1), tmptime(2)+1/EEG.srate);
end
EEG.xmin = tmptime(1);
EEG.xmax = tmptime(2);
EEG.pnts = size(EEG.data,2);
EEG.trials = size(EEG.data,3);
EEG.icaact = [];
if ~isempty(EEG.setname)
	if ~isempty(EEG.comments)
		EEG.comments = strvcat(['Parent dataset "' EEG.setname '": ----------'], EEG.comments);
	end
	EEG.comments = strvcat(['Parent dataset: ' EEG.setname ], ' ', EEG.comments);
end
EEG.setname = g.newname;

% count the number of events to duplicate and duplicate them
% ----------------------------------------------------------
totlen = 0;
for index=1:EEG.trials, totlen = totlen + length(epochevent{index}); end
EEG.event(1).epoch = 0;          % create the epoch field (for assignment consistency afterwards)
if totlen ~= 0
    newevent(totlen) = EEG.event(1); % reserve array
else
    newevent = [];
end

% modify the event structure accordingly (latencies and add epoch field)
% ----------------------------------------------------------------------
allevents = [];
count = 1;
for index=1:EEG.trials
    for indexevent = epochevent{index}
		newevent(count)         = EEG.event(indexevent);
        newevent(count).epoch   = index;
	    newevent(count).latency = newevent(count).latency ...
			- alllatencies(index) - tmptime(1)*EEG.srate + 1 + EEG.pnts*(index-1);
		count = count + 1;
    end
end
EEG.event = newevent;
EEG.epoch = [];
EEG.saved = 'no';
EEG = eeg_checkset(EEG, 'eventconsistency');

% check for boundary events
% -------------------------
disp('pop_epoch(): checking epochs for data discontinuity');
if ~isempty(EEG.event) && ischar(EEG.event(1).type)
    tmpevent = EEG.event;
    boundaryindex = eeg_findboundaries( tmpevent );
	if ~isempty(boundaryindex)
		indexepoch = [];
		for tmpindex = boundaryindex
            if isfield(tmpevent, 'epoch')
    			indexepoch = [indexepoch tmpevent(tmpindex).epoch ];
            else
                indexepoch = 1; % only one epoch
            end
		end
		EEG = pop_select(EEG, 'notrial', indexepoch);
        % update the "indices of accepted events", too
        indices = indices(setdiff(1:length(indices),indexepoch));
	end
end

% generate text command
% ---------------------
if ischar(types)
    com = sprintf('EEG = pop_epoch( EEG, ''%s'', [%s]', types, num2str(lim));
else
    com = sprintf('EEG = pop_epoch( EEG, { ');
    for j=1:length(types)
        if ischar( types{j} )  com = sprintf('%s ''%s'' ', com, types{j} );
        else                    com = sprintf('%s [%s] ',   com, num2str(types{j}) );
        end
    end
    com = sprintf('%s }, [%s]', com, num2str(lim));
end
for i=1:2:length(args)
    if ~isempty( args{i+1} )
        if ischar( args{i+1} )   com = sprintf('%s, ''%s'', ''%s''', com, args{i}, args{i+1} );
        else                    com = sprintf('%s, ''%s'', [%s]', com, args{i}, num2str(args{i+1}) );
        end
    end
end
com = [com ');'];
return; % text command

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [iPlot,minDist] = BF_ClosestPoint_ginput(xy,inputPoint,doNormalize)
% ClosestPoint_ginput   The closest point in a dataset to the input co-ordinates given.
%
%---INPUTS:
% xy, A Nx2 vector of x-y co-ordinates.
% inputPoint, The output of a ginput.
% doNormalize, whether to normalize the input space so that relative differences
%               in x and y are treated (rather than raw Euclidean distances).
%
%---OUTPUT:
% plotMe, the closest point.

% ------------------------------------------------------------------------------
% 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 work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 4.0 International License. To view a copy of
% this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/ or send
% a letter to Creative Commons, 444 Castro Street, Suite 900, Mountain View,
% California, 94041, USA.
% ------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
% Check inputs:
%-------------------------------------------------------------------------------
if nargin < 2 || isempty(inputPoint)
    inputPoint = ginput(1);
end
if nargin < 3
    doNormalize = true;
end

%-------------------------------------------------------------------------------
% Normalize
%-------------------------------------------------------------------------------
% Normalizing the space to a unit square treats distances in both axes similarly
% relative to their ranges.
if doNormalize
    minMaxx = [min(xy(:,1)),max(xy(:,1))];
    minMaxy = [min(xy(:,2)),max(xy(:,2))];
    xy(:,1) = (xy(:,1) - minMaxx(1))/(minMaxx(2)-minMaxx(1));
    inputPoint(1) = (inputPoint(1) - minMaxx(1))/(minMaxx(2)-minMaxx(1));
    xy(:,2) = (xy(:,2) - minMaxy(1))/(minMaxy(2)-minMaxy(1));
    inputPoint(2) = (inputPoint(2) - minMaxy(1))/(minMaxy(2)-minMaxy(1));
end

%-------------------------------------------------------------------------------
% Calculate distances from each point to input inputPoint
%-------------------------------------------------------------------------------
dpxy = sum((xy - repmat(inputPoint,size(xy,1),1)).^2,2); % Euclidean distances to the input point
[minDist,iPlot] = min(dpxy);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [idx,idy] = convertind(ind,m,n)
    if mod(ind,m) == 0
        idx = m;idy = ind/m;
    else
        idx = mod(ind,m);idy = floor(ind/m)+1;
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [feat, unique_box_id] = spp_features_convX_to_poolX(spp_pooler, feat, boxes, random_scale, dedup)
% [feat, unique_box_id] = spp_features_convX_to_poolX(feat, boxes, random_scale, dedup)
%
%   boxes           in rows with (l, t, r, b)
%   random_scale    find best scale for each box or not
%   dedup           remove duplicate boxes on response map
%
% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Shaoqing Ren
%
% This file is part of the SPP code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

if nargin < 4
    random_scale = false;
end

if nargin < 5
    dedup = false;
end

assert(isfield(feat, 'rsp') && ~isempty(feat.rsp));

min_img_sz = min(feat.im_height, feat.im_width);

if isempty(boxes)
    feat = [];
    return;
end

feat.scale = feat.scale(:)';
if ~random_scale
    expected_scale = spp_pooler.expected_scale(min_img_sz, boxes, spp_pooler);

    [~, best_scale_ids] = min(abs(bsxfun(@minus, feat.scale, expected_scale(:))), [], 2);
else
    best_scale_ids = randi(length(feat.scale), size(boxes, 1), 1);
end

boxes_scales = feat.scale(best_scale_ids(:));
scaled_boxes = bsxfun(@times, (boxes - 1), (boxes_scales(:) - 1)) / (min_img_sz - 1) + 1;

if dedup
    rep_boxes = spp_pooler.response_boxes(scaled_boxes, spp_pooler);
    rsp_keys = [rep_boxes, best_scale_ids];
    [~, ia] = unique(rsp_keys, 'rows');
    unique_box_id = sort(ia);
else
    unique_box_id = 1:size(boxes, 2);
end

feat = spm_pool(feat.rsp, spp_pooler.spm_divs, scaled_boxes', best_scale_ids, ...
    spp_pooler.offset0, spp_pooler.offset, spp_pooler.step_standard);

end

author: ShaoqingRen


text:
function [orderedPoints,biggest]=orderMeshPerimeterPointsAll(mesh,perimeterEdges)
% Take a set of mesh points on a perimeter and order them by their connections
%
%  orderedPoints=orderMeshPerimeterPointsAll(perimeterEdges)
%
% If more than one (unlinked) perimeter exists, return a structure
% containing all the separate ordered perimeters.
%
% N.B.  Each vertex should appear in exactly two rows.
% Starting with an initial position on the perimeter, startNode, the routine finds
% the other node in the row. By selectively deleting entries in the connection matrix
% as we progress, we can use a simple find operation each time to give us the next node.
%
% AUTHOR:  Wade
% DATE: 200-06-21

foundAllPerimeters=0; % Flag to say we've finished
nVerts=length(mesh.connectionMatrix); % Total number of vertices in the full connection matrix

edgePoints=unique(perimeterEdges(:)); % List of nodes that are on the perimeter
[y x]=size(mesh.connectionMatrix); % = nVerts

% Generate a connection matrix using only the edges - make sure it's symmetric
ec=sparse(perimeterEdges(:,1),perimeterEdges(:,2),1,nVerts,nVerts);
ec=ec+ec';
edgeConMat=(ec~=0);
%nnz(edgeConMat)

% Used for checking which perimeter is biggest
biggest=-Inf;
biggestSize=-Inf;

thisPerim=1; % Index of the current perimeter
nextRow=-99999; % Dummy

while (~foundAllPerimeters)
    counter=1;
    [startRow startCol]=find(edgeConMat);
    if (~isempty(startRow)) % If there is one...

        startRow=startRow(1); % ...pick the first non-zero row in the edge connection matrix
        startCol=startCol(1); % Pick the first non-zero column.
        orderedPoints{thisPerim}.points=startRow;

        thisRow=startRow;
        deleteCol=startCol;
        nextRow=-Inf;
        foundEnd=0;
        edgeConMat(thisRow,deleteCol)=0; % Zero one of the entries in this row

        while ((~foundEnd)&(counter<10000))

            % Do the deletions
            edgeConMat(thisRow,deleteCol)=0; % Zero one of the entries in this row
            edgeConMat(deleteCol,thisRow)=0; % Zero one of the entries in this row

            [thisCol]=find(edgeConMat(thisRow,:));

            thisCol=unique(thisCol(:));
            if (length(thisCol)>1)
                disp ('Warning!! - thiscol has more than 1 entry:');
                disp (thisCol);
                disp ('Perim node is connected to many other perim nodes. Carrying on but expect trouble....');
                % Can we just ignore superfluous ones? We'll press ahead regardless...
                thisCol=thisCol(1);
            end

            if (~isempty(thisCol)) % If there >is< another point in this perimeter
                orderedPoints{thisPerim}.points=[orderedPoints{thisPerim}.points;thisCol];
                deleteCol=thisRow;
                thisRow=thisCol;

            else % We've found the end of the perimeter.
                % fprintf ('\n.');

                foundEnd=1;
                orderedPoints{thisPerim}.size=counter;

            end % endif
            counter=counter+1;

        end % end while not startRow

        if (orderedPoints{thisPerim}.size>biggestSize) % Have to return the index of the larges tperim
            biggestSize=orderedPoints{thisPerim}.size;
            biggest=thisPerim;
        end % end size check

        fprintf('%d points\n');

    end % End if anything found to start with

    %keyboard;
    %nnz(edgeConMat)
    foundAllPerimeters=(x*y)-nnz(edgeConMat); %(sum(sum(edgeConMat==0)));
    thisPerim=thisPerim+1;
    %     disp(thisPerim);
    %     disp(nnz(edgeConMat));

end % End while not all perims found

fprintf('Found %d perimeters\n',thisPerim-1);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [sinx, cosx] = SinCosNorm(sinx, cosx)
%SINCOSNORM  Normalize sinx and cosx
%
%   [SINX, COSX] = SINCOSNORM(SINX, COSX) normalize SINX and COSX so that
%   SINX^2 + COSX^2 = 1.  SINX and COSX can be any shape.

  r = hypot(sinx, cosx);
  sinx = sinx ./ r;
  cosx = cosx ./ r;
end

author: Sable

repo: mcbench-benchmarks


text:
function [mResult] = calc_GridMcBboot(vResults, nMethod)
% function [mResult] = calc_GridMcBboot(vResults, nMethod)
% --------------------------------------------------------
% Calculate Mc and b-value using the bootstrap mean value for an already
% existing grid with another Mc method
%
% Incoming variables:
% vResults : Struct array from grid calculated with sv_calcMc / sv_calc
% nMethod  : Method to calculate Mc
%
% Outgoing variables:
% mResult : [nNodeGridPoint fMc fStd_Mc fBvalue fStd_B fAvalue fStd_A]
%          Standard deviations by bootstrap
%          fMc, fBvalue, fAvalue are mean values from the bootstrap
%
% J. Woessner: woessner@seismo.ifg.ethz.ch
% last update: 25.11.03

mResult = [];

for nNodeGridPoint=1:length(vResults.mPolygon(:,1))
    % Get the data for the grid node
    mNodeCatalog_ = vResults.mCatalog(vResults.caNodeIndices{nNodeGridPoint}, :);
    % Create the frequency magnitude distribution
    [vFMD, vNonCFMD] = calc_FMD(mNodeCatalog_);
    [nY,nX]=size(mNodeCatalog_);
    if nY > vResults.nMinimumNumber
       [fMc, fStd_Mc, fBvalue, fStd_B, fAvalue, fStd_A, vMc, mBvalue] = calc_McBboot(mNodeCatalog_,vResults.fBinning, vResults.fBstnum, nMethod);
       mResult = [mResult; nNodeGridPoint fMc fStd_Mc fBvalue fStd_B fAvalue fStd_A];
    else
       mResult = [mResult; nNodeGridPoint NaN NaN NaN NaN NaN NaN];
    end
    if rem(nNodeGridPoint,500) == 0
        save(['result' num2str(nNodeGridPoint) '.mat'],'mResult');
    end
end
save(['result' num2str(nNodeGridPoint) '.mat'],'mResult');

author: CelsoReyes


text:
% POP_SUBCOMP - remove specified components from an EEG dataset.
%                 and subtract their activities from the data. Else,
%                 remove components already marked for rejection. When used
%                 with the options 'keepcomp', the function will retain [1]
%                 or reject[0] the components provided as input.
% Usage:
%   >> OUTEEG = pop_subcomp( INEEG ); % pop-up window mode
%   >> OUTEEG = pop_subcomp( INEEG, components, plotag);
%   >> OUTEEG = pop_subcomp( INEEG, components, plotag, keepcomp);
%
% Pop-up window interface:
%   "Component(s) to remove ..." - [edit box] Array of components to
%                remove from the data. Sets the 'components' parameter
%                in the command line call (see below).
%   "Component(s) to retain ..." - [edit box] Array of components to
%                to retain in the data. Sets the 'components' parameter in
%                the command line call. Then, comp_to_remove = ...
%                    setdiff([1:size(EEG.icaweights,1)], comp_to_keep). See
%                    option 'keepcomp' for command line call.
%                Overwrites "Component(s) to remove" (above).
% Command line inputs:
%   INEEG      - Input EEG dataset.
%   components - Array of components to remove from the data. If empty,
%                 remove components previously marked for rejection (e.g.,
%                 EEG.reject.gcompreject).
%   plotag     - [0|1] Display the difference between original and processed
%                dataset. 1 = Ask for confirmation. 0 = Do not ask. {Default: 0}
%   keepcomp   - [0|1] If [1] will retain the components provided by the
%                input variable 'components', [0] will reject them. Option
%                intended to be used only with the components provided in
%                'components', and not with components marked for rejection
%                {Default: 0}
% Outputs:
%   OUTEEG     - output dataset.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: COMPVAR

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 01-25-02 reformated help & license -ad
% 02-15-02 propagate ica weight matrix -ad sm jorn

function [EEG, com] = pop_subcomp( EEG, components, plotag, keepflag)

com='';
if nargin < 1
   help pop_subcomp;
   return;
end
if nargin < 3
	plotag = 0;
end
if nargin < 4
    keepflag = 0;
end

if nargin < 2
    components = [];
	% popup window parameters
	% -----------------------
    res = 'Manual rej.';
    if any(cellfun(@(x)any(x.gcompreject), { EEG.reject }))
        if length(EEG) == 1
            compStr = sprintf('%d,', find(EEG.reject.gcompreject == 1));
            msg = sprintf('Components [%s] flagged for rejection.', compStr(1:end-1));
            res = questdlg2(strvcat(msg, 'Do you want to remove these components?', 'Note: we recommend removing components in STUDY instead'), 'Remove components from data', 'Cancel', 'Manual rej.', 'Yes', 'Cancel');
        else
            msg = 'Components flagged for rejection detected in some datasets.';
            res = questdlg2(strvcat(msg, 'Do you want to remove these components?', 'Note: we recommend removing components in STUDY instead'), 'Remove components from data', 'Cancel', 'Yes', 'Cancel');
        end
        if strcmpi(res, 'Cancel')
            return;
        end
        if strcmpi(res, 'Yes')
            components = '';
        end
    else
        if length(EEG) > 1
            warndlg2(strvcat('You have multiple datasets selected and no components', ...
                'flagged for rejection. Flag components first.'));
            return;
        end
    end

    if strcmpi(res, 'Manual rej.')
        if ~isempty(EEG.reject.gcompreject)
            components = find(EEG.reject.gcompreject == 1);
            components = components(:)';
            %promptstr    = { ['Components to subtract from data' 10 '(default: pre-labeled components to reject):'] };
        else
            components = [];
        end
        uilist    = { { 'style' 'text' 'string' 'Note: for group level analysis, remove components in STUDY' } ...
                      { 'style' 'text' 'string' 'List of component(s) to remove from data' } ...
                      { 'style' 'edit' 'string' int2str(components) } ...
                      { 'style' 'text' 'string' 'Or list of component(s) to retain' } ...
                      { 'style' 'edit' 'string' '' } ...
                      };
        geom = { 1 [2 0.7] [2 0.7] };
        result       = inputgui( 'uilist', uilist, 'geometry', geom, 'helpcom', 'pophelp(''pop_subcomp'')', ...
                                         'title', 'Remove components from data -- pop_subcomp()');
        if length(result) == 0 return; end
        components   = eval( [ '[' result{1} ']' ] );
        if ~isempty(result{2})
            componentsOld = components;
            components   = eval( [ '[' result{2} ']' ] );
            if isequal(components, componentsOld)
                components = [];
            end
            keepflag = 1; %components  = setdiff_bc([1:size(EEG.icaweights,1)], components);
        end
    end
end

% process multiple datasets
% -------------------------
if length(EEG) > 1
    if nargin < 2
        [ EEG, com ] = eeg_eval( 'pop_subcomp', EEG, 'params', { components, plotag, keepflag }, 'warning', 'on' );
    else
        [ EEG, com ] = eeg_eval( 'pop_subcomp', EEG, 'params', { components, plotag, keepflag } );
    end
    if isempty( components )
        com = [ com ' % [] or '' means removing components flagged for rejection' ];
    end
    return;
end

componentsOri = components;
if isempty(components)
    if ~isempty(EEG.reject.gcompreject)
        components = find(EEG.reject.gcompreject == 1);
    else
        fprintf('Warning: no components specified, no rejection performed\n');
        return;
    end
else
    if keepflag == 1; components  = setdiff_bc([1:size(EEG.icaweights,1)], components); end
    if (max(components) > size(EEG.icaweights,1)) || min(components) < 1
        error('Component index out of range');
    end
end

fprintf('Computing projection and removing %d components ....\n', length(components));
component_keep = setdiff_bc(1:size(EEG.icaweights,1), components);
compproj = EEG.icawinv(:, component_keep)*eeg_getdatact(EEG, 'component', component_keep, 'reshape', '2d');
compproj = reshape(compproj, size(compproj,1), EEG.pnts, EEG.trials);

%fprintf( 'The ICA projection accounts for %2.2f percent of the data\n', 100*varegg);

if nargin < 2 || plotag ~= 0

    ButtonName = 'continue';
    while ~strcmpi(ButtonName, 'Cancel') && ~strcmpi(ButtonName, 'Accept')
        ButtonName=questdlg2( [ 'Please confirm your choice. Are you sure you want to remove the selected components from the data?' ], ...
                             'Confirmation', 'Cancel', 'Plot ERPs', 'Plot single trials', 'Accept', 'Accept');
        if strcmpi(ButtonName, 'Plot ERPs')
            if EEG.trials > 1
                tracing  = [ squeeze(mean(compproj,3)) squeeze(mean(EEG.data(EEG.icachansind,:,:),3))];
                figure;
                plotdata(tracing, EEG.pnts, [EEG.xmin*1000 EEG.xmax*1000 0 0], ...
                    'Trial ERPs (blue) with and (red) without these components');
            else
                warndlg2('Cannot plot ERPs for continuous data');
            end
        elseif strcmpi(ButtonName, 'Plot single trials')
        	eegplot( EEG.data(EEG.icachansind,:,:), 'srate', EEG.srate, 'title', 'Black = channel before rejection; red = after rejection -- eegplot()', ...
            	 'limits', [EEG.xmin EEG.xmax]*1000, 'data2', compproj);
        end
    end
    switch ButtonName
        case 'Cancel'
        	disp('Operation cancelled');
        	return;
        case 'Accept'
       		disp('Components removed');
    end % switch
end
EEG.data(EEG.icachansind,:,:) = compproj;
EEG.setname = [ EEG.setname ' pruned with ICA'];
EEG.icaact  = [];
goodinds    = setdiff_bc(1:size(EEG.icaweights,1), components);
EEG.icawinv     = EEG.icawinv(:,goodinds);
EEG.icaweights  = EEG.icaweights(goodinds,:);
EEG.specicaact  = [];
EEG.specdata    = [];
EEG.reject      = [];

if isfield(EEG.etc, 'ic_classification')
    if isfield(EEG.etc.ic_classification, 'ICLabel')
        if isfield(EEG.etc.ic_classification.ICLabel, 'classifications')
            if ~isempty(EEG.etc.ic_classification.ICLabel.classifications)
                EEG.etc.ic_classification.ICLabel.classifications = EEG.etc.ic_classification.ICLabel.classifications(goodinds,:);
            end
        end
    end
end

try
    EEG.dipfit.model = EEG.dipfit.model(goodinds);
catch, end

com = sprintf('EEG = pop_subcomp( EEG, [%s], %d);', int2str(componentsOri(:)'), plotag);
if isempty( components )
    com = [ com ' % [] means removing components flagged for rejection' ];
end
return;


text:
function diffS = comparePyradWithIBSI1OrigImgWithInterp
% Compare pyradiomics features against IBSI benchmark for config C.
%--------------------------------------------------------------------------
% AI 7/1/2020

%% Calc. features for configuration 'C' using Pyradiomics
fpath = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing/data_for_cerr_tests/IBSI1_CT_phantom/IBSILungCancerCTImage.mat.bz2');
planC = loadPlanC(fpath,tempdir);
planC = updatePlanFields(planC);
planC = quality_assure_planC(fpath,planC);
strName='GTV-1';

configPath_avg = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing/settings_for_comparisons/pyRadConfigC_avg.yaml');
pyFeat1S = calcRadiomicsFeatUsingPyradiomics(planC,strName,configPath_avg);
configPath_merge = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing/settings_for_comparisons/pyRadConfigC_merge.yaml');
pyFeat2S = calcRadiomicsFeatUsingPyradiomics(planC,strName,configPath_merge);

%% Compare with IBSI

% Get IBSI bechmark
ibsiConfigCResult = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing/data_for_cerr_tests/IBSI1_CT_phantom/IBSI_results_configC.mat');
temp = load(ibsiConfigCResult);
IBSIfeatS = temp.IBSIfeatS;

% Shape features
pyShapeS = getPyradFeatDict(pyFeat1S,{'original_shape'});
pyShapeS = mapPyradFieldnames(pyShapeS,'original','shape');
shapeFeatC = fieldnames(IBSIfeatS.shapeS);
for n = 1:length(shapeFeatC)
    ibsiVal = IBSIfeatS.shapeS.(shapeFeatC{n});
    if isfield(pyShapeS,shapeFeatC{n})
        pyRadVal = pyShapeS.(shapeFeatC{n});
        pctDiff = (pyRadVal-ibsiVal)*100/ibsiVal;
        diffS.Shape.(shapeFeatC{n}) = pctDiff;
    end
end

% First order features
pyFirstOrdS = getPyradFeatDict(pyFeat1S,{'original_firstorder'});
pyFirstOrdS = mapPyradFieldnames(pyFirstOrdS,'original','firstorder');
%Convert kurtosis to excess kurtosis
pyFirstOrdS.kurtosis = pyFirstOrdS.kurtosis -3;
firstOrdFeatC = fieldnames(IBSIfeatS.Original.firstOrderS);
for n = 1:length(firstOrdFeatC)
    ibsiVal = IBSIfeatS.Original.firstOrderS.(firstOrdFeatC{n});
    if isfield(pyFirstOrdS,firstOrdFeatC{n})
        pyRadVal = pyFirstOrdS.(firstOrdFeatC{n});
        pctDiff = (pyRadVal-ibsiVal)*100/ibsiVal;
        diffS.FirstOrder.(firstOrdFeatC{n}) = pctDiff;
    end
end

% GLCM
%Avg
pyGlcmS = getPyradFeatDict(pyFeat1S,{'original_glcm'});
pyGlcmS = mapPyradFieldnames(pyGlcmS,'original','glcm');
glcmFeatC = fieldnames(IBSIfeatS.Original.glcmFeatS.AvgS);
for n = 1:length(glcmFeatC)
    ibsiVal = IBSIfeatS.Original.glcmFeatS.AvgS.(glcmFeatC{n});
    if isfield(pyGlcmS,glcmFeatC{n})
        pyRadVal = pyGlcmS.(glcmFeatC{n});
        pctDiff = (pyRadVal-ibsiVal)*100/ibsiVal;
        diffS.GLCM.Avg.(glcmFeatC{n}) = pctDiff;
    end
end
%Merge
pyGlcm2S = getPyradFeatDict(pyFeat2S,{'original_glcm'});
pyGlcm2S = mapPyradFieldnames(pyGlcm2S,'original','glcm');
glcmFeat2C = fieldnames(IBSIfeatS.Original.glcmFeatS.CombS);
for n = 1:length(glcmFeat2C)
    ibsiVal = IBSIfeatS.Original.glcmFeatS.CombS.(glcmFeat2C{n});
    if isfield(pyGlcm2S,glcmFeat2C{n})
        pyRadVal = pyGlcm2S.(glcmFeat2C{n});
        pctDiff = (pyRadVal-ibsiVal)*100/ibsiVal;
        diffS.GLCM.Merge.(glcmFeat2C{n}) = pctDiff;
    end
end

% GLRLM
%Avg
pyGlrlmS = getPyradFeatDict(pyFeat1S,{'original_glrlm'});
pyGlrlmS = mapPyradFieldnames(pyGlrlmS,'original','glrlm');
glrlmFeatC = fieldnames(IBSIfeatS.Original.rlmFeatS.AvgS);
for n = 1:length(glrlmFeatC)
    ibsiVal = IBSIfeatS.Original.rlmFeatS.AvgS.(glrlmFeatC{n});
    if isfield(pyGlrlmS,glrlmFeatC{n})
        pyRadVal = pyGlrlmS.(glrlmFeatC{n});
        pctDiff = (pyRadVal-ibsiVal)*100/ibsiVal;
        diffS.GLRLM.Avg.(glrlmFeatC{n}) = pctDiff;
    end
end
%Merge
pyGlrlm2S = getPyradFeatDict(pyFeat2S,{'original_glrlm'});
pyGlrlm2S = mapPyradFieldnames(pyGlrlm2S,'original','glrlm');
glrlm2FeatC = fieldnames(IBSIfeatS.Original.rlmFeatS.CombS);
for n = 1:length(glrlm2FeatC)
    ibsiVal = IBSIfeatS.Original.rlmFeatS.CombS.(glrlm2FeatC{n});
    if isfield(pyGlrlm2S,glrlm2FeatC{n})
        pyRadVal = pyGlrlm2S.(glrlm2FeatC{n});
        pctDiff = (pyRadVal-ibsiVal)*100/ibsiVal;
        diffS.GLRLM.Merge.(glrlm2FeatC{n}) = pctDiff;
    end
end

% NGTDM
pyNgtdmS = getPyradFeatDict(pyFeat1S,{'original_ngtdm'});
pyNgtdmS = mapPyradFieldnames(pyNgtdmS,'original','ngtdm');
ngtdmFeatC = fieldnames(IBSIfeatS.Original.ngtdmFeatS);
for n = 1:length(ngtdmFeatC)
    ibsiVal = IBSIfeatS.Original.ngtdmFeatS.(ngtdmFeatC{n});
    if isfield(pyNgtdmS,ngtdmFeatC{n})
        pyRadVal = pyNgtdmS.(ngtdmFeatC{n});
        pctDiff = (pyRadVal-ibsiVal)*100/ibsiVal;
        diffS.NGTDM.(ngtdmFeatC{n}) = pctDiff;
    end
end

% SZM features
pyGlszmS = getPyradFeatDict(pyFeat1S,{'original_glszm'});
pyGlszmS = mapPyradFieldnames(pyGlszmS,'original','glszm');
glszmFeatC = fieldnames(IBSIfeatS.Original.szmFeatS);
for n = 1:length(glszmFeatC)
    ibsiVal = IBSIfeatS.Original.szmFeatS.(glszmFeatC{n});
    if isfield(pyGlszmS,glszmFeatC{n})
        pyRadVal = pyGlszmS.(glszmFeatC{n});
        pctDiff = (pyRadVal-ibsiVal)*100/ibsiVal;
        diffS.GLSZM.(glszmFeatC{n}) = pctDiff;
    end
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function c = minus(c,I)

   % C = MINUS(C,I)
   % The function subtracts the waveform specified by trace index I from all
   % other traces. The resulting correlation object has the same number of
   % traces as the original though trace I is essentially zero amplitude.
   %
   % C = MINUS(C)
   % Short hand version of MINUS where I is assumed to be the last trace in
   % the list. That is, MINUS(C, C.ntraces ). This use is common
   % following the STACK routine.
   %
   % Before using MINUS, time windows need to be cropped to the same interval
   % relative to the triggers. This is most easily accomplished using the CROP
   % function prior to MINUS. If the traces have not been cropped, MINUS will
   % attempt to figure out which time window of data should be differenced.
   % This can result in unanticipated signals at the ends of each trace.
   % Caveat Emptor!
   %
   % MINUS will erase the CORR, LAG, STAT, LINK and CLUST fields from the
   % objects. After differencing, these fields are no longer valid.
   %
   % The traces are not normalized prior to subtracting. For most
   % applications the user will first want to normalize the trace amplitudes
   % using the NORM function. Typically, plots of the residual data will use
   % the 'raw' plotting option.
   %
   % In order to be differenced, the "phase shift" of the samples within the
   % trace must be the same. If it is not, MINUS will make a call to ALIGN to
   % resample all traces such that one sample falls directly on the trigger
   % time.
   %
   % Example - subtract the stacked waveform from each trace
   %   c = xcorr(c)
   %   c = adjusttrig(c)
   %   c = crop(c,-3,5)
   %   c = stack(c)
   %   c = norm(c)      % stacked waveform is the last trace
   %   c_residual = minus(c)
   %   plot(c_residual,'raw');

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   % READ & CHECK ARGUMENTS
   if ~exist('I','var')
      I = c.ntraces;
   end;
   assert(numel(I)==1, 'Only one trace can be subtracted from the others');

   if check(c,'SCALE')
      disp('Warning: Traces appear to have very different overall amplitudes');
   end

   % CREATE CORREATION OBJECT FOR INTERNAL MANIPULATIONS
   c1 = NewCorrelation;
   c1.traces = c.traces;
   c1.trig = c.trig;
   if exist('index','var')               % Is this necessary?
      c1 = subset(c1,index);
   end

   % CHECK IF TRACES ARE DESCRITIZED ON THE SAME INVERVALS
   t = 86400 * (c1.trig - c1.traces.firstsampletime());
   Fs = c1.samplerate;
   sampleshift = mod(t,1/Fs);
   if (mean(sampleshift) ~= sampleshift)   % if all sampleshifts are the same
      c1 = align(c1);
   end

   % CHECK IF TRACES ARE EQUALLY CROPPED
   pretrig  = c1.traces.firstsampletime() - c1.trig; %ALL? just one?
   posttrig = c1.traces.lastsampletime()  - c1.trig;
   if (mean(pretrig)~=pretrig(1)) || (mean(posttrig)~=posttrig(1))
      center = mean([pretrig ; posttrig])*86400;
      pretrig = center - (0.5 * c1.nsamples * get(c1,'Period'));
      posttrig = pretrig + (c1.nsamples + 1) * get(c1,'Period');
      % is it the right length?
      c1 = crop(c1,pretrig(1),posttrig(1)); %should all have same pre & post
   end

   toSubtract = c1.traces(I).data; %this trace will be subtracted from all
   for i = 1:length(c.trig)
      d = c1.traces(i).data - toSubtract;
      c.traces(i).data = d;
      c.trig(i) = c1.trig(i);
   end
   c.corrmatrix = [];
   c.lags = [];
   c.stat = [];
   c.link = [];
   c.clust = [];
end

author: geoscience-community-codes

repo: GISMO


text:
function [W,CovW,X,CovX,R] = rotate_to_pca(W,CovW,X,CovX,weightsWW)
% [W,CovW,X,CovX,R] = rotate_to_pca(W,CovW,X,CovX,weightsW)

% Dimensionalities
[M,D] = size(W);
[D,N] = size(X);

% Convert cell covariance matrices to arrays
cellW = false;
if iscell(CovW)
  CovW = covcell_to_covarray(CovW);
  cellW = true;
end

cellX = false;
if iscell(CovX)
  CovX = covcell_to_covarray(CovX);
  cellX = true;
end

if nargin < 5 || isempty(weightsWW)
  weightsWW = ones(M,1);
else
  weightsWW = weightsWW(:) .* ones(M,1);
end

%% Find mixing rotation Xpca = R * Xgp

R = 1;

% 1) Whiten <XX'>

XX = X * X' + sum(CovX,3);
[V,A] = svd(XX/N);
R = diag(sqrt(1./diag(A))) * V';
% Rotate W
W = W / R;
for i=1:rows(W)
  CovW(:,:,i) = R' \ CovW(:,:,i) / R;
end

% 2) Orthogonalise weighted <W'W>

% Evaluate weighted second moment
WW = W' * diag(weightsWW) * W;
for m=1:M
  WW = WW + weightsWW(m) * CovW(:,:,m);
end

% Diagonalise it
[V,D] = svd(WW);
R = V' * R;

% Rotate W
W = W * V;
for m=1:M
  CovW(:,:,m) = V' * CovW(:,:,m) * V;
end

% $$$ % Debug the rotation!
% $$$ weights2 = cosd(data.coordinates(2,:));
% $$$ WW = W' * diag(weights2) * W;
% $$$ for m=1:M
% $$$   WW = WW + weights2(m) * CovW(:,:,m);
% $$$ end
% $$$ XX = R * (X*X'+sum(CovX,3)) * R';
% $$$ WW(1:10,1:10)
% $$$ XX(1:10,1:10) / size(X,2)

% Rotate X to PCA
X = R * X;
for n=1:N
  CovX(:,:,n) = R * CovX(:,:,n) * R';
end

% Convert back to cells
if cellW
  CovW = covarray_to_covcell(CovW);
end
if cellX
  CovX = covarray_to_covcell(CovX);
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% [PYR, INDICES] = buildLpyr(IM, HEIGHT, FILT1, FILT2, EDGES)
%
% Construct a Laplacian pyramid on matrix (or vector) IM.
%
% HEIGHT (optional) specifies the number of pyramid levels to build. Default
% is 1+maxPyrHt(size(IM),size(FILT));
%
% FILT1 (optional) can be a string naming a standard filter (see
% namedFilter), or a vector which will be used for (separable)
% convolution.  Default = 'binom5'.  FILT2 specifies the "expansion"
% filter (default = filt1).  EDGES specifies edge-handling, and
% defaults to 'reflect1' (see corrDn).
%
% PYR is a vector containing the N pyramid subbands, ordered from fine
% to coarse.  INDICES is an Nx2 matrix containing the sizes of
% each subband.  This is compatible with the MatLab Wavelet toolbox.

% Eero Simoncelli, 6/96.

function [pyr,pind] = buildLpyr(im, ht, filt1, filt2, edges)

if (nargin < 1)
  error('First argument (IM) is required');
end

im_sz = size(im);

%------------------------------------------------------------
%% OPTIONAL ARGS:

if (exist('filt1') ~= 1)
  filt1 = 'binom5';
end

if isstr(filt1)
  filt1 = namedFilter(filt1);
end

if ( (size(filt1,1) > 1) & (size(filt1,2) > 1) )
  error('FILT1 should be a 1D filter (i.e., a vector)');
else
  filt1 = filt1(:);
end

if (exist('filt2') ~= 1)
  filt2 = filt1;
end

if isstr(filt2)
  filt2 = namedFilter(filt2);
end

if ( (size(filt2,1) > 1) & (size(filt2,2) > 1) )
  error('FILT2 should be a 1D filter (i.e., a vector)');
else
  filt2 = filt2(:);
end

max_ht = 1 + maxPyrHt(im_sz, max(size(filt1,1), size(filt2,1)));
if (exist('ht') ~= 1)
  ht = max_ht;
else
  if (ht > max_ht)
    error(sprintf('Cannot build pyramid higher than %d levels.',max_ht));
  end
end

if (exist('edges') ~= 1)
  edges= 'reflect1';
end

%------------------------------------------------------------

if (ht <= 1)

  pyr = im(:);
  pind = im_sz;

else

  if (im_sz(2) == 1)
    lo2 = corrDn(im, filt1, edges, [2 1], [1 1]);
  elseif (im_sz(1) == 1)
    lo2 = corrDn(im, filt1', edges, [1 2], [1 1]);
  else
    lo = corrDn(im, filt1', edges, [1 2], [1 1]);
    int_sz = size(lo);
    lo2 = corrDn(lo, filt1, edges, [2 1], [1 1]);
  end

  [npyr,nind] = buildLpyr(lo2, ht-1, filt1, filt2, edges);

  if (im_sz(1) == 1)
    hi2 = upConv(lo2, filt2', edges, [1 2], [1 1], im_sz);
  elseif (im_sz(2) == 1)
    hi2 = upConv(lo2, filt2, edges, [2 1], [1 1], im_sz);
  else
    hi = upConv(lo2, filt2, edges, [2 1], [1 1], int_sz);
    hi2 = upConv(hi, filt2', edges, [1 2], [1 1], im_sz);
  end

  hi2 = im - hi2;

  pyr = [hi2(:); npyr];
  pind = [im_sz; nind];

end


text:
function par_clearworkers(endpoints,nokillall)
% Clear workers from the listed hostnames (or if those are unspecified, from the current pool).
% par_clearworkers(Hostnames,NoKillall)
%
% Note: this function erases all of the user's MATLAB processes on the given nodes; if there are
% unrelated MATLABs running, these will also be killed. This may also shoot down the master process
% itself if it is located on one of the unfortunate machines...
%
% In:
%   Hostnames : cell array of hostnames or 'pid@host' or 'pid@host:port' strings on which to kill
%               workers
%
%   NoKillall : if given, this function will not perform "killall" for hosts where no process id
%               is given (default: false)
%
% Example:
%   % clear all workers on the given machines (of this user)
%   par_clearworkers({'computer1','computer2','computer3'})
%
%   % clear only specific workers with given pid's, and all MATLAB instances on computer3
%   par_clearworkers({'31443@computer1','41443@computer2','computer3'})
%
%   % same as before, but ignore computer3 (no killall)
%   par_clearworkers({'31443@computer1','41443@computer2','computer3'},true)
%
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2011-02-15

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

% optionally look up endpoints from the global pool
if ~exist('endpoints','var') || isempty(endpoints)
    endpoints = par_globalsetting('pool'); end
if ~iscellstr(endpoints)
    error('The given Hostnames argument must be a cell array of strings.'); end
if ~exist('nokillall','var') || isempty(nokillall)
    nokillall = false; end

% drop any port assignments
for i=1:length(endpoints)
    colons = endpoints{i}==':';
    if any(colons)
        hostpid{i} = endpoints{i}(1:find(colons)-1);
    else
        hostpid{i} = endpoints{i};
    end
end

% also split off the pid parts
for i=1:length(hostpid)
    match = hostpid{i}=='@';
    if any(match)
        hosts{i} = hostpid{i}(find(match,1)+1:end);
        pids{i} = hostpid{i}(1:find(match,1)-1);
    else
        hosts{i} = hostpid{i};
        pids{i} = 'all';
    end
end

% remove duplicates
uniquehosts = unique(hosts);
pids_per_host = repmat({{}},1,length(uniquehosts));
killall = false(1,length(uniquehosts));

% for each unique host, accumulate the pids to kill (or record 'all' if no pid given)
for i=1:length(hosts)
    % position of this host in uniquehosts, pids_per_host, and killall
    pos = find(strcmp(uniquehosts,hosts{i}));
    if strcmp(pids{i},'all')
        killall(pos) = true;
    else
        pids_per_host{pos}{end+1} = pids{i};
    end
end

% if the own host is among them, move it last
myhost = hlp_hostname;
match = strncmp(myhost,uniquehosts, min(cellfun('length',[{myhost} uniquehosts])));
if any(match)
    uniquehosts = [uniquehosts(~match) uniquehosts(match)]; end

% now clear them
fprintf('Clearing workers on %i machines...',length(uniquehosts));
for h = 1:length(uniquehosts)
    if ~isempty(pids_per_host{h})
        [status,info] = system(sprintf('ssh -x %s kill %s',uniquehosts{h},sprintf('%s ',pids_per_host{h}{:}))); end
    if ~nokillall && killall(h)
        [status,info] = system(sprintf('ssh -x %s killall MATLAB',uniquehosts{h})); end
end %#ok<NASGU,ASGLU>
fprintf('done.\n');

repo: REST


text:
%--- help for data_prerequest ---
%
%  INTERNAL FUNCTION
%
%

author: jmaih


text:
function pde = Maxwellsaddledata
%% MAXWELLDATA2 non-homogenous Dirichlet/Neumann boundary condition
%
%                   curl(mu^(-1)curl u)      = J    in \Omega,
%                                  div u     = 0    in \Omega.
%                                 n \cross u = n � g_D  on \Gamma_D,
%                   n \cross (mu^(-1)curl u) = n � g_N  on \Gamma_N.
%
%   mu = 1;
%   u = [0 cos(x) cos(x)];
%   curlu = [0 sin(x) -sin(x)];
%   J   = u;
%   g_N = n \cross curl u;
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

pde.mu = 1;
pde.J = @J;
pde.exactu = @exactu;
pde.g_D = @exactu;
pde.curlu = @curlu;
pde.g_N = @g_N;

    function s = J(p)
        s = exactu(p);
    end

    function s = exactu(p)
        x = p(:,1); %y = p(:,2); z = p(:,3);
        s = [0*x, cos(x), cos(x)];
    end

    function s = curlu(p)
        x = p(:,1); %y = p(:,2); z = p(:,3);
        s = [0*x, sin(x), -sin(x)];
    end

    function s = g_N(p,n)
        s = curlu(p);
        s = cross(n,s,2);
    end
end

repo: ifem


text:
%
%  Context-Aware Correlation Filters
%
%  Written by Luca Bertinetto, 2016
%  Adapted by Matthias Mueller, 2016
%
%  This function takes care of setting up parameters, loading video
%  information and computing precisions. For the actual tracking code,
%  check out the TRACKERMAIN.m function.
%

function run_tracker(video, start_frame)
% RUN_TRACKER  is the external function of the tracker - does initialization and calls trackerMain

	%path to the videos (you'll be able to choose one with the GUI).
	base_path = 'sequences/';

	%default settings
	if nargin < 1, video = 'Skiing'; end
	if nargin < 2, start_frame = 1; end
    if nargin < 3, show_plots = 1; end

    %% Read params.txt
    params = readParams('params.txt');
	%% load video info
    sequence_path = fullfile(base_path,video);
    img_path = fullfile(sequence_path, 'img');
    %% Read files
    text_files = dir([sequence_path '*_frames.txt']);
    if(~isempty(text_files))
        f = fopen([sequence_path text_files(1).name]);
        frames = textscan(f, '%f,%f');
        fclose(f);
    else
        frames = {};
    end
    if exist('start_frame')
        frames{1} = start_frame;
    else
        frames{1} = 1;
    end

    params.bb_VOT = csvread(fullfile(sequence_path, 'groundtruth_rect.txt'));
    region = params.bb_VOT(frames{1},:);
    %%%%%%%%%%%%%%%%%%%%%%%%%
    % read all the frames in the 'imgs' subfolder
    dir_content = dir(fullfile(sequence_path, 'img'));
    % skip '.' and '..' from the count
    n_imgs = length(dir_content) - 2;
    img_files = cell(n_imgs, 1);
    for ii = 1:n_imgs
        img_files{ii} = dir_content(ii+2).name;
    end

    img_files(1:start_frame-1)=[];

    im = imread(fullfile(img_path, img_files{1}));
    % is a grayscale sequence ?
    if(size(im,3)==1)
        params.grayscale_sequence = true;
    end

    params.img_files = img_files;
    params.img_path = img_path;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if(numel(region)==8)
        % polygon format
        [cx, cy, w, h] = getAxisAlignedBB(region);
    else
        x = region(1);
        y = region(2);
        w = region(3);
        h = region(4);
        cx = x+w/2;
        cy = y+h/2;
    end

    % init_pos is the centre of the initial bounding box
    params.init_pos = [cy cx];
    params.target_sz = round([h w]);
    [params, bg_area, fg_area, area_resize_factor] = initializeAllAreas(im, params);
	if params.visualization
		params.videoPlayer = vision.VideoPlayer('Position', [100 100 [size(im,2), size(im,1)]+30]);
	end
    % in runTracker we do not output anything
	params.fout = -1;
	% start the actual tracking
	results = trackerMain(params, im, bg_area, fg_area, area_resize_factor);

    %calculate and show precision plot, as well as frames-per-second
    precisions = precision_plot(results.res, params.bb_VOT, video, show_plots);
    fprintf('%12s - Precision (20px):% 1.3f, FPS:% 4.2f\n', video, precisions(20), results.fps)
    fclose('all');

end

repo: Context-Aware-CF-Tracking


text:
function [d,dt]=formatts(d)
%
% Usage: [d,dt]=formatts(d)
%
% Helper function for CWT,XWT,WTC
%
% Brings a timeseries into a shape so that it has two columns: [time, value].
%
%
% (C) Aslak Grinsted 2002-2004
%

% -------------------------------------------------------------------------
%   Copyright (C) 2002-2004, Aslak Grinsted
%   This software may be used, copied, or redistributed as long as it is not
%   sold and this copyright notice is reproduced on each copy made.  This
%   routine is provided as is without any express or implied warranties
%   whatsoever.

if (ndims(d)>2)
    error('Input time series should be 2 dimensional.');
end
if (numel(d)==length(d))
    d=[(1:length(d))' d(:)];
end
if size(d,1)<size(d,2)
    d=d';
end
if (size(d,2)~=2)
    error('Time series must have 2 columns.')
end

if (d(2,1)-d(1,1)<0)
    d=flipud(d);
end

dt=diff(d(:,1));
if any(abs(dt-dt(1))>1e-1*dt(1))
    error('Time step must be constant.');
end
if (dt==0)
    error('Time step must be greater than zero.')
end

dt=dt(1);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function status = unitTest_BandpassFilt_LPF(datafmt, dirname, newval)
global procStreamStyle
global testidx
global logger

if isempty(procStreamStyle)
    procStreamStyle = datafmt;
end
if includes(procStreamStyle,'snirf')
    datafmt = 'snirf';
end
if isempty(testidx)
    testidx = 0;
end
testidx = testidx+1;

status = -1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inititialize and error check input argument
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~exist('datafmt','var')
    datafmt = 'nirs';
end
if ~exist('dirname','var')
    return;
end
if ~exist('newval','var')
    newval = [];
end
logger = InitLogger(logger, 'unitTest_BandpassFilt_LPF');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set up logger and other administrative paramaters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
logger.Write('################################################################');
logger.Write(sprintf('Running test #%d - unitTest_BandpassFilt_LPF(''%s'', ''%s'', %0.1f)\n', testidx, datafmt, dirname, newval));
fprintf('\n');

[~, currpath] = findRootFolder(dirname);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Change processing stream param values to newval
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[dataTree, procStreamConfigFile] = changeProcStream(datafmt, 'processOpt_default_homer3', 'hmrR_BandpassFilt_old', 'lpf', newval);
if isempty(dataTree)
    status = exitEarly(sprintf('#%d - unitTest_BandpassFilt_LPF(''%s'', ''%s'', %0.1f): SKIPPING - This test does not apply to %s.\n', ...
                               testidx, datafmt, dirname, newval, dirname), logger);
    return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate Homer3 output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
logger.Write(sprintf('Loaded processing stream from %s\n', procStreamConfigFile));
generateOutput(dataTree);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compare Homer3 output to the available Homer2 output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
groupResults_h2 = mydir('./groupResults_homer2_*.mat');
s = zeros(length(groupResults_h2), 2);
for iG=1:length(groupResults_h2)
    group_h2 = load(groupResults_h2(iG).name);
    [~, groupResults_h2(iG).rootdir] = fileparts(groupResults_h2(iG).rootdir);

    s(iG,1) = compareDcAvg(group_h2, dataTree, 'dcAvg');
    s(iG,2) = compareProcStreams(dataTree, groupResults_h2(iG));

    msgs = MatchMessages(sum(s(iG,:)));
    logger.Write(sprintf('Comparing output to %s  ==>  Outputs: %s,   Proc Streams: %s\n', groupResults_h2(iG).name, msgs{1}, msgs{2}));
end
iMatch = find(s(:,1)==0 & s(:,2)==0);
if ~isempty(iMatch)
    status=0;
else
    status=1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Report results of the comparison
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
logger.Write(sprintf('\n'));
if status==0
    logger.Write(sprintf('#%d - unitTest_BandpassFilt_LPF(''%s'', ''%s'', %0.1f): TEST PASSED - Homer3 output matches %s.\n', ...
             testidx, datafmt, dirname, newval, [groupResults_h2(iMatch).rootdir, '/', groupResults_h2(iMatch).name]));
else
    logger.Write(sprintf('#%d - unitTest_BandpassFilt_LPF(''%s'', ''%s'', %0.1f): TEST FAILED - Homer3 output does NOT match ANY Homer2 groupResults.\n', testidx, datafmt, dirname, newval));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Clean up before exiting
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
logger.Write('\n');
logger.Close('unitTest_BandpassFilt_LPF');

dataTree.ResetAllGroups();

cd(currpath);

author: BUNPC


text:
% test for optimization of the conditionning of the CS matrix

n = 100; m = 300;

[M,eOrth,eAngle,eNorm] = perform_cs_matrix_optimization(randn(n,m), options);

clf;
subplot(2,1,1);
plot(log(eOrth(3:end))); axis tight;
title('Orth');
subplot(2,1,2);
plot(log(eNorm(3:end))); axis tight;
title('Norm');

return;

%%% OLD CODE %%%

c = 0.15*2;
c = Inf;

niter = 20;
eOrth = [];
eNorm = [];
eAngle = [];
for i=1:niter
    % orthogonality constraint
    eOrth(end+1) = norm( M*M' - eye(n)*m/n );
    M = orth(M')';
    M = M(1:n,:) * sqrt(m/n);
    % norm constraint
    d = sqrt( sum( M.^2, 1 ) );
    eNorm(end+1) = sum( abs(d-1) );
    M = M./repmat( d, [n 1] );
    % angle constraint
    if 1
    G = M'*M;
    eAngle(end+1) = norm( G-eye(m) );
    G1 = min(max(G,-c),c);
    G1 = G1-diag(diag(G1))+diag(ones(1,m)); % remet ? 1 la diagonale
    [U,S,V] = svd(G1);
    M = U*diag(sqrt(diag(S)));
    M = U(:,1:n)';
    end
end

clf;
subplot(2,1,1);
plot(log(eNorm(3:end)));
subplot(2,1,2);
plot(log(eOrth(3:end)));

author: gpeyre

repo: matlab-toolboxes


text:
%PLOTSOM Plot the Self-Organizing Map in 2D
%
%    PRPLOTSOM(W)
%
% Plot the Self-Organizing Map W, trained by som.m. This is only
% possible if the map is 2D.
%
% SEE ALSO (PRTools Guide)
% SOM

% Copyright: D.M.J. Tax, davidt@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% Maybe I should introduce the possibility to set the linewidth and
% markersize...
%
% Changes:
% DR1 - Dick de Ridder, 02-03-2006:
%       If the mapped is trained on 3D data, plot it like this.

function h = prplotsom(W)

if ~ismapping(W)
	error('SOM mapping expected');
end
if ( ~strcmp(getmapping_file(W),'som') & ...
	 ~strcmp(getmapping_file(W),'som_dd') )
    error('I expect a SOM mapping!');
end
if size(W.data.neurons,2)~=2
	error('The SOM can only be plotted in 2D');
end

% Get the data:
W = +W;  w=W.neurons; k=W.k;
% Plot the bloody thing:
hold on;

% DR: Handle maps trained on 3D data.

if (size(w,2) == 3)
	% The 'horizontal' lines:
	for i=0:k(2)-1
			h=plot3(w(i*k(1)+(1:k(1)),1),w(i*k(1)+(1:k(1)),2),w(i*k(1)+(1:k(1)),3),'o-');
			set(h,'linewidth',2,'markersize',8);
	end
	I = reshape(1:k(1)*k(2),k(1),k(2))';
	% The 'vertical' lines:
	for i=0:k(1)-1
			h=plot3(w(I(i*k(2)+(1:k(2))),1),w(I(i*k(2)+(1:k(2))),2),w(I(i*k(2)+(1:k(2))),3),'o-');
			set(h,'linewidth',2,'markersize',8);
	end
	view(3);
else
	% The 'horizontal' lines:
	for i=0:k(2)-1
			h=plot(w(i*k(1)+(1:k(1)),1),w(i*k(1)+(1:k(1)),2),'o-');
			set(h,'linewidth',2,'markersize',8);
	end
	I = reshape(1:k(1)*k(2),k(1),k(2))';
	% The 'vertical' lines:
	for i=0:k(1)-1
			h=plot(w(I(i*k(2)+(1:k(2))),1),w(I(i*k(2)+(1:k(2))),2),'o-');
			set(h,'linewidth',2,'markersize',8);
	end
end;

% Return the handle only when it is required
if (nargout==0)
	clear h;
end

return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [F_maxSize, F_maxLoc] = getFeatMatHeatmap( jobID, taskID, objIDs, CONTIG_THR )
% Identify features covering a contig block w/ length at least CONTIG_THR
% Input:
%    jobID : integer id of job
%    taskID : integer id of which run of the job to examine
%    objIDs : vector of objects to consider
%               if empty, we look at all objects
%    CONTIG_THR : integer specifying minimum length of a contig block
%                  mocap data recorded at 120fps,
%                  block-averaged so more like 10 fps
%                  thus setting CONTIG_THR=10 means 1 sec
% OUTPUT
%   F_maxSize : N x K matrix
%                  where entry n,k gives integer length of *longest*
%                   contig block assigned to feature k in sequence n
%   F_maxLoc  : N x K matrix
%                  where entry n,k gives timestep id of the start of the
%                    where *longest* contig block of feat k in seq. n

if ~exist( 'CONTIG_THR', 'var' )
    CONTIG_THR = 10;
end

INFO = loadSamplerInfo( jobID, taskID );
Q = loadSamplerOutput( jobID, taskID );

data = INFO.data;

F = Q.Psi(end).F;
stateSeq = Q.Psi(end).stateSeq;

K = size(F,2);
if ~exist( 'objIDs', 'var') || isempty( objIDs )
    objIDs = 1:size(F,1);
end
N = length(objIDs);

F_maxSize = zeros( N, K );
F_maxLoc = zeros( N, K );
for aa = 1:length(objIDs )
    [bSiz,bLoc,bVal] = getContigBlocks( stateSeq( objIDs(aa) ).z );

    us = unique( bVal );
    for uu = 1:length( us )
        occurIDs = find( bVal == us(uu) );

        keepers =  bSiz( occurIDs ) > CONTIG_THR;
        if any( keepers )
            occurIDs = occurIDs(keepers);
            if isempty( occurIDs )
                continue;
            end
            [maxSize,maxID] = max( bSiz( occurIDs ) );
            F_maxSize( aa, us(uu) ) = maxSize;
            F_maxLoc( aa, us(uu) ) = bLoc( occurIDs(maxID) );
        end

    end

end

%           if strcmp( INFO.model.obsModel.type, 'Multinomial' )
%               % Make sure each interval has STIPS in >1/2 of the included tsteps
%               doKeep = true(1, length(occurIDs) );
%               for oo = 1:length( occurIDs )
%                   tstart = bLoc(occurIDs(oo));
%                   tstop  = tstart+bSiz( occurIDs(oo) )-1;
%                   if sum( data_struct(aa).nEmissions(tstart:tstop) ) < length( tstart:tstop )/2
%                       doKeep(oo) = false;
%                   end
%               end
%               occurIDs = occurIDs( doKeep );
%           end

author: michaelchughes

repo: NPBayesHMM


text:
function varargout = solvemp(P,x,y,options)
%SOLVEMP Computes solution to multi-parametric optimization problem
%
% min_z(x)   h(x,z)
% subject to
%            F(x,z) > 0
%
%
% [SOL, DIAGNOSTIC,Z,HPWF,ZPWF] = SOLVEMP(P,options,x,y)
%
% SOL        : Multi-parametric solution (see MPT toolbox)
%
% DIAGNOSTIC : struct with diagnostic information
%
% Z          : SDPVAR object with the detected decision variable z
%
% HPWF       : The value function as a pwf function
%
% ZPWF       : The optimal decision variable as a pfw function
%
% Input
%    P        : Optimization model
%    options  : solver options. See SDPSETTINGS. Can be [].
%    x        : Parametric variables
%    y        : Requested decision variables (subset of z)
%
% NOTE : If you are solving a problem leading to an mpMILP, the
% output SOL will be a set-valued map. To obtain the minimal
% solution (without so called overlaps), run removeOverlaps(SOL). If you
% have requested the 5th output ZPWF, overlaps are automatically removed.
% If your problem leads to an mpMIQP,  the output SOL will also be a
% set-valued map, but there is currently no way in MPT to obtain a
% non-overlapping solution. To use the solution in MPT, the command
% mpt_mergeCS(SOL) can be useful. Notice that the fifth output argument
% not will be available for mpMIQP problems.
%
% See also PARAMETRIC, SDPSETTINGS, YALMIPERROR

if nargin < 4
    options = sdpsettings;
else
    options = P.Options;
end

[sol,info] = solvemp(P.Constraints,P.Objective,options,x,y)

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function varargout = chol(varargin)
%CHOL   Cholesky factorization of a DISKFUN.
%
%   R = CHOL( F ), if F is a nonnegative definite DISKFUN then this returns
%   an upper triangular quasimatrix so that R'*R is a decomposition of F.
%   If F is not nonnegative definite then an error is thrown.
%
%   L = CHOL(F, 'lower'), if F is a nonnegative definite DISKFUN then this
%   produces a lower triangular quasimatrix so that L*L' is a decomposition
%   of F. If F is not nonnegative definite then an error is thrown.
%
%   [R, p] = CHOL( F ), with two outputs never throwns an error message. If
%   F is nonnegative definite then p is 0 and R is the same as above. If F
%   is symmetric but negative definite or semidefinite then p is a positive
%   integer such that R has p columns and R'*R is a rank p nonnegative
%   definite DISKFUN that approximates F. This is particular useful when F
%   is nonnegative definite, but rounding error have perturbed it to be
%   semidefinite.
%
%   [L, p] = CHOL(F, 'lower') same as above but the first argument is lower
%   triangular.
%
%   For more information about the factorization:
%   A. Townsend and L. N. Trefethen, Continuous analogues of matrix
%   factorizations, Proc. Royal Soc. A., 2015.
%
% See also DISKFUN/LU, and DISKFUN/QR.

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

% Requires evaluations of f in polar coords once passed to separableApprox,
% so we make a chebfun2:
f = cart2pol(varargin{1}, 'cdr');

[varargout{1:nargout}] = chol@separableApprox(f, varargin{2:end});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function prime_test ( )

%*****************************************************************************80
%
%% PRIME_TEST tests PRIME.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 December 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PRIME_TEST\n' );
  fprintf ( 1, '  PRIME returns primes from a table.\n' );

  n = -1;
  prime_max = prime ( n );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of primes stored is %d\n', prime_max );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I    Prime(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 10
    fprintf ( 1, '  %4d  %6d\n', i, prime(i) );
  end
  fprintf ( 1, '\n' );
  for i = prime_max - 10 : prime_max
    fprintf ( 1, '  %4d  %6d\n', i, prime(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m