text:
function [varargout] = xpower(x,varargin)
% [z80_ind,z80_grp,OUT] = xpower(x,[c],[true],[N],[gsd],[V],[S],[trueX])
% Calculates individual and group power as a function of design matrix,
% true effect magnitudes, and contrasts
% using central limit theorem approximations/theory for estimates
%
% tor wager, 2 / 22 / 04
%
% Defaults
% gsd = .28;          % group standard deviation, relative to true effect;
%                       default from VNL experiment (could be overestimate, bec. btwn + within
%                       not separated)
% N = 10;             % number of subjects, for group stats
% true                true effects; 1 is an increase of 1 std. dev. of scanner noise per event
%                     default is .47 SNR everywhere the first contrast has positive
%                     weights
%                     VNL estimates of signal / noise ratio for visual impulse response, 16 Hz
%                     flashing checkerboard for .5 s
%                     0.83% signal change response, 1.67% signal change error std. dev., .47 SNR
% c = [1 1 -1 -1];    % contrast matrix; one row per contrast, one column
% per condition; default is one contrast for each effect
%
% V                   is noise autocorrelation matrix,
%                       symmetric and ones on the diagonal, n x n (n is number of samples)
%                       default is white noise; see also
%                       canonical_autocorrelation
%
% S                   is filtering (high/low pass) matrix; see getSmoothing
%
%
% see power_map1.m for examples and scripts
% see xzpower for a simulation-based approach
%
% outputs:
% z80_ind and z80_grp are the power estimates for individual and group
% designs.  These correspond to the z-value that is expected to be reached
% in 80% of the realizations of the design.

% Another way of creating the matrices:
%E = inline('inv(x'' * x) * x'' * Vi','x','Vi'); % E is error cov matrix
% matrix of orthogonal projection onto model space * Vi
% test: demo to verify that E is constant for constant SNR, regardless of
% sigma and t (as long as their ratio is equivalent)
% sigma = 1; t = .5; e = E(S*X,S* sigma * Vi);ee = e * e', t ./ trace(ee).^.5

% -----------------------------------------------------------------
% Set up input arguments and defaults
% -----------------------------------------------------------------

% notes: 2/7/07:
% keep true signal |h| = 1, or input as entered
% wsd and gsd are error st. deviations, so SNR = true / wsd (individual
% level) and true / ...

wsd = .47;          % within-subjects d, signal / noise == wsd / 1
gsd = .28;            % group standard deviation, relative to true effect; default from VNL
N = 10;             % number of subjects, for group stats

for i = 1:length(varargin)
    switch i
        case 1
            if isempty(varargin{i})
                c = [eye(size(x,2)-1) zeros(size(x,2)-1,1)];
                % contrast matrix; one row per contrast, one column per condition
            else
                c = varargin{i};
            end
        case 2
            if isempty(varargin{i})
                true = c(1,:) > 0;   % true effects; 1 is an increase of magnitude 1; SNR depends on wsd and gsd
                % % % true = true .* .47;  % estimated SNR from visual experiment (VNL), n = 11
            else
                true = varargin{i};
            end
        case 3
            N = varargin{i};
        case 4
            gsd = varargin{i};
        case 5
            V = varargin{i};
        case 6
            S = varargin{i};
        case 7
            trueX = varargin{i};
    end
end

% defaults, if not entered

if ~(exist('c') == 1) || isempty(c), c = [eye(size(x,2)-1) zeros(size(x,2)-1,1)];   end
if ~(exist('true') == 1) || isempty(true),true = 1 .* (c(1,:) > 0); end

if ~(exist('V') == 1) || isempty(V), V = eye(size(x,1)); xc = [1 0 0];
    %[xc,V] = canonical_autocorrelation
end
if ~(exist('S') == 1), S = []; end

if ~(exist('trueX') == 1) || isempty(trueX), trueX = x; end

% defaults, if empty
%if isempty(true),true = 0.47 .* (c(1,:) > 0);,end
if isempty(N),N = 10; end
if isempty(gsd),gsd = 0.28; end
%if isempty(V), V = eye(size(x,1)); ,xc = [1 0 0];,end

% -----------------------------------------------------------------
% Do the power calculation - start with noncentrality parameters
% -----------------------------------------------------------------

if ~isempty(S),
    if size(S,1) ~= size(S,2), error('S is not square.'), end
    if size(S,1) ~= size(x,1), error('S and X dims do not match.'), end
    if size(S,1) ~= size(V,1), error('S and Vi dims do not match.'), end

    %xo = x;     % save original x, for true signal estimate
    x = S * x; V = S * V * S';  %V = S * V;
end

if size(c,2) < size(x,2), c(:,end+1) = 0;, end
if size(true,2) < size(x,2), true(:,end+1) = 0;, end

c = c';

if ~isempty(S)
    %y = S * x * true'; % S * xorig * true and x * true are equivalent; no need for xo
    % tor modified may 06
    y = S * trueX * true';      % which = S * (xorig*true), filtered orig. data
else
    y = trueX * true';
end

xtxi = inv(x' * x);
px = xtxi * x'; %pinv(x);

%b = c' * xtxi * x' * y; % REDUCES to input vector true, regardless of filtering; b = true';
b = c' * px * y;

% i think this one below has V as the sqrt of the V in Friston
%sx = sqrt(diag(c' * xtxi * x' * V' * V * x * xtxi * c))   % design part -- std of x, sqrt(var(x))
sx = sqrt(diag(c' * px * V * px' * c));

% %n = b ./ sx;    % noncentrality
% mod 2/07
n = b ./ (wsd .* sx);

% GETTING THE NONCENTRALITY PARAMETER n
% noncentrality parameter n = snr * b/se
% snr = |h| / sigma, where |h| = scaling of impulse response to neural
% event.  see Zarahn 2001, "Reference Effect for Power Analysis"
% We can simplify the computation by always assuming |h| = 1 and varying
% sigma, in the simplest case setting snr = 1.
%
% t-values are directly related to snr and n.
% t = (|h| * n) / sigma
% in the simplest case, where snr = 1, t = n.
%
% n is contrast weight invariant.  contrasts can be normalized by sum of
% squared contrast weights.
%
% FROM n to POWER
% By the central limit theorem, n = t if snr = 1, and n is normally
% distributed with std 1.  Then we can integrate the normal distribution.
%
% but since the noncentrality param is essentially a t-value, as is the 20%
% level for the 80% power value, we have to do a small additional
% adjustment: convert t to Z scores based on the df.

% Degrees of freedom: Worsley and Friston (1995) way (effective df),
% df reduction due to autocorrelation.  This will be different if
% a prewhitening strategy is used.

R = eye(size(x,1)) - x * xtxi * x';
U = R * V * V';
df = (trace(U) .^ 2) ./ (trace(U * U));

%df = size(x,1) - size(x,2);  % standard way, but not correct with
%autocorrelated errors

%t80_ind = norminv(.2,n,1);
% mod 2/07
t80_ind = nctinv(.2,df,n);  % 20th% of noncentral t distribution
ztmp = spm_t2z(t80_ind,df);
ztmp = min(ztmp,t80_ind);    % use min, as conversion fails for high-power

% % tmp = norminv(tcdf(t80_ind,df));
% % tmp = min(tmp,t80_ind);             % use min, as conversion fails for high-power
% %
varargout{1} = ztmp;

% GROUP POWER
% group power is calculated by considering the contribution of the
% individual se to the group se, and assuming a fixed intersubject variance (set nominally to 1).  with a
% with group rrue intersubject standard deviation (gsd) of 1, then
% gse = (gsd + E(indiv_se)) ./ sqrt(2N)
% by the central limit theorem.
% we use sqrt(2N) because std dev of n rnd variables = (s1 + s2 ... sn) / sqrt(n)
% individual se (ise) is the standard deviation of the sampling distribution
% gse is the standard deviation
%
% take as ise the coefficient of variation (cv), sx / b, which is contrast
% scale independent.

if nargout > 1

    %gse = (gsd + (sx ./ b)) ./ sqrt(2 * N);

    %ng = 1 ./ gse;  % normalized b / gse.  this is contrast weight-scaling independent

    % group standard error.  sqrt of sum of within + between, over
    % sqrt(N-1) to get standard error from standard deviation
    gse = sqrt( (wsd.*sx).^2 + gsd.^2 ) ./ sqrt(N - 1);

    ng = true(1:end-1)' ./ gse;  %  b / gse.  this is contrast weight-scaling independent

    dfg = N - 1;
    t80_grp = nctinv(.2, dfg, ng);  % 20th% of noncentral t distribution
    ztmp = spm_t2z(t80_grp,dfg);
    ztmp = min(ztmp,t80_grp);    % use min, as conversion fails for high-power

% %     t80_grp = norminv(.2,ng,1);
% %
% %     tmp = norminv(tcdf(t80_grp,dfg));   % adjust to Z-score based on df
% %     tmp = min(tmp,t80_grp);             % use min, as conversion fails for high-power
    varargout{2} = ztmp;

    if nargout > 2

        OUT.x = x;
        OUT.gsd = gsd;
        OUT.N = N;
        OUT.true = true;
        OUT.c = c';
        OUT.y = y;
        OUT.b = b;
        OUT.sx = sx;
        OUT.n = n;
        OUT.df = df;
        OUT.ng = ng;

        OUT.t80_ind = t80_ind;
        OUT.t80_grp = t80_grp;
        OUT.expected_t = b ./ (wsd .* sx);

        varargout{3} = OUT;

    end

end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function out = nsamdf_nonlinear_measure(data,fs,winLenRel,shiftLenRel,lagRel,degree,doPlot)
% nsamdf_nonlinear_measure computes the nonlinearity measure L through nsAMDF
% (nonlinear average magnitude difference function), developed by:
%
% Ozkurt et al. (2020), "Identification of nonlinear features in cortical and
% subcortical signals of Parkinson's Disease patients via a novel efficient
% measure", Neuroimage.
%
% Please refer to and cite this paper if you utilize this function in your work.
%
%---INPUTS:
%
% data = One dimensional input time-series (it can be raw, but it is a good idea to low-pass filter it to get rid of high frequency nuisance)
%        For the neural data (LFP, MEG) in the aferomentioned paper, we
%        low-passed the raw data for 40 Hz. The data should be long enough for proper estimation of nonlinearity.
%
% winlen = window length (a long enough segment is important to estimate the nonlinearity)
%          We chose the window length as 14 sec., i.e., winlen = 14*fs
%
% shiftlen = This amounts to window length - overlap length btw windows
%            We chose it as half the window length, i.e., shiftlen = 0.5*winlen
%
% lag = TMaximum lag for nsAMDF, we chose it as 1 sec, i.e., lag = fs.
%
% fs = sampling frequency of the data
%
% degree = The chosen degree p should ideally be large enough to capture the highest order of nonlinearity within the data.
%          We chose p=7 for in our case of Parkinsonian data in the paper.
%
% doPlot = true to plot nsAMDF sequences, otherwise just assign it false
%
%---OUTPUTS:
%
% L: nonlinearity measure
%
% s2: normalized nsAMDF for the degree 2
%
% sd: normalized nsAMDF for the chosen degree greater than 2
%
%
% Required subfunctions are NormedSingleCurveLengthWindowed.m & NormedSingleCurveLength.m
%
%
%   Authored by Tolga Esat Ozkurt, 2020. (tolgaozkurt@gmail.com)

%-------------------------------------------------------------------------------
% Set defaults:

if nargin < 2
    fs = 1;
end
if nargin < 3
    winLenRel = 14;
end
windowLength = winLenRel*fs;
if nargin < 4
    shiftLenRel = 0.5;
end
shiftLength = shiftLenRel*windowLength;
if nargin < 5
    lagRel = 1;
end
lag = fs*lagRel;
if nargin < 6
    degree = 7;
end
if nargin < 7
    doPlot = false
end

%-------------------------------------------------------------------------------
% nsAMDF for p=2:
s2 = NormedSingleCurveLengthWindowed(data,windowLength,shiftLength,lag,fs,2);
out.s2 = s2 ./ max(s2); % normalized

% nsAMDF for p=degree:
sd = NormedSingleCurveLengthWindowed(data,windowLength,shiftLength,lag,fs,degree);
out.sd = sd ./ max(sd); % normalized

% If you like, you can bandpass filter s2 and sd for the specific frequency band
% of nonlinear effect both to compute L and plot them as such

out.L = norm(s2 - sd);

%-------------------------------------------------------------------------------
if doPlot
    figure
    plot(s2,'b')
    hold on
    plot(sd,'g')
end

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function r8row_variance_test ( )

%*****************************************************************************80
%
%% R8ROW_VARIANCE_TEST tests R8ROW_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  m = 3;
  n = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8ROW_VARIANCE_TEST\n' );
  fprintf ( 1, '  For a R8ROW (a matrix regarded as rows):\n' );
  fprintf ( 1, '  R8ROW_VARIANCE computes variances;\n' );

  k = 0;
  for i = 1 : m
    for j = 1 : n
      k = k + 1;
      a(i,j) = k;
    end
  end

  r8mat_print ( m, n, a, '  The original matrix:' );

  variance = r8row_variance ( m, n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Row variances:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : m
    fprintf ( 1, '  %3d  %10f\n', i, variance(i) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear, clc;

% This is an example for running the function tree_LogisticR
%
%  Problem:
%
%  min  f(x,c) = - sum_i weight_i * log (p_i) + z * sum_j w_j ||x_{G_j}||
%
%  a_i denotes a training sample,
%      and a_i' corresponds to the i-th row of the data matrix A
%
%  y_i (either 1 or -1) is the response
%
%  p_i= 1/ (1+ exp(-y_i (x' * a_i + c) ) ) denotes the probability
%
%  G_j's are nodes with tree structure
%
%  The tree structured group information is contained in
%  opts.ind, which is a 3 x nodes matrix, where nodes denotes the number of
%  nodes of the tree.
%
%  opts.ind(1,:) contains the starting index
%  opts.ind(2,:) contains the ending index
%  opts.ind(3,:) contains the corresponding weight (w_j)
%
%  Note:
%  1) If each element of x is a leaf node of the tree and the weight for
%  this leaf node are the same, we provide an alternative "efficient" input
%  for this kind of node, by creating a "super node" with
%  opts.ind(1,1)=-1; opts.ind(2,1)=-1; and opts.ind(3,1)=the common weight.
%
%  2) If the features are well ordered in that, the features of the left
%  tree is always less than those of the right tree, opts.ind(1,:) and
%  opts.ind(2,:) contain the "real" starting and ending indices. That is to
%  say, x( opts.ind(1,j):opts.ind(2,j) ) denotes x_{G_j}. In this case,
%  the entries in opts.ind(1:2,:) are within 1 and n.
%
%
%  If the features are not well ordered, please use the input opts.G for
%  specifying the index so that
%   x( opts.G ( opts.ind(1,j):opts.ind(2,j) ) ) denotes x_{G_j}.
%  In this case, the entries of opts.G are within 1 and n, and the entries of
%  opts.ind(1:2,:) are within 1 and length(opts.G).
%
%% Related papers
%
% [1] Jun Liu and Jieping Ye, Moreau-Yosida Regularization for
%     Grouped Tree Structure Learning, NIPS 2010
%
%%

cd ..
cd ..

root=cd;
addpath(genpath([root '/SLEP']));
                     % add the functions in the folder SLEP to the path

% change to the original folder
cd Examples/tree;

m=50;  n=100;       % The data matrix is of size m x n

% ---------------------- generate random data ----------------------
%randn('state',(randNum-1)*3+1);
A=randn(m,n);        % the data matrix

y=[-ones(25,1); ones(25,1)];      % the response

%% In this example, the tree is set as:
%
% root, 1:100, with weight 0
% its children nodes, 1:50, and 51:100
%
% For 1:50, its children are 1:20, 21:40, and 41:50
%
% For 51:100, its children are 51:70, and 71:100
%
% These nodes in addition have each individual features (they contain) as
% children nodes.
%
%%

%% One efficient way
% We make use of the fact that the indices of the left nodes of the tree
% are smaller than the right nodes.

%----------------------- Set optional items -----------------------
opts=[];

% Starting point
opts.init=2;        % starting from a zero point

% Termination
opts.tFlag=5;       % run .maxIter iterations
opts.maxIter=100;   % maximum number of iterations

% regularization
opts.rFlag=1;       % use ratio

% Normalization
opts.nFlag=0;       % without normalization

% Group Property
opts.ind=[[-1, -1, 1]',... % leave nodes (each node contains one feature)
    [1, 20, sqrt(20)]', [21, 40, sqrt(20)]',... % the layer above the leaf
    [41, 50, sqrt(10)]', [51, 70, sqrt(20)]', [71,100, sqrt(30)]',...
    [1, 50, sqrt(50)]', [51, 100, sqrt(50)]']; % the higher layer

%----------------------- Run the code mc_cgLassoLeast -----------------------
z=0.1;
tic;
[x, c, funVal, ValueL]= tree_LogisticR(A, y, z, opts);
toc;

%% An alternative way
% We make use of the fact that the indices of the left nodes of the tree
% are smaller than the right nodes.
%%

%----------------------- Set optional items -----------------------
opts=[];

% Starting point
opts.init=2;        % starting from a zero point

% Termination
opts.tFlag=5;       % run .maxIter iterations
opts.maxIter=100;   % maximum number of iterations

% regularization
opts.rFlag=1;       % use ratio

% Normalization
opts.nFlag=0;       % without normalization

% Group Property
opts.ind=[[-1, -1, 1]',... % leave nodes (each node contains one feature)
    [1, 20, sqrt(20)]', [21, 40, sqrt(20)]',... % the layer above the leaf
    [41, 50, sqrt(10)]', [51, 70, sqrt(20)]', [71,100, sqrt(30)]',...
    [101, 150, sqrt(50)]', [151, 200, sqrt(50)]']; % the higher layer
opts.G=[1:100, 1:100];

%----------------------- Run the code mc_cgLassoLeast -----------------------
z=0.1;
tic;
[x2, c2, funVal2, Value2L]= tree_LogisticR(A, y, z, opts);
toc;

author: zzstefan

repo: BrainNetClass


text:
function [offsetsRows, offsetsCols, distances] = templateMatchingNaive(row, col,...
    patchSize, searchWindowSize, image)
% This function should for each possible offset in the search window
% centred at the current row and col, save a value for the offsets and
% patch distances, e.g. for the offset (-1,-1)
% offsetsRows(1) = -1;
% offsetsCols(1) = -1;
% distances(1) = 0.125;

% The distance is simply the SSD over patches of size patchSize between the
% 'template' patch centred at row and col and a patch shifted by the
% current offset

% NOTE : row and col are the coordinates of the pixel we want our pach to be centered.
% We need the delta to see when we're sliding the patch if we're still inside the image
% or not
delta = floor(patchSize/2);

[rows columns dimensions] = size(image);

% Let's grab the reference patch centered at (row, col)
reference_patch = double(image(row-delta:row+delta, col-delta:col+delta));

% We'll have as many patches as the window Area which is
% (searchWindowSize*searchWindowSize)
% Therefore we'll have the same number for the distances and for the offsets
% in the X and Y direction
distances = zeros(1,searchWindowSize*searchWindowSize);
offsetsRows = zeros(1,searchWindowSize*searchWindowSize);
offsetsCols = zeros(1,searchWindowSize*searchWindowSize);

distances_index = 1;

% We use the delta window to center the window to the pixel we want to
% denoise.
delta_window = floor(searchWindowSize/2);

% CLIP the search window agains the image
% Suppose we want to denoise the pixel x and center the
%
% |------------|
% |     #######|####################
% |     #x     |                   #
% |-----#------|                   #
%       #                          #
%       #                          #
%       #                          #
%       ############################
%
%

start_rows = max(row-delta_window, 1+delta);
start_columns = max(col-delta_window, 1+delta);

end_rows = min(row+delta_window, rows-delta);
end_columns = min(col+delta_window, columns-delta);

% Loop through the clipped window
for row_searchWindow = start_rows:end_rows
    for column_searchWindow = start_columns:end_columns

        %obtain the patch we want to compare with
        patch = double(image(row_searchWindow-delta:row_searchWindow+delta, column_searchWindow-delta:column_searchWindow+delta));

        %compute the difference
        sum_squared_distance = sum(sum(double(reference_patch - patch).^2));

        %store the results
        distances(distances_index) = sum_squared_distance;
        offsetsRows(distances_index) = row_searchWindow - row;
        offsetsCols(distances_index) = column_searchWindow - col;
        distances_index = distances_index + 1;

    end
end

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [new,messages] = searchbk (w, i, thres,freq,thres_alpha,messages)
%performs the searching back when no qrs has been detected for
% 1.5 times the RR before.  i is the beginning sample.  thres
% is the threshold for scale 4;
% Last update: Rute Almeida 27Jan2012

if nargin<6 || (isempty(freq) && (~isfield(messages.setup.wavedet,'freq') || ~(messages.setup.wavedet.freq>0)))
    messages.errors=[messages.errors {'Fatal error in wavedet_3D: no sampling frequency variable found.'}];
    warning(char(messages.errors(end)))
    messages.errors_desc=[messages.errors_desc {'No sampling frequency variable give to seachbk from fiducialf.m.'}];
    new=0;
    return
elseif ~isempty(freq) && isfield(messages.setup.wavedet,'freq') && messages.setup.wavedet.freq>0 && freq~=messages.setup.wavedet.freq
    messages.warnings=[messages.warnings {'sampling frequency used in searchbk is different from messages.setup.wavedet.freq.'}];
elseif isempty(freq)
    freq=messages.setup.wavedet.freq;
end

if ~isfield(messages.setup.wavedet,'refrper') || ~isfield(messages.setup.wavedet,'peakcriteria') || ~isfield(messages.setup.wavedet,'timelapthr') || ~isfield(messages.setup.wavedet,'timelapthr') || ~isfield(messages.setup.wavedet,'intvlthr1') || ~isfield(messages.setup.wavedet,'intvlthr2') || ~isfield(messages.setup.wavedet,'nghbhd')
    messages.errors=[messages.errors {'Fatal error in wavedet_3D: missing parameters on seachbk.'}];
    warning(char(messages.errors(end)))
    messages.errors_desc=[messages.errors_desc {'Setup value is missing on seachbk from fiducialf.m.'}];
    new=0;
    return
end
if ~isfield(messages.setup.wavedet,'intvlthr1_2_sbk')
    messages.setup.wavedet.intvlthr1_2_sbk=2;%reduction on threshold interval for considering redundancy if it subsists
end
if ~isfield(messages.setup.wavedet,'thfraction_sbk')
    if isfield(messages.setup.wavedet,'thfraction')
        messages.setup.wavedet.thfraction_sbk=messages.setup.wavedet.thfraction;
    else
        messages.errors=[messages.errors {'Fatal error in wavedet_3D: missing parameters on seachbk.'}];
        warning(char(messages.errors(end)))
        messages.errors_desc=[messages.errors_desc {'Setup value thfraction is missing on seachbk from fiducialf.m.'}];
        new=0;
        return
    end
end

if ~isfield(messages.setup.wavedet,'nghbhd')
    messages.setup.wavedet.nghbhd=0.025; % neighbourwood fro maximum across scales in sec
end

thfraction_sbk=messages.setup.wavedet.thfraction_sbk;
peakcriteria=messages.setup.wavedet.peakcriteria;
timelap = ceil(messages.setup.wavedet.timelapthr*freq);

thres = thres/2;  % When searching back we take half
thres(4) = thres(4)/2; % so new thres(4) = old thres (4)/4

n = modmax(w(:,4),2,thres(4),0);  % Search maximum moduli at scale 4
neighb = ceil(freq*messages.setup.wavedet.nghbhd);        % Neighbourhood = 25 ms

m = zeros(4,length(n));
signo = sign(w(n,4))';
if ~isempty(n),
    m(4,:) = n';
end
new = []; %#ok<NASGU>

for k = 1:size(m,2),
    window = [max(n(k)-neighb,1), min(n(k)+neighb,size(w,1))];
    n3 = modmax(w(window(1):window(2),3),2,thres(3),signo(k));
    n3 =window(1)-1+n3;
    num = length(n3);
    if num>0,
        if num==1,
            m(3,k)= n3;
        elseif num>1
            if length(find(max(abs(w(n3,3)))./abs(w(n3,3))<peakcriteria))==1,
                [aux1,ind]=max(abs(w(n3,3))); %#ok<ASGLU> % greatest modulus
                m(3,k)= n3(ind);
            else                              % minimum distance
                [aux1,ind]=min(abs(m(4,k)-n3)); %#ok<ASGLU>
                m(3,k)= n3(ind);
            end
        end
    end
end
ind = find(m(3,:)==0);
m(:,ind) = [];
signo(ind) = [];

% Search for maximum moduli in the neighborhood at scale 2
for k = 1:size(m,2),
    window = [max(m(3,k)-neighb,1), min(m(3,k)+neighb,size(w,1))];
    n2 =modmax(w(window(1):window(2),2),2,thres(2),signo(k));
    n2 =window(1)-1+n2;
    num = length(n2);
    if num>0,
        if num==1,
            m(2,k)= n2;
        elseif num>1
            if length(find(max(abs(w(n2,2)))./abs(w(n2,2))<peakcriteria))==1,
                [aux1,ind]=max(abs(w(n2,2))); %#ok<ASGLU> % greatest modulus
                m(2,k)= n2(ind);
            else                            % shortest distance
                [aux1,ind]=min(abs(m(3,k)-n2)); %#ok<ASGLU>
                m(2,k)= n2(ind);
            end
        end
    end
end
ind = find(m(2,:)==0);
m(:,ind) = [];
signo(ind) = [];

for k = 1:size(m,2),
    window = [max(m(2,k)-neighb,1), min(m(2,k)+neighb,size(w,1))];
    n1 =modmax(w(window(1):window(2),1),2,thres(1),signo(k));
    n1 =window(1)-1+n1;
    num = length(n1);
    if num>0,
        if num==1,
            m(1,k)= n1;
        elseif num>1
            if length(find(max(abs(w(n1,1)))./abs(w(n1,1))<peakcriteria))==1,
                [aux1,ind]=max(abs(w(n1,1)));%#ok<ASGLU> % greatest modulus
                m(1,k)= n1(ind);
            else                             % shortest distance
                [aux1,ind]=min(abs(m(2,k)-n1)); %#ok<ASGLU>
                m(1,k)= n1(ind);
            end
        end
    end
end

ind = find(m(1,:)==0);               %Discard all maximum lines with no
m(:,ind) = [];                       % associated maximum at scale 1
signo(ind) = [];

% Regularity Exponent Validation
% alpha proportional to log(a3(nk3))-log(a1(nk1))
alpha = log(abs(w(m(3,:),3)));  %%!!!!
% alpha = log(abs(w(m(3,:),3))) - log (abs(w(m(1,:),1)));  !!!
ind = find(alpha <= thres_alpha-thfraction_sbk);   %%%% !!!
m(:,ind) = [];
signo(ind) = [];

thresinterval = ceil(messages.setup.wavedet.intvlthr2 * freq);         % 120 ms. (Li)

if size(m,2)>2,
    ind = find( ((m(1,2:end-1)-m(1,1:end-2))>thresinterval) ...
        &      (m(1,3:end)- m(1,2:end-1))>thresinterval)+1;
    if (m(1,2)-m(1,1))>thresinterval,
        ind = [1 ind];
    end
    if (m(1,end)-m(1,end-1))>thresinterval,
        ind = [ind size(m,2)];
    end
    m(:,ind) =[];                     % Discard isolated maximum lines
    signo(ind) = [];

elseif size(m,2)==2,                 % If only two lines
    if (m(1,2)-m(1,1))>thresinterval, % discard them if too separated
        m(:,1:2)=[];
        signo(1:2)=[];
    end
elseif size(m,2)==1,                % If only one, discard it
    m(:,1)=[];
    signo(1)=[];
end

% Threshold interval for considering redundancy
redundant= [];
thresinterval = ceil(messages.setup.wavedet.intvlthr1*freq);           % 120 ms. (li)

for l = find(signo>0),                      % For each positive maximum line
    if ~any(redundant ==l),                   % If it has not been declared redundant yet
        ind=find((m(3,:)>m(3,l)-messages.setup.wavedet.intvlthr1_2_sbk*thresinterval)&(m(3,:)<m(3,l)+messages.setup.wavedet.intvlthr1_2_sbk*thresinterval)&signo>0);
        % index of positive lines near the present one (including it)
        if length(ind)>1,                        % If more than one --> redundancy
            [mx,ind2]= max(abs(w(m(3,ind),3))); %#ok<ASGLU>
            ind(ind2)=[];
            redundant = [redundant ind];             %#ok<AGROW> % All but the greatest are redundant
        end
    end
end

m(:,redundant)=[];                         % Discard redundant lines
signo(redundant)=[];
redundant= [];

for l = find(signo>0),              % For each remaining positive maximum line
    ind = find((m(3,:)>m(3,l)-thresinterval)&(m(3,:)<m(3,l)+thresinterval)&signo<0);
    % Search for negative minima near it
    if length(ind)>1,               % If more than one ---> redundancy
        aux = abs(w(m(3,ind),3)./(m(3,ind)'-m(3,l)));
        % auxiliary variable: height over distance
        [mx,indmx]=max(aux);
        ind2=ind;
        aux(indmx)=[]; ind2(indmx)=[];
        if all((mx./aux)>peakcriteria),        % RULE 2 (see PFC or Li's paper)
            redundant = [redundant ind2]; %#ok<AGROW>
        else
            [aux,aux2]=min(abs(m(3,l)-m(3,ind))); %#ok<ASGLU>
            ind(aux2)=[];              % RULE 1 (see PFC or Li's paper)
            redundant = [redundant ind]; %#ok<AGROW>
        end
    end
end
m(:,redundant)=[];                  % Discard redundant lines
signo(redundant)=[];
redundant = [];
for l = find(signo<0),              % For each remaining negative minimum line
    ind = find((m(3,:)>m(3,l)-thresinterval)&(m(3,:)<m(3,l)+thresinterval)&signo>0);
    % Search for positive maxima near it
    if length(ind)>1,                % If more than one ----> redundancy
        aux = abs(w(m(3,ind),3)./(m(3,ind)'-m(3,l)));
        % auxiliary variable: height over distance
        [mx,indmx]=max(aux);
        ind2=ind;
        aux(indmx)=[]; ind2(indmx)=[];
        if all((mx./aux)>peakcriteria),
            redundant = [redundant ind2]; %#ok<AGROW>   % RULE 2
        else
            [aux,aux2]=min(abs(m(3,l)-m(3,ind))); %#ok<ASGLU>
            ind(aux2)=[];
            redundant = [redundant ind]; %#ok<AGROW>    % RULE 1
        end
    end
end
m(:,redundant)=[];
signo(redundant)=[];

%%%%%%%%%%%%%% isolated maximum lines resulting from Discarded
%%%%%%%%%%%%%% redundant OUT1011

ind = find( ((m(1,2:end-1)-m(1,1:end-2))>thresinterval) ...
    &      (m(1,3:end)- m(1,2:end-1))>thresinterval)+1;

if size(m,2)<2 && ~isempty(ind), ind=1; end  % when only one maximum %16DEZ08
m(:,ind) = [];                     %Discard isolated maximum lines
signo(ind)=[];
if size(m,2)<2
    m=[];
    signo=[];
end

eliminar=[];
%%%%extra protection% OUT2011
for ii=1:size(m,2)
    pa = picant(w(max(1,m(2,ii)-round(messages.setup.wavedet.pictime*freq)):m(2,ii),2),m(2,ii));
    % first peak before detected qrs position at scale 2
    pp = picpost(w(m(2,ii):min(size(w,1),m(2,ii)+round(messages.setup.wavedet.pictime*freq)),2),m(2,ii));

    if isempty(pa) || isempty(pp)
        eliminar=[eliminar ii]; %#ok<AGROW>
    end
end

m(:,eliminar)=[];
signo(eliminar)=[];

% QRS peak detection / wavelet Zero cross detection
time = [];
aux=[]; %Rute26Jun09
if length(signo)>1,
    for l = find(signo>0),
        if (l==1)
            if (signo(2)<0)&&(m(1,2)-m(1,1)<thresinterval),
                ind = zerocros(w(m(1,l):min(m(1,l)+timelap,size(w,1)),1));
                % Zero crossing at scale 1
                time = [time ind+m(1,l)-1]; %#ok<AGROW>
                aux=[aux abs(w(m(1,l))-w(m(1,l+1)))];  %#ok<AGROW> %Rute26Jun09
            end
        elseif (l==size(m,2))                % Special case: last line
            if (signo(l-1)<0)&&(m(1,end)-m(1,end-1)<thresinterval),
                ind = zerocros(w(m(1,l-1):min(m(1,l-1)+timelap,size(w,1)),1));
                time = [time ind+m(1,l-1)-1]; %#ok<AGROW>
                aux=[aux abs(w(m(1,l-1))-w(m(1,l)))];  %#ok<AGROW>  %Rute26Jun09
            end
        elseif signo(l+1)<0 && ((signo(l-1)>0) ...
                || ((m(1,l+1)-m(1,l))<(m(1,l)-m(1,l-1))))
            ind = zerocros(w(m(1,l):min(m(1,l)+timelap,size(w,1)),1));
            time = [time ind+m(1,l)-1]; %#ok<AGROW>
            aux=[aux abs(w(m(1,l))-w(m(1,l+1)))]; %#ok<AGROW>  %Rute26Jun09
        elseif signo(l-1)<0,
            ind = zerocros(w(m(1,l-1):min(m(1,l-1)+timelap,size(w,1)),1));
            time = [time ind+m(1,l-1)-1]; %#ok<AGROW>
            aux=[aux abs(w(m(1,l-1))-w(m(1,l)))]; %#ok<AGROW> %Rute 26Jun09
        end
    end
end

% Refractary period after a QRS detection (200 ms).
rr = (time(2:end)-time(1:end-1));
ind = find(rr<ceil(messages.setup.wavedet.refrper(end)*freq));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%RUTE 26Jun09
%time(ind+1)=[]; % always the second one is eliminated!!! to be changed
for auxi=1:length(ind) %RUTE 27Jun11
    [M,ii]=min([aux(ind(auxi)) aux(ind(auxi)+1)]); %#ok<ASGLU> % 18JUL2011
    if ii==2, ind(auxi)=ind(auxi)+1; end
end
time(ind)=[]; %RUTE 27Jun11
aux(ind)=[]; %#ok<NASGU>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
new = time +i -1;

author: marianux

repo: ecg-kit

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


text:
function varargout = chebfuneps(val)
%CHEBFUNEPS   Set the default value of the chebfuneps preference.
%   CHEBFUNEPS VAL, or CHEBFUNEPS(VAL) sets the default value of the
%   chebfuneps preference to the specified value VAL. CHEBFUNEPS(VAL) is
%   equivalent to CHEBFUNPREF.SETDEFAULTS('chebfuneps', VAL).
%
%   CHEBFUNEPS factory, or CHEBFUNEPS('factory') sets the default
%   chebfuneps preference to the factory value. This is equivalent to
%   chebfunpref.setDefaults('chebfuneps', 'factory').
%
%   CHEBFUNEPS prints the current default value of the chebfuneps
%   preference.
%
%   If changing the preference is needed only for a single construction,
%   calling constructor with the 'eps' flag is a better option.
%
% See also CHEBFUN2EPS and CHEBFUN3EPS.

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

if ( nargin == 0 )
    % Return current chebfuneps:
    varargout{1} = chebfunpref().techPrefs.chebfuneps;
else
    if ( strcmpi(val, 'factory') )
        chebfunpref.setDefaults('chebfuneps', 'factory');
    else
        eval(['val=',val,';'])
        if ( isnumeric(val) )
            chebfunpref.setDefaults('chebfuneps', val);
        else
            error('CHEBFUN:chebfuneps:unknownOption',...
                'Unknown CHEBFUNEPS option.')
        end
    end
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% function [Lbdi, Rbdi] = synsq_bd_loc(t, fs)
%
% For a given time series and set of frequency values, calculate the left
% and right boundaries for each frequency.  Beyond these boundaries,
% boundary effects affect the signal.  Output is in indices.
%   If length(t)==n, then length(Lbdi)==length(Rbdi)==n
%
% t is assumed to be sampled uniformly
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------
function [Lbdi, Rbdi] = synsq_bd_loc(t, fs)
    if (fs(2)>fs(1))
        Ts = 1./fs;
    else
        Ts = fs;
    end

    dt = t(2)-t(1);
    Tsn = ceil(1.5 * Ts / dt);

    Lbdi = min(length(t), Tsn);
    Rbdi = max(1, length(t)-Tsn);
end

author: HeLiangHIT


text:
function [ bvec3, overflow ] = bvec_add ( n, bvec1, bvec2 )

%*****************************************************************************80
%
%% BVEC_ADD adds two binary vectors.
%
%  Discussion:
%
%    A BVEC is an integer vector of binary digits, intended to
%    represent an integer.  BVEC(1) is the units digit, BVEC(N-1)
%    is the coefficient of 2**(N-2), and BVEC(N) contains sign
%    information.  It is 0 if the number is positive, and 1 if
%    the number is negative.
%
%  Example:
%
%    N = 4
%
%    BVEC1    dec  BVEC2    dec  BVEC3    dec
%    -------  ---  -------  ---  -------  ---
%    1 0 0 0   1   1 1 0 0   3   0 0 1 0   4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the length of the vectors.
%
%    Input, integer BVEC1(N), BVEC2(N), the vectors to be added.
%
%    Output, integer BVEC3(N), the sum of the two input vectors.
%
%    Output, logical OVERFLOW, is true if the sum overflows.
%
  base = 2;
  overflow = 0;

  bvec3(1:n) = floor ( bvec1(1:n) ) + floor ( bvec2(1:n) );

  for i = 1 : n
    while ( base <= bvec3(i) )
      bvec3(i) = bvec3(i) - base;
      if ( i < n )
        bvec3(i+1) = bvec3(i+1) + 1;
      else
        overflow = 1;
      end
    end
  end

  return
end


text:
function [J, grad] = cofiCostFunc(params, Y, R, num_users, num_movies, ...
                                  num_features, lambda)
%COFICOSTFUNC Collaborative filtering cost function
%   [J, grad] = COFICOSTFUNC(params, Y, R, num_users, num_movies, ...
%   num_features, lambda) returns the cost and gradient for the
%   collaborative filtering problem.
%

% Unfold the U and W matrices from params
X = reshape(params(1:num_movies*num_features), num_movies, num_features);
Theta = reshape(params(num_movies*num_features+1:end), ...
                num_users, num_features);

% You need to return the following values correctly
J = 0;
X_grad = zeros(size(X));
Theta_grad = zeros(size(Theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost function and gradient for collaborative
%               filtering. Concretely, you should first implement the cost
%               function (without regularization) and make sure it is
%               matches our costs. After that, you should implement the
%               gradient and use the checkCostFunction routine to check
%               that the gradient is correct. Finally, you should implement
%               regularization.
%
% Notes: X - num_movies  x num_features matrix of movie features
%        Theta - num_users  x num_features matrix of user features
%        Y - num_movies x num_users matrix of user ratings of movies
%        R - num_movies x num_users matrix, where R(i, j) = 1 if the
%            i-th movie was rated by the j-th user
%
% You should set the following variables correctly:
%
%        X_grad - num_movies x num_features matrix, containing the
%                 partial derivatives w.r.t. to each element of X
%        Theta_grad - num_users x num_features matrix, containing the
%                     partial derivatives w.r.t. to each element of Theta
%

H_Theta = X * Theta';
Error = H_Theta - Y;
Error_logical = Error .* R;
J = 1/2 * sum(Error_logical(:).^2) + ...
    lambda/2*Theta(:)'*Theta(:) + ...
    lambda/2*X(:)'*X(:);

X_grad = Error_logical*Theta + lambda*X;
Theta_grad = Error_logical'*X + lambda*Theta;

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

grad = [X_grad(:); Theta_grad(:)];

end

author: zlotus


text:
function [ODS,OIS,AP] = stEvalBsds( model, varargin )
% Evaluate sketch token edge detector on BSDS500.
%
% USAGE
%  [ODS,OIS,AP] = stEvalBsds( model, parameters )
%
% INPUTS
%  model      - sketch token model or model filename
%  parameters - parameters (struct or name/value pairs)
%   .dataType   - ['test'] should be either 'test' or 'val'
%   .nThresh    - [99] number of thresholds for evaluation
%   .cleanup    - [0] if true delete temporary files
%   .show       - [0] figure for displaying results (or 0)
%   .modelDir   - [] directory for storing models
%   .bsdsDir    - [] directory of BSDS dataset
%   .name       - [''] name to append to evaluation
%   .pDistr     - [{'type','parfor'}] parameters for fevalDistr
%
% OUTPUTS
%  ODS        - standard error measure on BSDS500
%  OIS        - standard error measure on BSDS500
%  AP         - standard error measure on BSDS500
%
% EXAMPLE
%
% See also stDetect, stTrain
%
% Sketch Token Toolbox     V0.95
% Copyright 2013 Joseph Lim [lim@csail.mit.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see bsd.txt]

    % get default parameters
    dfs={'dataType','test', 'nThresh',99, 'cleanup',0, 'show',0, ...
      'modelDir',[], 'bsdsDir',[], 'name','', 'pDistr',{{'type','parfor'}} };
    p=getPrmDflt(varargin,dfs,1);
    if( ischar(model) ),
        model=load(model);
        model=model.model;
    end
    if( isempty(p.modelDir )),
        p.modelDir=model.opts.modelDir;
    end
    if( isempty(p.bsdsDir )),
        p.bsdsDir=model.opts.bsdsDir;
    end
    p.modelDir = [p.modelDir '/' p.dataType '/'];

    % eval on either validation set or test set
    imgDir = [p.bsdsDir '/images/' p.dataType '/'];
    gtDir = [p.bsdsDir '/groundTruth/' p.dataType '/'];
    evalDir = [p.modelDir model.opts.modelFnm p.name '-eval/'];
    resDir = [p.modelDir model.opts.modelFnm p.name '/'];
    assert(exist(imgDir,'dir')==7);
    assert(exist(gtDir,'dir')==7);

    % if evaluation exists collect results and display
    if(exist([evalDir '/eval_bdry.txt'],'file'))
        [ODS,~,~,~,OIS,~,~,AP]=collect_eval_bdry(evalDir);
        fprintf('ODS=%.3f OIS=%.3f AP=%.3f\n',ODS,OIS,AP);
        if( p.show ),
            plot_eval(evalDir,'r');
        end;
        return;
    end

    % get image ids
    ids=dir([imgDir '*.jpg']);
    ids={ids.name};
    n=length(ids);
    for i=1:n,
        ids{i}=ids{i}(1:end-4);
    end

    % extract sketch tokens and save results
    if(~exist(resDir,'dir')),
        mkdir(resDir);
    end;
    do=false(1,n);

    for i=1:n,
        do(i)=~exist([resDir ids{i} '.png'],'file');
    end
    do=find(do);
    m=length(do);
    parfor i=1:m,
        id=ids{do(i)}; %#ok<PFBNS>

        I = imread([imgDir id '.jpg']);
        st = max(min(1,stDetect(I,model)),0);
        E = stToEdges(st,1);

        imwrite(uint8(E*255),[resDir id '.png']);
    end

    % perform evaluation on each image (Linux only, slow)
    if(ispc),
        error('Evaluation code runs on Linux ONLY.');
    end
    do=false(1,n);
    jobs=cell(1,n);
    for i=1:n,
        do(i)=~exist([evalDir ids{i} '_ev1.txt'],'file');
    end
    for i=1:n,
        id=ids{i};
        jobs{i}={[resDir id '.png'],...
        [gtDir id '.mat'],[evalDir id '_ev1.txt'],p.nThresh};
    end
    if(~exist(evalDir,'dir')),
        mkdir(evalDir);
    end
    fevalDistr('evaluation_bdry_image',jobs(do),p.pDistr{:});

    % collect results and display
    [ODS,~,~,~,OIS,~,~,AP]=collect_eval_bdry(evalDir);
    fprintf('ODS=%.3f OIS=%.3f AP=%.3f\n',ODS,OIS,AP);
    if( p.show ),
        plot_eval(evalDir,'r');
    end
    if( p.cleanup ),
        delete([evalDir '/*_ev1.txt']);
        delete([evalDir '/eval_bdry_img.txt']);
        delete([evalDir '/eval_bdry_thr.txt']);
        delete([resDir '/*.png']); rmdir(resDir);
    end

end

author: zouchuhang

repo: LayoutNet


text:
function visible_probability = hidden_state_to_visible_probabilities(rbm_w, hidden_state)
  % <rbm_w> is a matrix of size <number of hidden units> by <number of visible units>
  % <hidden_state> is a binary matrix of size <number of hidden units> by <number of configurations that we're handling in parallel>.
  % The returned value is a matrix of size <number of visible units> by <number of configurations that we're handling in parallel>.
  % This takes in the (binary) states of the hidden units, and returns the activation probabilities of the visible units, conditional on those states.
  m = size(hidden_state, 2);
  number_hidden = size(rbm_w, 1);
  number_visible = size(rbm_w, 2);
  visible_probability = zeros(number_visible, m);

  for k = 1:m
    hidden_state_k = hidden_state(:, k);
    visible_probability(:, k) = 1 ./ (1 + exp(-(rbm_w' * hidden_state_k)));
  end
end

author: BradNeuberg

repo: hinton-coursera


text:
function [A_hat E_hat iter] = exact_alm_rpca(D, lambda, tol, maxIter)

% Oct 2009
% This matlab code implements the augmented Lagrange multiplier method for
% Robust PCA.
%
% D - m x n matrix of observations/data (required input)
%
% lambda - weight on sparse error term in the cost function
%
% tol - tolerance for stopping criterion.
%     - DEFAULT 1e-7 if omitted or -1.
%
% maxIter - maximum number of iterations
%         - DEFAULT 1000, if omitted or -1.
%
% Initialize A,E,Y,u
% while ~converged
%   minimize
%     L(A,E,Y,u) = |A|_* + lambda * |E|_1 + <Y,D-A-E> + mu/2 * |D-A-E|_F^2;
%   Y = Y + \mu * (D - A - E);
%   \mu = \rho * \mu;
% end
%
% Minming Chen, October 2009. Questions? v-minmch@microsoft.com ;
% Arvind Ganesh (abalasu2@illinois.edu)
%
% Copyright: Perception and Decision Laboratory, University of Illinois, Urbana-Champaign
%            Microsoft Research Asia, Beijing

% addpath PROPACK;

[m n] = size(D);

if nargin < 2
    lambda = 1 / sqrt(m);
end

if nargin < 3
    tol = 1e-7;
elseif tol == -1
    tol = 1e-7;
end

if nargin < 4
    maxIter = 1000;
elseif maxIter == -1
    maxIter = 1000;
end

% initialize
Y = sign(D);
%norm_two = lansvd(Y, 1, 'L');
norm_two = svdex(Y, 1, struct('svd_solver','svds'));
norm_inf = norm( Y(:), inf) / lambda;
dual_norm = max(norm_two, norm_inf);
Y = Y / dual_norm;

A_hat = zeros( m, n);
E_hat = zeros( m, n);
dnorm = norm(D, 'fro');
tolProj = 1e-6 * dnorm;
total_svd = 0;
mu = .5/norm_two; % this one can be tuned
rho = 6;          % this one can be tuned

iter = 0;
converged = false;
stopCriterion = 1;
sv = 5;
svp = sv;

% excelly
temp_T = D - A_hat + (1/mu)*Y;
temp_E = max( temp_T - lambda/mu,0) + min( temp_T + lambda/mu,0);
[U dummy dummy] = svdex(D - temp_E + (1/mu)*Y, 1, struct('svd_solver','svds'));
while ~converged
    iter = iter + 1;

    % solve the primal problem by alternative projection
    primal_converged = false;
    primal_iter = 0;
    sv = sv + round(n * 0.1);
    while primal_converged == false

        temp_T = D - A_hat + (1/mu)*Y;
        temp_E = max( temp_T - lambda/mu,0) + min( temp_T + lambda/mu,0);

        [U diagS V] = svdex(D - temp_E + (1/mu)*Y, sv, struct('init',U(:,1)));
        % if choosvd(n, sv) == 1
        %     [U S V] = lansvd(D - temp_E + (1/mu)*Y, sv, 'L');
        % else
        %     [U S V] = svd(D - temp_E + (1/mu)*Y, 'econ');
        % end
        % diagS = diag(S);
        svp = length(find(diagS > 1/mu));
        if svp < sv
            sv = min(svp + 1, n);
        else
            sv = min(svp + round(0.05*n), n);
        end
        temp_A = U(:,1:svp)*diag(diagS(1:svp)-1/mu)*V(:,1:svp)';

        if norm(A_hat - temp_A, 'fro') < tolProj && norm(E_hat - temp_E, 'fro') < tolProj
            primal_converged = true;
        end
        A_hat = temp_A;
        E_hat = temp_E;
        primal_iter = primal_iter + 1;
        total_svd = total_svd + 1;

    end

    Z = D - A_hat - E_hat;
    Y = Y + mu*Z;
    mu = rho * mu;

    %% stop Criterion
    stopCriterion = norm(Z, 'fro') / dnorm;
    if stopCriterion < tol
        converged = true;
    end

    disp(['Iteration' num2str(iter) ' #svd ' num2str(total_svd) ' r(A) ' num2str(svp)...
        ' |E|_0 ' num2str(length(find(abs(E_hat)>0)))...
        ' stopCriterion ' num2str(stopCriterion)]);

    if ~converged && iter >= maxIter
        disp('Maximum iterations reached') ;
        converged = 1 ;
    end
end

if nargin == 5
    fclose(fid);
end

author: andrewssobral

repo: lrslibrary


text:
function writetoPAJ(CIJ, fname, arcs)
%WRITETOPAJ         Write to Pajek
%
%   writetoPAJ(CIJ, fname, arcs);
%
%   This function writes a Pajek .net file from a MATLAB matrix
%
%   Inputs:     CIJ,        adjacency matrix
%               fname,      filename minus .net extension
%               arcs,       1 for directed network
%                           0 for an undirected network
%
%   Chris Honey, Indiana University, 2007

N = size(CIJ,1);
fid = fopen(cat(2,fname,'.net'), 'w');

%%%VERTICES
fprintf(fid, '*vertices %6i \r', N);
for i = 1:N
    fprintf(fid, '%6i "%6i" \r', [i i]);
end

%%%ARCS/EDGES
if arcs
    fprintf(fid, '*arcs \r');
else
    fprintf(fid, '*edges \r');
end

for i = 1:N
    for j = 1:N
        if CIJ(i,j) ~= 0
            fprintf(fid, '%6i %6i %6f \r', [i j CIJ(i,j)]);
        end
    end
end

fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [net,stats] = cnn_train_dag(net, imdb, getBatch, varargin)
%CNN_TRAIN_DAG Demonstrates training a CNN using the DagNN wrapper
%    CNN_TRAIN_DAG() is similar to CNN_TRAIN(), but works with
%    the DagNN wrapper instead of the SimpleNN wrapper.

% Copyright (C) 2014-16 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).
addpath(fullfile(vl_rootnn, 'examples'));

opts.expDir = fullfile('data','exp') ;
opts.continue = true ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.train = [] ;
opts.val = [] ;
opts.gpus = [] ;
opts.prefetch = false ;
opts.epochSize = inf;
opts.numEpochs = 300 ;
opts.learningRate = 0.001 ;
opts.weightDecay = 0.0005 ;

opts.solver = [] ;  % Empty array means use the default SGD solver
[opts, varargin] = vl_argparse(opts, varargin) ;
if ~isempty(opts.solver)
  assert(isa(opts.solver, 'function_handle') && nargout(opts.solver) == 2,...
    'Invalid solver; expected a function handle with two outputs.') ;
  % Call without input arguments, to get default options
  opts.solverOpts = opts.solver() ;
end

opts.momentum = 0.9 ;
opts.saveSolverState = true ;
opts.nesterovUpdate = false ;
opts.randomSeed = 0 ;
opts.profile = false ;
opts.parameterServer.method = 'mmap' ;
opts.parameterServer.prefix = 'mcn' ;

opts.derOutputs = {'objective', 1} ;
opts.extractStatsFn = @extractStats ;
opts.plotStatistics = true;
opts.postEpochFn = [] ;  % postEpochFn(net,params,state) called after each epoch; can return a new learning rate, 0 to stop, [] for no change
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isscalar(opts.train) && isnumeric(opts.train) && isnan(opts.train)
  opts.train = [] ;
end
if isscalar(opts.val) && isnumeric(opts.val) && isnan(opts.val)
  opts.val = [] ;
end

% -------------------------------------------------------------------------
%                                                            Initialization
% -------------------------------------------------------------------------

evaluateMode = isempty(opts.train) ;
if ~evaluateMode
  if isempty(opts.derOutputs)
    error('DEROUTPUTS must be specified when training.\n') ;
  end
end

% -------------------------------------------------------------------------
%                                                        Train and validate
% -------------------------------------------------------------------------

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('%s: resuming by loading epoch %d\n', mfilename, start) ;
  [net, state, stats] = loadState(modelPath(start)) ;
else
  state = [] ;
end

for epoch=start+1:opts.numEpochs

  % Set the random seed based on the epoch and opts.randomSeed.
  % This is important for reproducibility, including when training
  % is restarted from a checkpoint.

  rng(epoch + opts.randomSeed) ;
  prepareGPUs(opts, epoch == start+1) ;

  % Train for one epoch.
  params = opts ;
  params.epoch = epoch ;
  params.learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  params.train = opts.train(randperm(numel(opts.train))) ; % shuffle
  params.train = params.train(1:min(opts.epochSize, numel(opts.train)));
  params.val = opts.val(randperm(numel(opts.val))) ;
  params.imdb = imdb ;
  params.getBatch = getBatch ;

  if numel(opts.gpus) <= 1
    [net, state] = processEpoch(net, state, params, 'train') ;
    [net, state] = processEpoch(net, state, params, 'val') ;
    if ~evaluateMode
      saveState(modelPath(epoch), net, state) ;
    end
    lastStats = state.stats ;
  else
    spmd
      [net, state] = processEpoch(net, state, params, 'train') ;
      [net, state] = processEpoch(net, state, params, 'val') ;
      if labindex == 1 && ~evaluateMode
        saveState(modelPath(epoch), net, state) ;
      end
      lastStats = state.stats ;
    end
    lastStats = accumulateStats(lastStats) ;
  end

  stats.train(epoch) = lastStats.train ;
  stats.val(epoch) = lastStats.val ;
  clear lastStats ;
  saveStats(modelPath(epoch), stats) ;

  if opts.plotStatistics
    switchFigure(1) ; clf ;
    plots = setdiff(...
      cat(2,...
      fieldnames(stats.train)', ...
      fieldnames(stats.val)'), {'num', 'time'}) ;
    for p = plots
      p = char(p) ;
      values = zeros(0, epoch) ;
      leg = {} ;
      for f = {'train', 'val'}
        f = char(f) ;
        if isfield(stats.(f), p)
          tmp = [stats.(f).(p)] ;
          values(end+1,:) = tmp(1,:)' ;
          leg{end+1} = f ;
        end
      end
      subplot(1,numel(plots),find(strcmp(p,plots))) ;
      plot(1:epoch, values','o-') ;
      xlabel('epoch') ;
      title(p) ;
      legend(leg{:}) ;
      grid on ;
    end
    drawnow ;
    print(1, modelFigPath, '-dpdf') ;
  end

  if ~isempty(opts.postEpochFn)
    if nargout(opts.postEpochFn) == 0
      opts.postEpochFn(net, params, state) ;
    else
      lr = opts.postEpochFn(net, params, state) ;
      if ~isempty(lr), opts.learningRate = lr; end
      if opts.learningRate == 0, break; end
    end
  end
end

% With multiple GPUs, return one copy
if isa(net, 'Composite'), net = net{1} ; end

% -------------------------------------------------------------------------
function [net, state] = processEpoch(net, state, params, mode)
% -------------------------------------------------------------------------
% Note that net is not strictly needed as an output argument as net
% is a handle class. However, this fixes some aliasing issue in the
% spmd caller.

% initialize with momentum 0
if isempty(state) || isempty(state.solverState)
  state.solverState = cell(1, numel(net.params)) ;
  state.solverState(:) = {0} ;
end

% move CNN  to GPU as needed
numGpus = numel(params.gpus) ;
if numGpus >= 1
  net.move('gpu') ;
  for i = 1:numel(state.solverState)
    s = state.solverState{i} ;
    if isnumeric(s)
      state.solverState{i} = gpuArray(s) ;
    elseif isstruct(s)
      state.solverState{i} = structfun(@gpuArray, s, 'UniformOutput', false) ;
    end
  end
end
if numGpus > 1
  parserv = ParameterServer(params.parameterServer) ;
  net.setParameterServer(parserv) ;
else
  parserv = [] ;
end

% profile
if params.profile
  if numGpus <= 1
    profile clear ;
    profile on ;
  else
    mpiprofile reset ;
    mpiprofile on ;
  end
end

num = 0 ;
epoch = params.epoch ;
subset = params.(mode) ;
adjustTime = 0 ;

stats.num = 0 ; % return something even if subset = []
stats.time = 0 ;

start = tic ;
for t=1:params.batchSize:numel(subset)
  fprintf('%s: epoch %02d: %3d/%3d:', mode, epoch, ...
          fix((t-1)/params.batchSize)+1, ceil(numel(subset)/params.batchSize)) ;
  batchSize = min(params.batchSize, numel(subset) - t + 1) ;

  for s=1:params.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+params.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
    num = num + numel(batch) ;
    if numel(batch) == 0, continue ; end

    inputs = params.getBatch(params.imdb, batch) ;

    if params.prefetch
      if s == params.numSubBatches
        batchStart = t + (labindex-1) + params.batchSize ;
        batchEnd = min(t+2*params.batchSize-1, numel(subset)) ;
      else
        batchStart = batchStart + numlabs ;
      end
      nextBatch = subset(batchStart : params.numSubBatches * numlabs : batchEnd) ;
      params.getBatch(params.imdb, nextBatch) ;
    end

    if strcmp(mode, 'train')
      net.mode = 'normal' ;
      net.accumulateParamDers = (s ~= 1) ;
      net.eval(inputs, params.derOutputs, 'holdOn', s < params.numSubBatches) ;
    else
      net.mode = 'test' ;
      net.eval(inputs) ;
    end
  end

  % Accumulate gradient.
  if strcmp(mode, 'train')
    if ~isempty(parserv), parserv.sync() ; end
    state = accumulateGradients(net, state, params, batchSize, parserv) ;
  end

  % Get statistics.
  time = toc(start) + adjustTime ;
  batchTime = time - stats.time ;
  stats.num = num ;
  stats.time = time ;
  stats = params.extractStatsFn(stats,net) ;
  currentSpeed = batchSize / batchTime ;
  averageSpeed = (t + batchSize - 1) / time ;
  if t == 3*params.batchSize + 1
    % compensate for the first three iterations, which are outliers
    adjustTime = 4*batchTime - time ;
    stats.time = time + adjustTime ;
  end

  fprintf(' %.1f (%.1f) Hz', averageSpeed, currentSpeed) ;
  for f = setdiff(fieldnames(stats)', {'num', 'time'})
    f = char(f) ;
    fprintf(' %s: %.3f', f, stats.(f)) ;
  end
  fprintf('\n') ;
end

% Save back to state.
state.stats.(mode) = stats ;
if params.profile
  if numGpus <= 1
    state.prof.(mode) = profile('info') ;
    profile off ;
  else
    state.prof.(mode) = mpiprofile('info');
    mpiprofile off ;
  end
end
if ~params.saveSolverState
  state.solverState = [] ;
else
  for i = 1:numel(state.solverState)
    s = state.solverState{i} ;
    if isnumeric(s)
      state.solverState{i} = gather(s) ;
    elseif isstruct(s)
      state.solverState{i} = structfun(@gather, s, 'UniformOutput', false) ;
    end
  end
end

net.reset() ;
net.move('cpu') ;

% -------------------------------------------------------------------------
function state = accumulateGradients(net, state, params, batchSize, parserv)
% -------------------------------------------------------------------------
numGpus = numel(params.gpus) ;
otherGpus = setdiff(1:numGpus, labindex) ;

for p=1:numel(net.params)

  if ~isempty(parserv)
    parDer = parserv.pullWithIndex(p) ;
  else
    parDer = net.params(p).der ;
  end

  switch net.params(p).trainMethod

    case 'average' % mainly for batch normalization
      thisLR = net.params(p).learningRate ;
      net.params(p).value = vl_taccum(...
          1 - thisLR, net.params(p).value, ...
          (thisLR/batchSize/net.params(p).fanout),  parDer) ;

    case 'gradient'
      thisDecay = params.weightDecay * net.params(p).weightDecay ;
      thisLR = params.learningRate * net.params(p).learningRate ;

      if thisLR>0 || thisDecay>0
        % Normalize gradient and incorporate weight decay.
        parDer = vl_taccum(1/batchSize, parDer, ...
                           thisDecay, net.params(p).value) ;

        if isempty(params.solver)
          % Default solver is the optimised SGD.
          % Update momentum.
          state.solverState{p} = vl_taccum(...
            params.momentum, state.solverState{p}, ...
            -1, parDer) ;

          % Nesterov update (aka one step ahead).
          if params.nesterovUpdate
            delta = params.momentum * state.solverState{p} - parDer ;
          else
            delta = state.solverState{p} ;
          end

          % Update parameters.
          net.params(p).value = vl_taccum(...
            1,  net.params(p).value, thisLR, delta) ;

        else
          % call solver function to update weights
          [net.params(p).value, state.solverState{p}] = ...
            params.solver(net.params(p).value, state.solverState{p}, ...
            parDer, params.solverOpts, thisLR) ;
        end
      end
    otherwise
      error('Unknown training method ''%s'' for parameter ''%s''.', ...
        net.params(p).trainMethod, ...
        net.params(p).name) ;
  end
end

% -------------------------------------------------------------------------
function stats = accumulateStats(stats_)
% -------------------------------------------------------------------------

for s = {'train', 'val'}
  s = char(s) ;
  total = 0 ;

  % initialize stats stucture with same fields and same order as
  % stats_{1}
  stats__ = stats_{1} ;
  names = fieldnames(stats__.(s))' ;
  values = zeros(1, numel(names)) ;
  fields = cat(1, names, num2cell(values)) ;
  stats.(s) = struct(fields{:}) ;

  for g = 1:numel(stats_)
    stats__ = stats_{g} ;
    num__ = stats__.(s).num ;
    total = total + num__ ;

    for f = setdiff(fieldnames(stats__.(s))', 'num')
      f = char(f) ;
      stats.(s).(f) = stats.(s).(f) + stats__.(s).(f) * num__ ;

      if g == numel(stats_)
        stats.(s).(f) = stats.(s).(f) / total ;
      end
    end
  end
  stats.(s).num = total ;
end

% -------------------------------------------------------------------------
function stats = extractStats(stats, net)
% -------------------------------------------------------------------------
sel = find(cellfun(@(x) isa(x,'dagnn.Loss'), {net.layers.block})) ;
for i = 1:numel(sel)
  if net.layers(sel(i)).block.ignoreAverage, continue; end;
  stats.(net.layers(sel(i)).outputs{1}) = net.layers(sel(i)).block.average ;
end

% -------------------------------------------------------------------------
function saveState(fileName, net_, state)
% -------------------------------------------------------------------------
net = net_.saveobj() ;
save(fileName, 'net', 'state') ;

% -------------------------------------------------------------------------
function saveStats(fileName, stats)
% -------------------------------------------------------------------------
if exist(fileName)
  save(fileName, 'stats', '-append') ;
else
  save(fileName, 'stats') ;
end

% -------------------------------------------------------------------------
function [net, state, stats] = loadState(fileName)
% -------------------------------------------------------------------------
load(fileName, 'net', 'state', 'stats') ;
net = dagnn.DagNN.loadobj(net) ;
if isempty(whos('stats'))
  error('Epoch ''%s'' was only partially saved. Delete this file and try again.', ...
        fileName) ;
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

% -------------------------------------------------------------------------
function switchFigure(n)
% -------------------------------------------------------------------------
if get(0,'CurrentFigure') ~= n
  try
    set(0,'CurrentFigure',n) ;
  catch
    figure(n) ;
  end
end

% -------------------------------------------------------------------------
function clearMex()
% -------------------------------------------------------------------------
clear vl_tmove vl_imreadjpeg ;

% -------------------------------------------------------------------------
function prepareGPUs(opts, cold)
% -------------------------------------------------------------------------
numGpus = numel(opts.gpus) ;
if numGpus > 1
  % check parallel pool integrity as it could have timed out
  pool = gcp('nocreate') ;
  if ~isempty(pool) && pool.NumWorkers ~= numGpus
    delete(pool) ;
  end
  pool = gcp('nocreate') ;
  if isempty(pool)
    parpool('local', numGpus) ;
    cold = true ;
  end

end
if numGpus >= 1 && cold
  fprintf('%s: resetting GPU\n', mfilename)
  clearMex() ;
  if numGpus == 1
    gpuDevice(opts.gpus)
  else
    spmd
      clearMex() ;
      gpuDevice(opts.gpus(labindex))
    end
  end
end

repo: LapSRN


text:
% ==============================================================================
%
%    Software License Agreement (BSD License)
%    Copyright (c) 2019
%    (www.aimlab.wpi.edu)
%
%    All rights reserved.
%
%    Redistribution and use in source and binary forms, with or without
%    modification, are permitted provided that the following conditions
%    are met:
%
%    * Redistributions of source code must retain the above copyright
%    notice, this list of conditions and the following disclaimer.
%
%    * 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.
%
%    * Neither the name of authors nor the names of its contributors may
%    be used to endorse or promote products derived from this software
%    without specific prior written permission.
%
%    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 OWNER 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.
%
%    \author:    <http://www.aimlab.wpi.edu>
%    \author:    <amunawar@wpi.edu>
%    \author:    Adnan Munawar
%    \version:   0.1$
% ==============================================================================

function [its,sizePath,run_time] =  RRTstar3D(dim,segmentLength,radius,random_world,show_output,samples)
% if samples < 4000
%     disp('ERROR! SPECIFY ATLEAST 4000 SAMPLES')
%     return
% end

% dim = 2;
% radius =0;
% segmentLength = 5;
% random_world = 0;
% n_its = 1000;
% standard length of path segments
if dim ==2
    start_cord = [5,5];
    goal_cord = [95,95];

else

    start_cord = [5,5,5];
    goal_cord = [95,95,95];
end

% create random world
Size = 100;
NumObstacles = 100;

if random_world ==1
    world = createWorld(NumObstacles,ones(1,dim)*Size,zeros(1,dim),dim);
else
    [world NumObstacles] = createKnownWorld(ones(1,dim)*Size,[0;0;0],dim);
end
% randomly select start and end nodes
%start_node = generateRandomNode(world,dim)
%end_node   = generateRandomNode(world,dim)
start_node = [start_cord,0,0,0];
end_node = [goal_cord,0,0,0];
% establish tree starting with the start node
tree = start_node;

numPaths = 0;
a = clock;
% check to see if start_node connects directly to end_node
if ( (norm(start_node(1:dim)-end_node(1:dim))<segmentLength )...
        &&(collision(start_node,end_node,world,dim)==0) )

    disp('CASE 1')
    path = [start_node; end_node];
else
    disp('CASE 2')
    if samples >0
        draw = floor(samples/8);
        its = 0;
        for i = 1:samples
            flag = 0;
            [tree,flag] = extendTree(tree,end_node,segmentLength,radius,world,flag,dim);
            numPaths = numPaths + flag;
            its = its+1;

            if its == draw
                tree_1 = tree;
            elseif its == draw*2
                tree_2 = tree;
            elseif its == draw*3
                tree_3 = tree;
            elseif its == draw*4
                tree_4 = tree;
            elseif its == draw*5
                tree_5 = tree;
            elseif its == draw*6
                tree_6 = tree;
            elseif its == draw*7
                tree_7 = tree;
            elseif its == samples
                tree_8 = tree;
            end
        end

    else
        its = 0;
        numPaths = 0;
        flag = 0;
        while numPaths < 1,
            [tree,flag] = extendTree(tree,end_node,segmentLength,radius,world,flag,dim);
            numPaths = numPaths + flag;
            its = its+1;
        end
    end

end
numPaths

% find path with minimum cost to end_node
path = findMinimumPath(tree,end_node,dim);

b = clock;
run_time = 3600*(b(4)-a(4)) + 60 * (b(5)-a(5)) + (b(6) - a(6));

path_1 = findMinimumPath(tree_1,end_node,dim);

path_2 = findMinimumPath(tree_2,end_node,dim);

path_3 = findMinimumPath(tree_3,end_node,dim);

path_4 = findMinimumPath(tree_4,end_node,dim);

path_5 = findMinimumPath(tree_5,end_node,dim);

path_6 = findMinimumPath(tree_6,end_node,dim);

path_7 = findMinimumPath(tree_7,end_node,dim);

path_8 = findMinimumPath(tree_8,end_node,dim);

sizePath = size(path,1);

if show_output == 1

    if size(path_1, 1) > 0
        figure;
        plotExpandedTree(world,tree_1,dim);
        plotWorld(world,path_1,dim);
    else
        disp('COULD NOT FIND A CONNECTING TREE TILL 1/8th SAMPLES SO NOT DRAWING THAT PATH')
    end
    if size(path_2, 1) > 0
        figure;
        plotExpandedTree(world,tree_2,dim);
        plotWorld(world,path_2,dim);
    else
        disp('COULD NOT FIND A CONNECTING TREE TILL 2/8th SAMPLES SO NOT DRAWING THAT PATH')
    end
    if size(path_3, 1) > 0
    figure;
    plotExpandedTree(world,tree_3,dim);
    plotWorld(world,path_3,dim);
    else
        disp('COULD NOT FIND A CONNECTING TREE TILL 3/8th SAMPLES SO NOT DRAWING THAT PATH')
    end
    if size(path_4, 1) > 0
    figure;
    plotExpandedTree(world,tree_4,dim);
    plotWorld(world,path_4,dim);
    else
        disp('COULD NOT FIND A CONNECTING TREE TILL 4/8th SAMPLES SO NOT DRAWING THAT PATH')
    end
    if size(path_5, 1) > 0
    figure;
    plotExpandedTree(world,tree_5,dim);
    plotWorld(world,path_5,dim);
    else
        disp('COULD NOT FIND A CONNECTING TREE TILL 5/8th SAMPLES SO NOT DRAWING THAT PATH')
    end
    if size(path_6, 1) > 0
    figure;
    plotExpandedTree(world,tree_6,dim);
    plotWorld(world,path_6,dim);
    else
        disp('COULD NOT FIND A CONNECTING TREE TILL 6/8th SAMPLES SO NOT DRAWING THAT PATH')
    end
    if size(path_7, 1) > 0
    figure;
    plotExpandedTree(world,tree_7,dim);
    plotWorld(world,path_7,dim);
    else
        disp('COULD NOT FIND A CONNECTING TREE TILL 7/8th SAMPLES SO NOT DRAWING THAT PATH')
    end
    if size(path_8, 1) > 0
    figure;
    plotExpandedTree(world,tree_8,dim);
    plotWorld(world,path_8,dim);
    else
        disp('COULD NOT FIND A CONNECTING TREE TILL 8/8th SAMPLES SO NOT DRAWING THAT PATH')
    end
    if size(path, 1) > 0
    figure;
    plotExpandedTree(world,tree,dim);
    plotWorld(world,path,dim);
    else
        disp('COULD NOT FIND A CONNECTING TREE FOR THE SPECIFIED SAMPLES. PLEASE INCREASE THE NUMBER OF SAMPLES')
    end
end
end

function world = createWorld(NumObstacles, endcorner, origincorner,dim)

if dim == 2

    % check to make sure that the region is nonempty
    if (endcorner(1) <= origincorner(1)) | (endcorner(2) <= origincorner(2))
        disp('Not valid corner specifications!')
        world=[];

        % create world data structure
    else
        world.NumObstacles = NumObstacles;
        world.endcorner = endcorner;
        world.origincorner = origincorner;

        % create NumObstacles
        maxRadius = min(endcorner(1)- origincorner(1), endcorner(2)-origincorner(2));
        maxRadius = 5*maxRadius/NumObstacles/2;
        for i=1:NumObstacles,
            % randomly pick radius
            world.radius(i) = maxRadius*rand;
            % randomly pick center of obstacles
            cx = origincorner(1) + world.radius(i)...
                + (endcorner(1)-origincorner(1)-2*world.radius(i))*rand;
            cy = origincorner(2) + world.radius(i)...
                + (endcorner(2)-origincorner(2)-2*world.radius(i))*rand;
            world.cx(i) = cx;
            world.cy(i) = cy;
        end
    end

elseif dim ==3;
    % check to make sure that the region is nonempty
    if (endcorner(1) <= origincorner(1)) || (endcorner(2) <= origincorner(2)) || (endcorner(3) <= origincorner(3))
        disp('Not valid corner specifications!')
        world=[];

        % create world data structure
    else
        world.NumObstacles = NumObstacles;
        world.endcorner = endcorner;
        world.origincorner = origincorner;

        % create NumObstacles
        bounds = [endcorner(1)- origincorner(1), endcorner(2)-origincorner(2), endcorner(3)-origincorner(3)];
        maxRadius = min(bounds);
        maxRadius = 5*maxRadius/NumObstacles;
        for i=1:NumObstacles,
            % randomly pick radius
            world.radius(i) = maxRadius*rand;
            % randomly pick center of obstacles
            cx = origincorner(1) + world.radius(i)...
                + (endcorner(1)-origincorner(1)-2*world.radius(i))*rand;
            cy = origincorner(2) + world.radius(i)...
                + (endcorner(2)-origincorner(2)-2*world.radius(i))*rand;
            cz = origincorner(2) + world.radius(i)...
                + (endcorner(2)-origincorner(2)-2*world.radius(i))*rand;
            world.cx(i) = cx;
            world.cy(i) = cy;
            world.cz(i) = cz;
        end
    end
end
end

function [world NumObstacles] = createKnownWorld(endcorner, origincorner,dim)
NumObstacles = 5;
if dim == 2
    % check to make sure that the region is nonempty
    if (endcorner(1) <= origincorner(1)) | (endcorner(2) <= origincorner(2)),
        disp('Not valid corner specifications!')
        world=[];
        % create world data structure
    else
        world.NumObstacles = NumObstacles;
        world.endcorner = endcorner;
        world.origincorner = origincorner;

        % create NumObstacles
        maxRadius = 10;

        world.radius(1) = maxRadius;
        cx = 50;
        cy = 50;
        world.cx(1) = cx;
        world.cy(1) = cy;

        world.radius(2) = maxRadius;
        cx = 75;
        cy = 25;
        world.cx(2) = cx;
        world.cy(2) = cy;

        world.radius(3) = maxRadius;
        cx = 25;
        cy = 75;
        world.cx(3) = cx;
        world.cy(3) = cy;

        world.radius(4) = maxRadius;
        cx = 25;
        cy = 25;
        world.cx(4) = cx;
        world.cy(4) = cy;

        world.radius(5) = maxRadius;
        cx = 75;
        cy = 75;
        world.cx(5) = cx;
        world.cy(5) = cy;
    end

elseif dim == 3

    NumObstacles = 9;
    % check to make sure that the region is nonempty
    if (endcorner(1) <= origincorner(1)) | (endcorner(2) <= origincorner(2)) | (endcorner(3) <= origincorner(3)),
        disp('Not valid corner specifications!')
        world=[];

        % create world data structure
    else
        world.NumObstacles = NumObstacles;
        world.endcorner = endcorner;
        world.origincorner = origincorner;

        % create NumObstacles
        maxRadius = 10;

        world.radius(1) = maxRadius;
        cx = 50;
        cy = 50;
        cz = 50;
        world.cx(1) = cx;
        world.cy(1) = cy;
        world.cz(1) = cz;

        world.radius(2) = maxRadius;
        cx = 25;
        cy = 25;
        cz = 25;
        world.cx(2) = cx;
        world.cy(2) = cy;
        world.cz(2) = cz;

        world.radius(3) = maxRadius;
        cx = 75;
        cy = 75;
        cz = 75;
        world.cx(3) = cx;
        world.cy(3) = cy;
        world.cz(3) = cz;

        world.radius(4) = maxRadius;
        cx = 25;
        cy = 25;
        cz = 75;
        world.cx(4) = cx;
        world.cy(4) = cy;
        world.cz(4) = cz;

        world.radius(5) = maxRadius;
        cx = 75;
        cy = 75;
        cz = 25;
        world.cx(5) = cx;
        world.cy(5) = cy;
        world.cz(5) = cz;

        world.radius(6) = maxRadius;
        cx = 25;
        cy = 75;
        cz = 25;
        world.cx(6) = cx;
        world.cy(6) = cy;
        world.cz(6) = cz;

        world.radius(7) = maxRadius;
        cx = 75;
        cy = 25;
        cz = 25;
        world.cx(7) = cx;
        world.cy(7) = cy;
        world.cz(7) = cz;

        world.radius(8) = maxRadius;
        cx = 75;
        cy = 25;
        cz = 75;
        world.cx(8) = cx;
        world.cy(8) = cy;
        world.cz(8) = cz;

        world.radius(9) = maxRadius;
        cx = 25;
        cy = 75;
        cz = 75;
        world.cx(9) = cx;
        world.cy(9) = cy;
        world.cz(9) = cz;
    end
end
end

function node=generateRandomNode(world,dim)

if dim ==2;
    % randomly pick configuration
    px       = (world.endcorner(1)-world.origincorner(1))*rand;
    py       = (world.endcorner(2)-world.origincorner(2))*rand;

    chi      = 0;
    cost     = 0;
    node     = [px, py, chi, cost, 0];

    % check collision with obstacle
    while collision(node, node, world,dim),
        px       = (world.endcorner(1)-world.origincorner(1))*rand;
        py       = (world.endcorner(2)-world.origincorner(2))*rand;

        chi      = 0;
        cost     = 0;
        node     = [px, py, chi, cost, 0];
    end

elseif dim ==3;
    % randomly pick configuration
    px       = (world.endcorner(1)-world.origincorner(1))*rand;
    py       = (world.endcorner(2)-world.origincorner(2))*rand;
    pz       = (world.endcorner(3)-world.origincorner(3))*rand;

    chi      = 0;
    cost     = 0;
    node     = [px, py, pz, chi, cost, 0];

    % check collision with obstacle
    while collision(node, node, world,dim),
        px       = (world.endcorner(1)-world.origincorner(1))*rand;
        py       = (world.endcorner(2)-world.origincorner(2))*rand;
        pz       = (world.endcorner(3)-world.origincorner(3))*rand;

        chi      = 0;
        cost     = 0;
        node     = [px, py, pz, chi, cost, 0];
    end

end

end

function collision_flag = collision(node, parent, world,dim)

collision_flag = 0;

for i=1:dim
    if (node(i)>world.endcorner(i))|(node(i)<world.origincorner(i))
        collision_flag = 1;
    end
end

if collision_flag == 0 && dim ==2
    for sigma = 0:.2:1,
        p = sigma*node(1:dim) + (1-sigma)*parent(1:dim);
        % check each obstacle
        for i=1:world.NumObstacles,
            if (norm([p(1);p(2)]-[world.cx(i); world.cy(i)])<=1*world.radius(i)),
                collision_flag = 1;
                break;
            end
        end
    end

elseif collision_flag == 0 && dim ==3
    for sigma = 0:.2:1,
        p = sigma*node(1:dim) + (1-sigma)*parent(1:dim);
        % check each obstacle
        for i=1:world.NumObstacles,
            if (norm([p(1);p(2);p(3)]-[world.cx(i); world.cy(i); world.cz(i)])<=1*world.radius(i)),
                collision_flag = 1;
                break;
            end
        end
    end
end
end

function flag = canEndConnectToTree(tree,end_node,minDist,world,dim)
flag = 0;
% check only last node added to tree since others have been checked
if ( (norm(tree(end,1:dim)-end_node(1:dim))<minDist)...
        & (collision(tree(end,1:dim), end_node(1:dim), world,dim)==0) ),
    flag = 1;
end

end

function [new_tree,flag] = extendTree(tree,end_node,segmentLength,r,world,flag_chk,dim)

flag1 = 0;
while flag1==0,
    % select a random point
    randomPoint = ones(1,dim);
    for i=1:dim
        randomPoint(1,i) = (world.endcorner(i)-world.origincorner(i))*rand;
    end

    % find leaf on node that is closest to randomPoint
    tmp = tree(:,1:dim)-ones(size(tree,1),1)*randomPoint;
    sqrd_dist = sqr_eucl_dist(tmp,dim);
    [min_dist,idx] = min(sqrd_dist);
    min_parent_idx = idx;

    new_point = (randomPoint-tree(idx,1:dim));
    new_point = tree(idx,1:dim)+(new_point/norm(new_point))*segmentLength;

    min_cost  = cost_np(tree(idx,:),new_point,dim);
    new_node  = [new_point, 0, min_cost, idx];

    if collision(new_node, tree(idx,:), world,dim)==0

        tmp_dist = tree(:,1:dim)-(ones(size(tree,1),1)*new_point);
        dist = sqr_eucl_dist(tmp_dist,dim);
        near_idx = find(dist <= r^2);

        if size(near_idx,1)>1
            size_near = size(near_idx,1);

            for i = 1:size_near
                if collision(new_node, tree(near_idx(i),:), world,dim)==0

                    cost_near = tree(near_idx(i),dim+2)+line_cost(tree(near_idx(i),:),new_point,dim);

                    if  cost_near < min_cost
                        min_cost = cost_near;
                        min_parent_idx = near_idx(i);
                    end

                end
            end
        end

        new_node = [new_point, 0 , min_cost, min_parent_idx];
        new_tree = [tree; new_node];
        new_node_idx = size(new_tree,1);

        if size(near_idx,1)>1
            reduced_idx = near_idx;
            for j = 1:size(reduced_idx,1)
                near_cost = new_tree(reduced_idx(j),dim+2);
                lcost = line_cost(new_tree(reduced_idx(j),:),new_point,dim);
                if near_cost > min_cost + lcost ...
                        && collision(new_tree(reduced_idx(j),:),new_node,world,dim)
                    before = new_tree(reduced_idx(j),dim+3)
                    new_tree(reduced_idx(j),dim+3) = new_node_idx;
                    after = new_tree(reduced_idx(j),dim+3)
                end

            end
        end
        flag1=1;
    end
end

if flag_chk == 0
    % check to see if new node connects directly to end_node
    if ( (norm(new_node(1:dim)-end_node(1:dim))<segmentLength )...
            && (collision(new_node,end_node,world,dim)==0) )
        flag = 1;
        new_tree(end,dim+1)=1;  % mark node as connecting to end.
    else
        flag = 0;
    end

else flag = 1;
end
end

function e_dist = sqr_eucl_dist(array,dim)

sqr_e_dist = zeros(size(array,1),dim);
for i=1:dim

    sqr_e_dist(:,i) = array(:,i).*array(:,i);

end
e_dist = zeros(size(array,1),1);
for i=1:dim

    e_dist = e_dist+sqr_e_dist(:,i);

end

end

%calculate the cost from a node to a point
function [cost] = cost_np(from_node,to_point,dim)

diff = from_node(:,1:dim) - to_point;
eucl_dist = norm(diff);
cost = from_node(:,dim+2) + eucl_dist;

end

%calculate the cost from a node to a node
function [cost] = cost_nn(from_node,to_node,dim)

diff = from_node(:,1:dim) - to_node(:,1:dim);
eucl_dist = norm(diff);
cost = from_node(:,dim+2) + eucl_dist;

end

function [cost] = line_cost(from_node,to_point,dim)
diff = from_node(:,1:dim) - to_point;
cost = norm(diff);
end

function path = findMinimumPath(tree,end_node,dim)

% find nodes that connect to end_node
connectingNodes = [];
for i=1:size(tree,1),
    if tree(i,dim+1)==1,
        connectingNodes = [connectingNodes ; tree(i,:)];
    end
end

if size(connectingNodes, 1) > 0

    % find minimum cost last node
    [tmp,idx] = min(connectingNodes(:,dim+2));

    % construct lowest cost path
    path = [connectingNodes(idx,:); end_node];
    parent_node = connectingNodes(idx,dim+3);
    while parent_node>1,
        parent_node = tree(parent_node,dim+3);
        path = [tree(parent_node,:); path];
    end

else
    path = [];
end

end

function plotExpandedTree(world,tree,dim)
ind = size(tree,1);
while ind>0
    branch = [];
    node = tree(ind,:);
    branch = [ branch ; node ];
    parent_node = node(dim+3);
    while parent_node > 1
        cur_parent = parent_node;
        branch = [branch; tree(parent_node,:)];
        parent_node = tree(parent_node,dim+3);
    end
    ind = ind - 1;

    if dim == 2
        X = branch(:,1);
        Y = branch(:,2);

        p = plot(X,Y);
        set(p,'Color','r','LineWidth',0.5,'Marker','.','MarkerEdgeColor','g');
        hold on;

    elseif dim == 3
        X = branch(:,1);
        Y = branch(:,2);
        Z = branch(:,3);

        p = plot3(X,Y,Z);
        set(p,'Color','r','LineWidth',0.5,'Marker','.','MarkerEdgeColor','g');
        hold on;
    end
end
end

function plotWorld(world,path,dim)
% the first element is the north coordinate
% the second element is the south coordinate
if dim ==2

    N = 10;
    th = 0:2*pi/N:2*pi;
    axis([world.origincorner(1),world.endcorner(1),...
        world.origincorner(2), world.endcorner(2)]);
    hold on

    for i=1:world.NumObstacles,
        X = world.radius(i)*sin(th) + world.cx(i);
        Y = world.radius(i)*cos(th) + world.cy(i);
        fill(X,Y,'blue');
    end

    X = path(:,1);
    Y = path(:,2);
    p = plot(X,Y);

elseif dim ==3
    axis([world.origincorner(1),world.endcorner(1),...
        world.origincorner(2), world.endcorner(2),...
        world.origincorner(3), world.endcorner(3)]);
    hold on

    for i=1:world.NumObstacles,
        [X Y Z] = sphere(10);
        X = (X*world.radius(i));
        Y = (Y*world.radius(i));
        Z = (Z*world.radius(i));
        surf(X+world.cx(i),Y+world.cy(i),Z+world.cz(i));
        colormap([0.5 0.2 0.3]);
    end

    X = path(:,1);
    Y = path(:,2);
    Z = path(:,3);
    p = plot3(X,Y,Z);
end
set(p,'Color','black','LineWidth',3)
xlabel('X axis');
ylabel('Y axis');
zlabel('Z axis');
title('RRT Star Algorithm');
end

author: adnanmunawar

repo: matlab-rrt-variants


text:
function varargout = process_fix_headcoils(varargin)
% PROCESS_FIX_HEADCOILS: Estimate the position of a missing or bad head localization coil.

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

eval(macro_method);
end

function sProcess = GetDescription() %#ok<DEFNU>
    % Description of the process
    sProcess.Comment     = 'Fix bad head coil (CTF)';
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/HeadMotion#Dealing_with_a_missing_or_bad_head_coil';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = {'Import', 'Channel file'};
    sProcess.Index       = 54;
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw', 'data'};
    sProcess.OutputTypes = {'raw', 'data'};
    sProcess.nInputs     = 2;
    sProcess.nOutputs    = 1;
    sProcess.nMinFiles   = 1;
    sProcess.isSeparator = 0;
    sProcess.isPaired    = 1;
    sProcess.FileTag     = 'fix';

    % Options
    % === Sensor types
    sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
    sProcess.options.sensortypes.Type    = 'text';
    sProcess.options.sensortypes.Value   = {'HLU'};
    sProcess.options.sensortypes.Hidden  = 1;
    %     sProcess.options.read_all.Comment    = ''; % Automatically added to filter-type processes.
    %     sProcess.options.read_all.Type       = 'checkbox';
    %     sProcess.options.read_all.Value      = 1;
    %     sProcess.options.read_all.Hidden     = 1;

    sProcess.processDim = 2; % Only split by column blocks, if needed.

    sProcess.options.BadCoil.Type    = 'radio_line';
    sProcess.options.BadCoil.Comment = {'Nasion', 'Left ear', 'Right ear', 'Bad coil: '};
    sProcess.options.BadCoil.Value   = 1;
    %     sProcess.options.label1.Type    = 'label';
    %     sProcess.options.label1.Comment = 'How will you provide a "good and similar" head position?';
    %     sProcess.options.GoodSource.Type    = 'radio';
    %     sProcess.options.GoodSource.Comment = {'Recording already in Brainstorm', 'Raw recording file', 'Manual entry of dewar coordinates'};
    %     sProcess.options.GoodSource.Value   = 1;
    %     sProcess.options.display.Type    = 'checkbox';
    %     sProcess.options.display.Comment = 'Display "before" and "after" alignment figures.';
    %     sProcess.options.display.Value   = 0;

end

function Comment = FormatComment(sProcess)
    Comment = [sProcess.Comment, ': ', ...
        sProcess.options.BadCoil.Comment{sProcess.options.BadCoil.Value}];

end

function OutputFiles = Run(sProcess, sInputA, sInputB)

    % If called from ProcessFilter
    if isfield(sProcess, 'FromProcessFilter') && sProcess.FromProcessFilter
            OutputFiles = [];
            ChannelMat = in_bst_channel(sInputA.ChannelFile);

           % Check the input is CTF.
            DataMat = in_bst_data(sInputA.FileName, 'Device');
            if ~strcmp(DataMat.Device, 'CTF')
                bst_report('Error', sProcess, sInputA, ...
                    'Fix bad head coil is currently only available for CTF data.');
                return;
            end

            % Difficult to find the initial head sample in each epoch. And
            % if wrong, fixed HLU channel will have offset "steps". So
            % process redundant samples, which we would need to recreate
            % anyway.
            %             % Load head coil locations, in m.
            %             bst_progress('text', 'Loading head coil locations...');
            %             [Locations, HeadSamplePeriod] = process_evt_head_motion('LoadHLU', sInputA, [], false); % [nChannels, nSamples, nEpochs]
            %             if isempty(Locations)
            %                 % No HLU channels. Error already reported. Skip this file.
            %                 return;
            %             end
            bst_progress('text', 'Correcting head position...');
            % Sort in case coils were swapped.
            [Unused, iSortHlu] = sort({ChannelMat.Channel(strcmp({ChannelMat.Channel.Type}, 'HLU')).Name});

            % Compute
            GoodAxis = sInputA.A(iSortHlu(sProcess.Ext.iGood(1:3)), :, :) - sInputA.A(iSortHlu(sProcess.Ext.iGood(4:6)), :, :); %[3, nS, nE]
            GoodAxis = bsxfun(@rdivide, GoodAxis, sqrt(sum(GoodAxis.^2, 1)));
            FixedCoil = (sInputA.A(iSortHlu(sProcess.Ext.iGood(4:6)), :, :) + sInputA.A(iSortHlu(sProcess.Ext.iGood(1:3)), :, :)) ./ 2 + ...
                GoodAxis * sProcess.Ext.BadAlongGoodAxis + bsxfun(@cross, GoodAxis, sProcess.Ext.Normal) * sProcess.Ext.BadFromGoodAxis;
            % If a collection was aborted, the channels will be filled with
            % zeros. The "fixed" data will contain NaN, change them back to 0.
            FixedCoil(isnan(FixedCoil)) = 0;

            %             % Re-upsample (duplicate) and save new data.
            %             [n3, nHeadSamples, nEpochs] = size(FixedCoil);
            %             FixedCoil = reshape(permute(FixedCoil, [2, 3, 1]), [], 1);
            %             FixedCoil = permute(reshape(FixedCoil(:, ones(1,HeadSamplePeriod))', [], nEpochs, n3), [3, 1, 2]);

            %             [nChan, nSamples, nEpochs] = size(sInputA.A);
            %             sInputA.A(iSortHlu(sProcess.Ext.iBad), :, :) = FixedCoil(:, 1:nSamples, 1:nEpochs);
            sInputA.A(iSortHlu(sProcess.Ext.iBad), :, :) = FixedCoil;
            OutputFiles = sInputA;

    % Initial call from pipeline
    else
        OutputFiles = {sInputA.FileName};
        % Flag to indicate Run will be called by ProcessFilter.
        sProcess.FromProcessFilter = 1;

        % Bad coil
        sProcess.Ext.iBad = (1:3) + (sProcess.options.BadCoil.Value - 1) * 3;
        % Good coils
        sProcess.Ext.iGood = setdiff(1:9, sProcess.Ext.iBad);

        bst_progress('start', 'Fix bad head coil', ' ', 0, numel(sInputA));

        % Process each input file pair.
        for iFile = 1:length(sInputA)
            % [Possibly to add: options to load from external file or prompt for manual entry, showing current (bad) dewar coordinates.]
            % Get good head coil Dewar coordinates from second channel file.
            ChannelMatB = in_bst_channel(sInputB(iFile).ChannelFile);
            sProcess.Ext.RefLoc = process_adjust_coordinates('ReferenceHeadLocation', ChannelMatB);
            % Sort in case coils were swapped.
            [Unused, iSortHlu] = sort({ChannelMatB.Channel(strcmp({ChannelMatB.Channel.Type}, 'HLU')).Name});

            % Get relative coordinates of the bad coil from the external source.
            sProcess.Ext.CoilA = sProcess.Ext.RefLoc(iSortHlu(sProcess.Ext.iGood(1:3)));
            sProcess.Ext.CoilB = sProcess.Ext.RefLoc(iSortHlu(sProcess.Ext.iGood(4:6)));
            sProcess.Ext.Orig = (sProcess.Ext.CoilA + sProcess.Ext.CoilB) / 2;
            sProcess.Ext.GoodAxis = sProcess.Ext.CoilA - sProcess.Ext.CoilB;
            sProcess.Ext.GoodAxis = sProcess.Ext.GoodAxis ./ norm(sProcess.Ext.GoodAxis);
            sProcess.Ext.BadFromO = sProcess.Ext.RefLoc(iSortHlu(sProcess.Ext.iBad)) - sProcess.Ext.Orig;
            sProcess.Ext.BadAlongGoodAxis = sProcess.Ext.BadFromO' * sProcess.Ext.GoodAxis;
            sProcess.Ext.BadFromGoodAxis = norm(sProcess.Ext.BadFromO - sProcess.Ext.BadAlongGoodAxis * sProcess.Ext.GoodAxis);
            sProcess.Ext.Normal = cross(sProcess.Ext.BadFromO, sProcess.Ext.GoodAxis);
            sProcess.Ext.Normal = sProcess.Ext.Normal ./ norm(sProcess.Ext.Normal);
            % Capture process crashes
            try
                OutputFiles{iFile} = bst_process('ProcessFilter', sProcess, sInputA(iFile));
                bst_progress('inc', 1);
            catch
                strError = bst_error();
                bst_report('Error', sProcess, [sInputA(iFile), sInputB(iFile)], strError);
                continue;
            end

        end
        bst_progress('stop');

        % Attempt to adjust the initial/reference head position.
        %             ChannelMat = in_bst_channel(sInputA(iFile).ChannelFile);
        %             sOutput = in_bst_data(
        OutputFiles = bst_process('CallProcess', 'process_adjust_coordinates', OutputFiles, [], ...
            'reset', 0, 'head', 1, 'bad', 1, 'points', 0, 'remove', 0, 'display', 1);
        %             [ChannelMat, Failed] = process_adjust_coordinates('AdjustHeadPosition', ...
        %                 ChannelMat, OutputFiles(iFile), sProcess);
        %             if ~Failed
        %                 % Save channel file.
        %                 bst_save(file_fullpath(OutputFiles(iFile).ChannelFile), ChannelMat, 'v7');
        %                 % Show new alignment.
        %                 channel_align_manual(OutputFiles(iFile).ChannelFile, 'MEG', 0);
        %             % else already noted in report.
        %             end
    end % "called from" if
end

author: brainstorm-tools


text:
function cc_levels_minmax_animate ( )

%*****************************************************************************80
%
%% CC_LEVELS_MINMAX_ANIMATE displays the sequence of grids generated by CC_LEVELS_MINMAX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 November 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_LEVELS_MINMAX_ANIMATE:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Display the sequence of nested 2D Clenshaw-Curtis grids\n' );
  fprintf ( 1, '  generated by CC_LEVELS_MINMAX.\n' );

  dim_num = 2;

  while ( 1 )
%
%  Get user input.
%
    fprintf ( 1, '\n' );
    level_max = input ( 'Enter LEVEL_MAX or RETURN to exit;' );

    if ( isempty ( level_max ) )
      break
    end

    grid_points_old = [];
    grid_points_old_num = 0;

    grid_points = [];
    grid_points_num = 0;

    for level = 0 : level_max

      if ( 0 < grid_points_num )

        grid_points_old(1:2,grid_points_old_num+1:grid_points_old_num+grid_points_num) = ...
          grid_points(1:2,1:grid_points_num);

        grid_points_old_num = grid_points_old_num + grid_points_num;

      end
%
%  Compute data.
%
      [ grid_num, point_num ] = cc_levels_minmax_size ( dim_num, ...
        level, level );

      [ grid_level, grid_order, grid_points ] = cc_levels_minmax ( dim_num, ...
        level, level,  grid_num, point_num );

      [ dim_num, grid_points_num ] = size ( grid_points );

      clf
%
%  We have to name the axes in order to control the grid.
%
      axes_handle = axes;
%
%  Plot the points.
%
      handle_new = scatter ( grid_points(1,:), grid_points(2,:), 'r', 'filled' );

      hold on

      if ( 0 < grid_points_old_num )
        handle_old = scatter ( grid_points_old(1,:), grid_points_old(2,:), 'b', 'filled' );
      end
%
%  Force the plotting region to be square, not rectangular.
%
      axis square
%
%  Request grid lines.
%
      grid on
%
%  Specify the location of the grid lines, and suppress labeling.
%
      set ( axes_handle, 'xtick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
      set ( axes_handle, 'xticklabel', [] );
      set ( axes_handle, 'ytick', [ -1, -.75, -.5, -.25, 0, .25, .50, .75, 1] );
      set ( axes_handle, 'yticklabel', [] );
%
%  Make the plotting region slightly bigger than the data.
%
      axis ( [ -1.1, 1.1, -1.1, 1.1 ] )
%
%  Title
%
      s = sprintf ( '+LEVEL %d', level );
      title ( s );

      fprintf ( 1, '+LEVEL %d, Press return\n', level );
      pause

      hold off

    end

    fprintf ( 1, 'Press return to CLEAR the grid!\n' );
    pause
    clf

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CC_LEVELS_MINMAX_ANIMATE:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  return
end