text:
function scimatsk = scimat_skeleton_prune(scimatsk, scimat, minlen, lratio, alphamax, p)
% SCIMAT_SKELETON_PRUNE  Prune branches in a segmentation's skeletonization
%
% This function prunes the leaves of a segmentation skeleton in three
% steps:
%
%   1. Erode clumps of voxels to look like branches
%
%   2. Prune very short leaves (iteratively until no more leaves can be
%      pruned)
%
%   3. Prune leaves created by artifacts in the segmentation. Spurious
%      leaves are those roughly as long as the local radius of the main
%      branch they are attached to (iteratively until no more leaves can
%      be pruned)
%
%
% SCIMATPR = scimat_skeleton_prune(SCIMATSK)
% SCIMATPR = scimat_skeleton_prune(SCIMATSK, [], MINLEN)
%
%   This syntax runs steps 1 and 2 only.
%
%   SCIMATSK is an SCIMAT struct. SCIMATSK.data contains the result of
%   running a skeletonization algorithm on a binary segmentation SCIMAT,
%   e.g.
%
%     >> scimatsk = scimat;
%     >> scimatsk.data = itk_imfilter('skel', scimat);
%
%   MINLEN is a scalar with the minimum length in voxels for a leaf. Any
%   leaf shorter than MINLEN will be pruned. By default, MINLEN = 5 voxel.
%
%
% SCIMATPR = scimat_skeleton_prune(SCIMATSK, SCIMAT, MAXCLUMP, MINLEN, LRATIO, ALPHAMAX, P)
%
%   This syntax runs steps 1, 2 and 3.
%
%   SCIMAT is the binary segmentation mentioned above.
%
%   LRATIO is a scalar. Leaves with L/R < LRATIO will be pruned, where
%   L is the length from the bifurcation to the tip of the leaf. By
%   default, LRATIO=1.2.
%
%   ALPHAMAX and P are merging parameters. See the help of function
%   skeleton_label for details. If merging is not enabled, then no leaves
%   will be pruned in step 3.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2011, 2014 University of Oxford
% Version: 0.5.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
narginchk(1, 6);
nargoutchk(0, 1);

if ~isstruct(scimatsk)
    error('SCIMATSK must be an SCIMAT struct')
end

% defaults
if (nargin < 2)
    scimat = [];
end
if (nargin < 3 || isempty(minlen))
    minlen = 5;
end
if (nargin < 4 || isempty(lratio))
    lratio = 1.2;
end
if (nargin < 5 || isempty(alphamax))
    alphamax = -Inf; % don't merge by default
end
if (nargin < 6 || isempty(p))
    p = 1.0; % don't smooth by default
end

%% Step 1: removal of big clumps of voxels

% if we have no voxels to prune, don't waste time processing
if (nnz(scimatsk.data) == 0)
    return
end

% get sparse matrix of distances between voxels. To label the skeleton we
% don't care about the actual distances, just need to know which voxels are
% connected to others. Actual distances are needed to parameterize the
% branches, though
[dsk, dictsk, idictsk] = seg2dmat(scimatsk.data, 'seg', ...
    [scimatsk.axis.spacing]);

% find bifurcation voxels

% compute degree of each skeleton voxel
deg = sum(dsk > 0, 2);

% get distance matrix index of the bifurcation voxels
bifidx = deg >= 3;

% matrix index => image index
bifidx = idictsk(bifidx);

% label connected components of skeleton branches

% remove bifurcation voxels from original skeleton
scimatsk.data(bifidx) = 0;

% get connected components in the image
cc = bwconncomp(scimatsk.data);

% make size vector always have size(3)
if length(cc.ImageSize) == 2
    cc.ImageSize = [cc.ImageSize 1];
end

% connectivity between branches and bifurcation clumps,
% and find which branches should be merged together

% tag each branch with its label
scimatsk.data = labelmatrix(cc);

% create empty image volume and add only bifurcation voxels
sk2 = zeros(cc.ImageSize, 'uint8');
sk2(bifidx) = 1;

% compute connected components to obtain clumps of bifurcation voxels
bifcc = bwconncomp(sk2);

% loop all the bifurcation clumps, to find which branches are neighbours of
% each other
for I = 1:bifcc.NumObjects

    %% find which branches are connected to which bifurcation clumps

    % get voxels in the bifurcation clump
    % linear index -> row, col, slice
    [r, c, s] = ind2sub(cc.ImageSize, bifcc.PixelIdxList{I});

    % get a box 1 voxel bigger than the clump
    r0 = max(1, min(r) - 1);
    c0 = max(1, min(c) - 1);
    s0 = max(1, min(s) - 1);
    rend = min(cc.ImageSize(1), max(r) + 1);
    cend = min(cc.ImageSize(2), max(c) + 1);
    send = min(cc.ImageSize(3), max(s) + 1);

    % extract that box from the volume with the branches
    boxbr = scimatsk.data(r0:rend, c0:cend, s0:send);

    % create a box for the bifurcation clump
    boxbif = zeros(size(boxbr), class(sk2));
    r = r - r0 + 1;
    c = c - c0 + 1;
    s = s - s0 + 1;
    boxbif(sub2ind(size(boxbif), r, c, s)) = 1;

    % create a box with the same size
    box  = zeros(size(boxbif));

    % tag as TODO=2 branch voxels
    box(boxbr ~= 0) = 2;

    % add to the vox the bifurcation clump voxels
    box(boxbif == 1) = 1;

    % dilate the clump 1 voxel
    box = bwregiongrow(box, 2, [], 1);

    % get the branches that the dilated bifurcation clump overlaps
    idx = double(boxbr(box == 1));

    % because in sk bifurcation voxels were removed, we can have "0" values
    % in idx. Remove them
    idx = idx(idx ~= 0);

    % number of branches
    N = length(idx);

    % in order to merge, the bifurcation needs to have at least 2
    % branches
    if (N < 2)
        continue
    end

    % get all combinations of pairs of branches that share this clump
    idx = nchoosek(idx, 2);

    % get voxel indices for bifurcation clump
    bif = bifcc.PixelIdxList{I};

    % duplicate the list. bif0 is going to be the list of superfluous
    % bifurcation voxels
    bif0 = bif;

    % loop each pair of branches combination
    for J = 1:size(idx, 1)

        % get voxel indices for each branch
        br0 = cc.PixelIdxList{idx(J, 1)};
        br1 = cc.PixelIdxList{idx(J, 2)};

        % merge and sort both branches and the bifurcation clump
        br = sort_branch([br0(:); bif(:); br1(:)], ...
            dsk, dictsk, idictsk);

        % remove from the list of superfluous bifurcation voxels those
        % that are needed to connect the two branches
        bif0 = setdiff(bif0, br);

    end

    % remove from the skeleton the list of superfluous bifurcation
    % voxels
    bifcc.PixelIdxList{I} = setdiff(bif, bif0);

end

% convert labelled skeleton to binary mask
scimatsk.data = uint8(scimatsk.data ~= 0);

% add bifurcation voxels to skeleton binary mask
scimatsk.data(cat(1, bifcc.PixelIdxList{:})) = 1;

%% Step 2: pruning of very short leaf branches
while (1)

    % compute skeleton labelling
    [~, cc] = skeleton_label(scimatsk.data, [], [scimatsk.axis.spacing]);

    % get number of voxels in each branch
    n = cellfun(@(x) length(x), cc.PixelIdxList);

    % find leaf-branches that are shorter than the minimum length
    idx1 = find(n < minlen & cc.IsLeaf);

    % remove short branches from the segmentation
    scimatsk.data(cat(1, cc.PixelIdxList{idx1})) = 0;

    % recompute the skeleton labelling
    [~, ~, bifcc, mcon] = skeleton_label(scimatsk.data, [], [scimatsk.axis.spacing]);

    % find bifurcation clusters that are connected to 0 or 1 branches
    idx2 = find(sum(mcon, 1) < 2);

    % remove those bifurcation clumps, because they are not connecting
    % branches, they are either floating alone in space, or terminating a
    % branch
    scimatsk.data(cat(1, bifcc.PixelIdxList{idx2})) = 0;

    % if no bifurcation clumps were found to be removed, stop the
    % algorithm, because that means that no new short leaf-braches can be
    % found either
    if (isempty(idx2))
        break;
    end

end

%% Step 3: pruning of leaf branches produced by segmentation artifacts

% skip if no full segmentation is provided
if isempty(scimat)
    return
end

% repeat the process until no branches are removed
atleastonepruning = true;
while (atleastonepruning)

    % if there are no prunings in this iterations, we stop
    atleastonepruning = false;

    % label the segmentation using multiple merging at every bifurcation
    % clump
    [scimat.data, cc, ~, mcon, ~, cc2] = ...
        skeleton_label(scimatsk.data, scimat.data, [scimat.axis.spacing], ...
        alphamax, p, false);

    % measure the stats of every merged branch
    stats2 = scimat_seg2label_stats(scimat, cc2, p);

    % loop each branch in the list of merged branches
    for I = 1:cc2.NumObjects

        % we consider the current branch a main branch, and any branch
        % coming out of it, a secondary branch

        % compute the major radius of the main branch
        r = sqrt(4 * stats2.Var(2, I));

        % loop consecutive pairs of segments in the main branch
        for J = 1:length(cc2.MergedBranches{I})-1

            % get the bifurcation clump between both segments
            bif = cc2.MergedBifClumps{I}(J);

            % get secondary branches attached to the bifurcation clump
            % (usually, there's only 1 secondary branch, but there can be more)
            bn = setdiff(find(mcon(:, bif)), cc2.MergedBranches{I}(J:J+1));

            % keep only secondary branches that are leaf-branches
            bn = bn(cc.IsLeaf(bn));

            % skip to next segment if there are no valid leaf secondary
            % branches
            if isempty(bn)
                continue
            end

            % compute distance between the first and last voxels in the
            % secondary branch, which gives a better estimate of whether the
            % branch protudes from the main vessel or not, than the branches
            % length (because the branch can be bent)
            len = zeros(size(bn));
            for K = 1:length(len)
                idx = cc.PixelIdxList{bn(K)}([1 end]);
                [r1, c1, s1] = ind2sub(size(scimatsk.data), idx);
                xyz = scimat_index2world([r1(:), c1(:), s1(:)], scimatsk);
                len(K) = sqrt(sum(diff(xyz).^2));
            end

            % get branches that have to be pruned
            toremove = (len / r) <= lratio;
            atleastonepruning = atleastonepruning || any(toremove);

            % if length of the secondary branch is not much larger than the
            % radius of the main branch, we assume that the secondary branch is
            % a segmentation artifact, and remove it
            scimatsk.data(cat(1, ...
                cc.PixelIdxList{bn(toremove)})) = 0;

        end

    end

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% auxiliary functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% sort a set of voxels so that they form a branch
% idx: voxel indices
% d: distance matrix
function [idx, param] = sort_branch(idx, d, dict, idict)

% number of voxels in the current branch
N = length(idx);

% degenerate case in which the branch has only one voxel
if (N == 1)
    param = 0;
    return
end

% create a "branch distance matrix" for only the voxels in the branch
idx2 = full(dict(idx));
dbr = d(idx2, idx2);

% compute Dijkstra's shortest path from an arbitrary voxel in the
% branch to every other voxel
[dp, ~] = dijkstra(dbr, 1);

% the furthest voxel should be one of the ends of the branch
[~, v0] = max(dp);

% compute shortest path to all other voxels in the branch
%
% note that even for apparently "wire" branches, sometimes we get small
% cycles of 1 voxel, and instead of a "linear" shortest-path, we have a
% tree, so in the next step, some of the voxels are going to be thrown
% away
[dp, parents] = dijkstra(dbr, v0);
[~, v1] = max(dp);

% backtrack the whole branch in order from the furthest point to
% the original extreme point
idx = [];
param = [];
J = 1;
while (v1 ~= 0)
    idx(J) = idict(idx2(v1));
    param(J) = dp(v1);
    v1 = parents(v1);
    J = J + 1;
end

% make sure that cc.PixelIdxList{I} is a column vector
idx = idx(:);

% reorder voxels so that the parameterization increases
% monotonically
idx = idx(end:-1:1);
param = param(end:-1:1);

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
% Wrapper for vl_nnlogScale2Scale block
% inputs{1} :   s     : 1 x 1 x 1 x b
% outputs{1}:   y     : 1 x 1 x 1 x b

classdef logScale2Scale < dagnn.Layer
    methods
        function outputs = forward(~, inputs, ~)

            useGPU = isa(inputs{1}, 'gpuArray');
            if useGPU
                outputs{1} = gpuArray( vl_nnlogScale2Scale(gather(inputs{1})) );
            else
                outputs{1} = vl_nnlogScale2Scale(inputs{1});
            end

        end

        function [derInputs, derParams] = backward(~, inputs, ~, derOutputs)

            useGPU = isa(inputs{1}, 'gpuArray');
            if useGPU
                derInputs{1} = gpuArray( vl_nnlogScale2Scale(gather(inputs{1}), gather(derOutputs{1})) );
            else
                derInputs{1} = vl_nnlogScale2Scale(inputs{1}, derOutputs{1});
            end

            derParams = {};
        end

        function outputSizes = getOutputSizes(~, inputSizes)
            outputSizes = inputSizes{1};
        end

        function obj = logScale2Scale(varargin)
            obj.load(varargin);
        end
    end
end

author: anilbas


text:
function [x, info] = umfpack_btf (A, b, Control)
%UMFPACK_BTF factorize A using a block triangular form
%
% Example:
%   x = umfpack_btf (A, b, Control)
%
% solve Ax=b by first permuting the matrix A to block triangular form via dmperm
% and then using UMFPACK to factorize each diagonal block.  Adjacent 1-by-1
% blocks are merged into a single upper triangular block, and solved via
% MATLAB's \ operator.  The Control parameter is optional (Type umfpack_details
% and umfpack_report for details on its use).  A must be square.
%
% See also umfpack, umfpack2, umfpack_details, dmperm

% Copyright 1995-2007 by Timothy A. Davis.

if (nargin < 2)
    help umfpack_btf
    error ('Usage: x = umfpack_btf (A, b, Control)') ;
end

[m n] = size (A) ;
if (m ~= n)
    help umfpack_btf
    error ('umfpack_btf:  A must be square') ;
end
m1 = size (b,1) ;
if (m1 ~= n)
    help umfpack_btf
    error ('umfpack_btf:  b has the wrong dimensions') ;
end

if (nargin < 3)
    Control = umfpack2 ;
end

%-------------------------------------------------------------------------------
% find the block triangular form
%-------------------------------------------------------------------------------

% dmperm built-in may segfault in MATLAB 7.4 or earlier; fixed in MATLAB 7.5
% since dmperm now uses CSparse
[p,q,r] = dmperm (A) ;
nblocks = length (r) - 1 ;

info = [0 0 0] ;    % [nnz(L), nnz(U), nnz(F)], optional 2nd output

%-------------------------------------------------------------------------------
% solve the system
%-------------------------------------------------------------------------------

if (nblocks == 1 | sprank (A) < n)					    %#ok

    %---------------------------------------------------------------------------
    % matrix is irreducible or structurally singular
    %---------------------------------------------------------------------------

    [x info2] = umfpack2 (A, '\', b, Control) ;
    info = [info2(78) info2(79) 0] ;

else

    %---------------------------------------------------------------------------
    % A (p,q) is in block triangular form
    %---------------------------------------------------------------------------

    b = b (p,:) ;
    A = A (p,q) ;
    x = zeros (size (b)) ;

    %---------------------------------------------------------------------------
    % merge adjacent singletons into a single upper triangular block
    %---------------------------------------------------------------------------

    [r, nblocks, is_triangular] = merge_singletons (r) ;

    %---------------------------------------------------------------------------
    % solve the system: x (q) = A\b
    %---------------------------------------------------------------------------

    for k = nblocks:-1:1

	% get the kth block
        k1 = r (k) ;
        k2 = r (k+1) - 1 ;

	% solve the system
        [x2 info2] = solver (A (k1:k2, k1:k2), b (k1:k2,:), ...
	    is_triangular (k), Control) ;
	x (k1:k2,:) = x2 ;

        % off-diagonal block back substitution
        F2 = A (1:k1-1, k1:k2) ;
        b (1:k1-1,:) = b (1:k1-1,:) - F2 * x (k1:k2,:) ;

        info (1:2) = info (1:2) + info2 (1:2) ;
        info (3) = info (3) + nnz (F2) ;

    end

    x (q,:) = x ;

end

%-------------------------------------------------------------------------------
% merge_singletons
%-------------------------------------------------------------------------------

function [r, nblocks, is_triangular] = merge_singletons (r)
%
% Given r from [p,q,r] = dmperm (A), where A is square, return a modified r that
% reflects the merger of adjacent singletons into a single upper triangular
% block.  is_triangular (k) is 1 if the kth block is upper triangular.  nblocks
% is the number of new blocks.

nblocks = length (r) - 1 ;
bsize = r (2:nblocks+1) - r (1:nblocks) ;
t = [0 (bsize == 1)] ;
z = (t (1:nblocks) == 0 & t (2:nblocks+1) == 1) | t (2:nblocks+1) == 0 ;
y = [(find (z)) nblocks+1] ;
r = r (y) ;
nblocks = length (y) - 1 ;
is_triangular = y (2:nblocks+1) - y (1:nblocks) > 1 ;

%-------------------------------------------------------------------------------
% solve Ax=b, but check for small and/or triangular systems
%-------------------------------------------------------------------------------

function [x, info] = solver (A, b, is_triangular, Control)
if (is_triangular)
    % back substitution only
    x = A \ b ;
    info = [nnz(A) 0 0] ;
elseif (size (A,1) < 4)
    % a very small matrix, solve it as a dense linear system
    x = full (A) \ b ;
    n = size (A,1) ;
    info = [(n^2+n)/2 (n^2+n)/2 0] ;
else
    % solve it as a sparse linear system
    [x info] = umfpack_solve (A, '\', b, Control) ;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [] = main(train_file, test_file)
    %===============================================
    %              Reading the files
    %===============================================
    [train_data_col1, train_data_col2, train_class, train_lenght] = load_data(train_file);
    [test_data_col1, test_data_col2, test_class, test_lenght] = load_data(test_file);

    %===============================================
    %                Perfrom DTW
    %===============================================
    dtw(train_data_col1, train_data_col2, train_class, train_lenght, test_data_col1, test_data_col2, test_class, test_lenght)
end

repo: Machine-Learning-Models


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function y = psabr_complex(a,b,r,n,f,k,t)
% sabr risk neutral density using Benhamou et al. pricing formula

    eps = 1e-006;

    y1 = sabr_complex(f,k+eps,t,a,b,r,n);
    y2 = sabr_complex(f,k,t,a,b,r,n);
    y3 = sabr_complex(f,k-eps,t,a,b,r,n);
    y = (y1-2*y2+y3)/eps^2;
    y(isnan(y)) = 0;
end

author: Sable

repo: mcbench-benchmarks


text:
function out = CO_Embed2_AngleTau(y,maxTau)
% CO_Embed2_AngleTau Angle autocorrelation in a 2-dimensional embedding space
%
% Investigates how the autocorrelation of angles between successive points in
% the two-dimensional time-series embedding change as tau varies from
% tau = 1, 2, ..., maxTau.
%
%---INPUTS:
% y, a column vector time series
% maxTau, the maximum time lag to consider

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

doPlot = false;
tauRange = (1:1:maxTau);
numTau = length(tauRange);

% Ensure y is a column vector
if size(y,2) > size(y,1);
	y = y';
end

stats_store = zeros(3,numTau);

for i = 1:numTau
	tau = tauRange(i);

	m = [y(1:end-tau), y(1+tau:end)];

	theta = diff(m(:,2))./diff(m(:,1));
	theta = atan(theta); % measured as deviation from the horizontal

	if isempty(theta)
		error('Time series (N=%u) too short for embedding',length(y));
	end

	stats_store(1,i) = CO_AutoCorr(theta,1,'Fourier');
	stats_store(2,i) = CO_AutoCorr(theta,2,'Fourier');
	stats_store(3,i) = CO_AutoCorr(theta,3,'Fourier');
end

if doPlot
    figure('color','w'); box('on');
    plot(stats_store');
end

% ------------------------------------------------------------------------------
% Compute lots of outputs statistics:
% ------------------------------------------------------------------------------
out.ac1_thetaac1 = CO_AutoCorr(stats_store(1,:),1,'Fourier');
out.ac1_thetaac2 = CO_AutoCorr(stats_store(2,:),1,'Fourier');
out.ac1_thetaac3 = CO_AutoCorr(stats_store(3,:),1,'Fourier');
out.mean_thetaac1 = mean(stats_store(1,:));
out.max_thetaac1 = max(stats_store(1,:));
out.min_thetaac1 = min(stats_store(1,:));
out.mean_thetaac2 = mean(stats_store(2,:));
out.max_thetaac2 = max(stats_store(2,:));
out.min_thetaac2 = min(stats_store(2,:));
out.mean_thetaac3 = mean(stats_store(3,:));
out.max_thetaac3 = max(stats_store(3,:));
out.min_thetaac3 = min(stats_store(3,:));
out.meanrat_thetaac12 = out.mean_thetaac1/out.mean_thetaac2;
out.diff_thetaac12 = sum(abs(stats_store(2,:)-stats_store(1,:)));

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [signal, state] = flt_beamforming(varargin)
% Return the current source density for a given head model and data using
% a simple beamforming method.
% The reconstructed CSD time-series (or source potential maps) will be
% stored in signal.srcpot. This matrix has dimension [num_voxels x num_samples].
%
% Author: Tim Mullen, Jan 2013, SCCN/INC/UCSD
%         Alejandro Ojeda, Jan 2013, SCCN/INC/UCSD
%         Christian Kothe, Jan 2013, SCCN/INC/UCSD

if ~exp_beginfun('filter'), return; end

declare_properties('name','Beamforming', 'experimental',true, 'independent_channels',false, 'independent_trials',false);

arg_define(varargin, ...
    arg_norep({'signal','Signal'}), ...
    arg_nogui({'lead_field','LeadField','K','ForwardModel'},[],[],'Forward model (matrix)','shape','matrix'), ...
    arg({'cov_range','CovRange'},'whole',{'whole','epoch'},'Covariance data range. This is the range over which the covariance shall be calculated.'), ...
    arg({'cov_shrinkage','CovShrinkage'},0,[],'Covariance shrinkage parameter. For better conditioning.'), ...
    arg({'cov_robust','CovRobust'},false,[],'Use robust covariance estimate.'), ...
    arg({'rescale_sources','RescaleSources','Rescale'},false,[],'Whether to rescale the solution. If false, the source activity scale will be incorrect, but downstream methods may be scale-invariant.'), ...
    arg({'standardize_trials','StandardizeTrials'},false,[],'Perform per-trial standardization'), ...
    arg_nogui({'state','State'},[],[],'State object. When provided, hyperparameters will be estimated adaptively from prior state'));

[dummy,S,T] = size(signal.data); %#ok<ASGLU>

if strcmp(cov_range,'whole')
    % use whole-data covariance matrix
    if isempty(state)
        LF = lead_field;
        if cov_robust
            C = cov_blockgeom(signal.data(:,:)');
        else
            C = cov(signal.data(:,:)');
        end
        C = (1-cov_shrinkage)*C + cov_shrinkage*mean(trace(C))*eye(length(C));
        state.srcweights = LF'/C;
        if rescale_sources
            if ndims(LF) == 3
                error('Rescaling not yet implemented for vectorial lead-field matrices.'); end
            for k=size(LF,2):-1:1
                scales(k) = 1./((LF(:,k)'/C)*LF(:,k)); end
            state.srcweights = bsxfun(@times,scales(:),state.srcweights);
        end
    end
    signal.srcpot = reshape(state.srcweights*signal.data(:,:),[],S,T);
elseif strcmp(cov_range,'epoch')
    % use per-epoch covariance matrix
    if isempty(signal.epoch)
        error('Trying to use epoch-wise covariance matrix, but data not epoched.'); end
    for t=T:-1:1
        LF = lead_field;
        if cov_robust
            C = cov_blockgeom(signal.data(:,:,t)');
        else
            C = cov(signal.data(:,:,t)');
        end
        C = (1-cov_shrinkage)*C + cov_shrinkage*mean(trace(C))*eye(length(C));
        srcweights = LF'/C;
        if rescale_sources
            if ndims(LF) == 3
                error('Rescaling not yet implemented for vectorial lead-field matrices.'); end
            for k=size(LF,2):-1:1
                scales(k) = 1./((LF(:,k)'/C)*LF(:,k)); end
            srcweights = bsxfun(@times,scales(:),srcweights);
        end
        signal.srcpot(:,:,t) = srcweights*signal.data(:,:,t);
    end
end
if standardize_trials
    signal.srcpot = bsxfun(@minus,signal.srcpot,mean(signal.srcpot,2));
    signal.srcpot = bsxfun(@times,signal.srcpot,1./(std(signal.srcpot,[],2)+eps));
end

exp_endfun;

repo: REST


text:
% A program of type function
% Program name is def611.m
% To accompany the textbook:
% Applications of MATLAB: Numerical Solutions.
% By Yasin A. Shiboul
% This function defines the differential equation required to be solved in  question 11 of chapter 6
% the differential equation  " dy/dt = y+t".
%
function dy=def611(t,y)
dy=[y+t]


text:
function [GHz] = kHz2GHz(kHz)
% Convert frequency from kilohertz to gigahertz.
% Chad A. Greene 2012
GHz = kHz*1e-6;


text:
function [occGrad,occlusionErrors] = occlusionGrad(S,state)

occGrad = zeros(size(S));
numpoints = size(S,1);

mask = double(state.mask);
[Y,X]=size(mask);
mkp = mean(state.kps,2);% we had shifted everything by this when computing projected points so shift back now

%% Finding points projected outside the image
points2d = state.cameraScale*state.cameraRot*S'; %% Verify this later
points2d = points2d' + repmat(mkp',size(S,1),1);
points2d = round(points2d(:,1:2)); %% Verify correctness later

badX = double(points2d(:,1)>X) + double(points2d(:,1)<1);
badY = double(points2d(:,2)>Y) + double(points2d(:,2)<1);

%% Visualization
%imagesc(mask);hold on;
%plot(points2d(:,1),points2d(:,2));pause();
%close;

%% Restructuring to ease computation to determine of projected points inside silhoutte
points2d = points2d*[Y;1]-Y;
points2d(points2d>X*Y)=X*Y;
points2d(points2d<1)=1;

mask = double(mask(:)==0);
occlusionErrors = double((mask(points2d)+badX+badY)>0);
occGrad = repmat(occlusionErrors,1,3).*(repmat(mean(S,1),numpoints,1)-S);

end

repo: CategoryShapes


text:
function CreateCenteredMIPanel(handles)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

    %CenteredMI Panel
    handles.centeredMI = uipanel('Parent',handles.toolPanel,'Units','normalize','FontSize',10,'Title','MI Optimizer Setup',...
                                    'Tag','centeredMIPanel','Clipping','off',...
                                    'Position',[0 0 1 1],...
                                    'Visible','on');
    dx = 0; dy = -90;
    h337 = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[8+dx 370-5+dy 133 29],...
    'String','Histogram Bins:',...
    'Style','text',...
    'Tag','text68');

    handles.centeredMI_histbins = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[148+dx 368+dy 116 30],...
    'String','32',...
    'Style','edit',...
    'Tag','centeredMI_histbins');

    h339 = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[6+dx 332-8+dy 135 32],...
    'String','Initial Radius:',...
    'Style','text',...
    'Tag','text69');

    handles.centeredMI_initradius = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[148+dx 329+dy 116 30],...
    'String','0.05',...
    'Style','edit',...
    'Tag','centeredMI_initradius');

    h341 = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[6+dx 295-5+dy 135 30],...
    'String','Epsilon:',...
    'Style','text',...
    'Tag','text70');

    handles.centeredMI_epsilon = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[148+dx 291+dy 116 30],...
    'String','0.001',...
    'Style','edit',...
    'Tag','centeredMI_epsilon');

    h343 = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[6+dx 260-5+dy 135 28],...
    'String','Max Iteration Number:',...
    'Style','text',...
    'Tag','text71');

    handles.centeredMI_iternum = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[148+dx 256+dy 116 30],...
    'String','1500',...
    'Style','edit',...
    'Tag','centeredMI_iternum');

%     h345 = uicontrol(...
%     'Parent',handles.centeredMI,...
%     'FontSize',10,...
%     'HorizontalAlignment','right',...
%     'Position',[8+dx 219+dy 133 30],...
%     'String','Initial Angle:',...
%     'Style','text',...
%     'Tag','text72');
%
%     h346 = uicontrol(...
%     'Parent',handles.centeredMI,...
%     'FontSize',10,...
%     'HorizontalAlignment','right',...
%     'Position',[7+dx 188+dy 134 27],...
%     'String','Initial TranslationX:',...
%     'Style','text',...
%     'Tag','text73');
%
%     handles.centeredMI_initangle = uicontrol(...
%     'Parent',handles.centeredMI,...
%     'FontSize',10,...
%     'HorizontalAlignment','left',...
%     'Position',[148+dx 223+dy 79 25],...
%     'String','0',...
%     'Style','edit',...
%     'Tag','centeredMI_initangle');
%
%     handles.centeredMI_initTx = uicontrol(...
%     'Parent',handles.centeredMI,...
%     'FontSize',10,...
%     'HorizontalAlignment','left',...
%     'Position',[148+dx 189+dy 79 26],...
%     'String','0',...
%     'Style','edit',...
%     'Tag','centeredMI_initTx');
%
%     h349 = uicontrol(...
%     'Parent',handles.centeredMI,...
%     'FontSize',10,...
%     'HorizontalAlignment','right',...
%     'Position',[6+dx 156+dy 135 26],...
%     'String','Initial TranslationY:',...
%     'Style','text',...
%     'Tag','text74');
%
%     handles.centeredMI_initTy = uicontrol(...
%     'Parent',handles.centeredMI,...
%     'FontSize',10,...
%     'HorizontalAlignment','left',...
%     'Position',[148+dx 157+dy 79 25],...
%     'String','0',...
%     'Style','edit',...
%     'Tag','centeredMI_initTy');

    h351 = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[6+dx 153-3+dy 135 28],...
    'String','Rotation Scale:',...
    'Style','text',...
    'Tag','text75');

    handles.centeredMI_rs = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[148+dx 153+dy 116 30],...
    'String','0.6',...
    'Style','edit',...
    'Tag','centeredMI_rs');

    h353 = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[8+dx 180+dy 133 30],...
    'String','Center Scale X:',...
    'Style','text',...
    'Tag','text76');

    h354 = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[7+dx 213+dy 134 27],...
    'String','Center Scale Y:',...
    'Style','text',...
    'Tag','text77');

    handles.centeredMI_rcsx = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[148+dx 189+dy 79 25],...
    'String','0',...
    'Style','edit',...
    'Tag','centeredMI_rcsx');

    handles.centeredMI_rcsy = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[148+dx 219+dy 79 26],...
    'String','0',...
    'Style','edit',...
    'Tag','centeredMI_rcsy');

    h357 = uicontrol(...
    'Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','right',...
    'Position',[6+dx 121+dy 135 26],...
    'String','Translation Scale:',...
    'Style','text',...
    'Tag','text78');

    handles.centeredMI_ts = uicontrol('Parent',handles.centeredMI,...
    'FontSize',10,...
    'HorizontalAlignment','left',...
    'Position',[148+dx 123+dy 79 25],...
    'String','0.001',...
    'Style','edit',...
    'Tag','centeredMI_ts');

    guidata(handles.mainframe, handles);

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function plotPredictionResults_STS(pathResults,fSetName,metrics,maxOrder)
% -------------------------------------------------------------------------
% function plotPredictionResults_STS(pathResults,fSetName,metrics,maxOrder)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function plots prediction performance estimation results for all the
% different feature set types entered as inputs. See ref. [1] for more
% details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% - pathResults: Full path to the 'RESULTS' folder where prediction results
%                are saved.
% - fSetName: Cell of strings specifying the name of the type of feature
%                 sets to analyze.
%                 Example: fSetName = {'PET','SEPARATE','FUSED'}
% - metrics: Name of the metrics to show in the plots.
%            Example: metrics = {'AUC632','Sensitivity632','Specificity632'}
% - maxOrder: Integer specifying the maximal multivariable model order.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: May 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
cd(pathResults)

signs = {'-r',':b','--g'};
nFSET = numel(fSetName);
nMetrics = numel(metrics);
figure
for i = 1:nFSET
    val = zeros(maxOrder,numel(metrics));
    val_SE = zeros(maxOrder,nMetrics);
    results = load(['RESULTS_',fSetName{i},'_BEST']); results = struct2cell(results); results = results{1};
    for j = 1:maxOrder
        orderName = ['Order',num2str(j)];
        for k =1:nMetrics
            val(j,k) = results.(orderName).(metrics{k});
            val_SE(j,k) = results.(orderName).(['SE_',metrics{k}]);
        end
    end
    subplot(1,nFSET,i)
    for k = 1:nMetrics
        errorbar(1:maxOrder,val(:,k),val_SE(:,k),signs{k},'LineWidth',2,'MarkerFaceColor',signs{k}(end),'MarkerSize',4)
        hold on
    end
    xlabel('Model Order','FontSize',24)
    ylabel('Prediction performance','FontSize',24)
    title(fSetName{i},'FontSize',30,'FontWeight','bold')
    legend(metrics,'Location','SouthEast')
    axis([0 maxOrder+1 0.5 1])
    set(gca,'FontSize',20)
    set(gca,'XTick',[1 2 3 4 5 6 7 8 9 10])
    hold off
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function [ offspring ] = DE_transfer(Problem, Population1, Population2, popsize)

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Kangjia Qiao

    Fm  = [0.6,0.8,1.0];
    CRm = [0.1,0.2,1.0];

    index = randi([1,length(Fm)],popsize,1);
    F     = Fm(index);
    F     = F';
    index = randi([1,length(CRm)],popsize,1);
    CR    = CRm(index);
    CR    = CR';

    index =randi(Problem.N, popsize,1);

    permutation = randperm(Problem.N);

    array = permutation(1:popsize);

    pop1 = Population1(array).decs;
    pop2 = Population2.decs;

    vi = pop2(index,:);

    mask  = rand(popsize, Problem.D) > CR(:, ones(1, Problem.D)); % mask is used to indicate which elements of ui comes from the parent
    rows  = (1 : popsize)'; cols = floor(rand(popsize, 1) * Problem.D)+1; % choose one position where the element of ui doesn't come from the parent
    jrand = sub2ind([popsize Problem.D], rows, cols); mask(jrand) = false;
    u     = vi; u(mask) = pop1(mask);

    offspring = Problem.Evaluation(u);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
close all;
clear;
clc
addpath(genpath('./'));

% Path planner
disp('Planning ...');
% This map is used to generate path and check collision
map_colli = load_map_inflated('maps/map1_full.txt', 0.5, 0.2, 0); % xy res, z res, margin
% This map is for visualization (without boundary obsatacles)
map = load_map('maps/map1.txt', 0.5, 0.15, 0.34); % 3 parameters are dummy

% start = {[5.0 -1 3.5]};
% stop  = {[5.0 19.0 .5]};

start = [5.0 -1 3.5];
stop  = [5.0 19.0 .5];
% stop  = [17.0 4.0 .5];

% Generate way points
path = dijkstra(map_colli, start, stop, true);

figure(1);
plot_path(map, path);
axis equal
grid on
title('Quad Simulator');xlabel('x');ylabel('y');zlabel('z');
hold on
plot3(start(1),start(2),start(3), '*r', 'markersize', 10)
plot3(stop(1),stop(2),stop(3), '*g', 'markersize', 10)
hold off
view(3)
%% Desired trajectory generation
% *Warning*
% Since demonstrated obstacles are dense, including narrow corridors, it
% may take up to a minute to get collision free trajectory.
% It is highly recommended to use simpler map if you just need a
% demonstration.
traj_obj = trajectoryGenerator(map_colli, path);

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ODE Simulation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tau_vec = traj_obj.tau_vec;
path = traj_obj.path;
ts = [0 cumsum(tau_vec)];

%%%%%%%%%%%%%%%%%%%%%%%%%%% Initial Condition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% State: [x, v, pitch, roll, yaw, w]'
x0 = [path(1,:) zeros(1,9)]';

% Disturbances on the initial condition
% x0(1:3) = x0(1:3)-[.1 .05 .15]';
% x0(4:6) = x0(4:6)+[.01 .01 .02]';

% Initial error of yaw angle must not exceed 90 degrees
% x0(7:9) = x0(7:9)+ones(3,1)*pi/18;
% x0(10:12) = x0(10:12) + .0;
%%%%%%%%%%%%%%%%%%%%%%%%%%% Initial Condition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Model parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model_param.grav = 9.807;
model_param.mass = 1.477;
model_param.I = [0.01152 0 0;0 0.01152 0;0 0 0.0218];
model_param.arm_length = 0.263;
model_param.c_tf = 8.004e-4;

% Gain
KK.Kp = 5;
KK.Kv = 5;
KK.KR = 5;
KK.K_omega = 5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Model parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% For real time video
target_fps = 30;
t_sim = 0:1/target_fps:sum(tau_vec);

% External force
start_T = 0;
duration = 0;
mag = 0;
direction = [1 0 0]';
Fmat = extForce_gen(t_sim, start_T, duration, mag, direction);

% Run ode45
options =  odeset('RelTol',1e-2,'AbsTol',1e-2);
[tsave, xsave] = quadSim(traj_obj, model_param, KK, t_sim, x0, Fmat, options);

%% Desired trajectory vs actual trajectory
desired_pos = zeros(length(tsave),3);
desired_vel = zeros(length(tsave),3);
desired_acc = zeros(length(tsave),3);
desired_jerk = zeros(length(tsave),3);
desired_yaw = zeros(length(tsave),1);
for i = 1:length(tsave)
    desired_s = desiredState(traj_obj, tsave(i));
    desired_pos(i,:) = desired_s.pos';
    desired_vel(i,:) = desired_s.vel';
    desired_acc(i,:) = desired_s.acc';
    desired_jerk(i,:) = desired_s.jerk';
    desired_yaw(i) = desired_s.yaw;
end

figure(1)
% Desired trajectory
plot3(desired_pos(:,1),desired_pos(:,2),desired_pos(:,3))
hold on
grid on
axis equal
% Actual trajectory
plot3(xsave(:,1),xsave(:,2),xsave(:,3),'--')
% legend('Desired','Actual')
%%
figure(5)
subplot(4,1,1)
plot(tsave,xsave(:,7)*180/pi,'-b','LineWidth',1.0);title('Euler angles');ylabel('roll, \phi')
grid on
subplot(4,1,2)
plot(tsave,xsave(:,8)*180/pi,'-b','LineWidth',1.0);ylabel('pitch, \theta')
grid on
subplot(4,1,3)
plot(tsave,xsave(:,9)*180/pi,'-b','LineWidth',1.0);ylabel('yaw, \psi');%xlabel('time, sec')
grid on

%%
m = model_param.mass;
g = model_param.grav;
% pos error
ep = xsave(:,1:3)-desired_pos;
% vel error
ev = xsave(:,4:6)-desired_vel;
% accel error
calc_acc = diff(xsave(:,4:6))./diff(tsave);
% plot(tsave(1:end-1),calc_acc)
%
ea = calc_acc-desired_acc(1:end-1,:);
%%
% desired force
Fd = -KK.Kp.*ep -KK.Kv.*ev + m*desired_acc;
Fd(:,3) = Fd(:,3) + m*g;
% Desired force derivative
% Fd_dot = -KK.Kp*ev -KK.Kv*ea + m*desired_jerk;
% desired rotation, Rd = [xbd ybd zbd]
zbd = zeros(length(tsave),3);
xcd = zeros(length(tsave),3);
zbdxcd = zeros(length(tsave),3);
norm_zbdxcd = zeros(length(tsave),1);
for i=1:length(tsave)
    zbd(i,:) = Fd(i,:)/norm(Fd(i,:));
    xcd(i,:) = [cos(desired_yaw(i)) sin(desired_yaw(i)) 0];
    dummy = hat_optr(zbd(i,:))*xcd(i,:)';
    zbdxcd(i,:) = dummy';
    norm_zbdxcd(i) = norm(zbdxcd(i,:));
end
figure(5)
subplot(4,1,4)
plot(tsave,norm_zbdxcd,'-','LineWidth',1.5)
grid on
%% Video generator
figure(1);
% plot_path(map, path{1});
filename = 'my_video.avi';
video_gen(tsave, xsave, filename, 15, Fmat)

author: yorgoon

repo: minimum-snap-geometric-control


text:
function t_opf_dc_ot(quiet)
%T_OPF_DC_OT  Tests for DC optimal power flow using Opt Tbx solvers.

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

if nargin < 1
    quiet = 0;
end

matlab = have_feature('matlab');
if matlab
    if have_feature('quadprog_ls')
        if have_feature('optimoptions')
            if have_feature('linprog_ds')
                if have_feature('quadprog', 'vnum') >= 7.005  %% R2016b and later
                    algs  = {'interior-point', 'dual-simplex'};
                else
                    algs  = {'interior-point', 'active-set', 'simplex', 'dual-simplex'};
                end
            else
                algs  = {'interior-point', 'simplex'};
            end
        else
            algs  = {'interior-point', 'active-set'};
        end
    else
        algs  = {''};
    end
else
    algs  = {'default'};
end

num_tests = 43 * length(algs);

t_begin(num_tests, quiet);

[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

casefile = 't_case9_opf';
if quiet
    verbose = 0;
else
    verbose = 0;
end
if have_feature('octave')
    if have_feature('octave', 'vnum') >= 4
        file_in_path_warn_id = 'Octave:data-file-in-path';
    else
        file_in_path_warn_id = 'Octave:load-file-in-path';
    end
    s1 = warning('query', file_in_path_warn_id);
    warning('off', file_in_path_warn_id);
end

mpopt = mpoption('out.all', 0, 'verbose', verbose);
mpopt = mpoption(mpopt, 'opf.dc.solver', 'OT');

warning off optim:linprog:IgnoreStartPoint;
warning off optim:linprog:AlgOptsWillError;

%% run DC OPF
if have_feature('quadprog')
  for k = 1:length(algs)
    if matlab
        if ~isempty(algs)
            mpopt = mpoption(mpopt, 'linprog.Algorithm', algs{k});
        else
            mpopt = mpoption(mpopt, 'linprog', []);
        end
    end
    if strcmp(algs{k}, 'dual-simplex') && strcmp(have_feature('fmincon', 'vstr'), '7.1')
        have_prices = 0;    %% dual-simplex did not return prices in MATLAB R2014b!?!
    else
        have_prices = 1;
    end
    t0 = sprintf('DC OPF (OT %s): ', algs{k});

    %% set up indices
    ib_data     = [1:BUS_AREA BASE_KV:VMIN];
    ib_voltage  = [VM VA];
    ib_lam      = [LAM_P LAM_Q];
    ib_mu       = [MU_VMAX MU_VMIN];
    ig_data     = [GEN_BUS QMAX QMIN MBASE:APF];
    ig_disp     = [PG QG VG];
    ig_mu       = (MU_PMAX:MU_QMIN);
    ibr_data    = (1:ANGMAX);
    ibr_flow    = (PF:QT);
    ibr_mu      = [MU_SF MU_ST];
    ibr_angmu   = [MU_ANGMIN MU_ANGMAX];

    %% get solved DC power flow case from MAT-file
    load soln9_dcopf;       %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF
    t = t0;
    [baseMVA, bus, gen, gencost, branch, f, success, et] = rundcopf(casefile, mpopt);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    if have_prices
        t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    else
        t_skip(1, [t 'bus lam: lambdas not computed']);
    end
    if have_prices
        t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    else
        t_skip(1, [t 'bus mu: lambdas not computed']);
    end
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    if have_prices
        t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    else
        t_skip(1, [t 'gen mu: lambdas not computed']);
    end
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    if have_prices
        t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
    else
        t_skip(1, [t 'branch mu: lambdas not computed']);
    end

    %%-----  test OPF with angle difference limits  -----
    t = [t0 'w/angle diff lims : '];
    mpc = loadcase(casefile);
    mpc.branch(4, ANGMAX) = 3;
    mpc.branch(7, ANGMIN) = -4.5;
    r = rundcopf(mpc, mpopt);
    [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success);
    t_ok(success, [t 'success']);
    t_is(   f, 6456.7213, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,PG        ),    [99.98497;89.35133;125.66371], 4, [t 'gen dispatch']);
    t_is(branch(:,ibr_data  ), mpc.branch(:,ibr_data   ), 10, [t 'branch data']);
    e = zeros(size(branch, 1), 1);
    e(4) = 297.83776;
    e(7) = -26.94788;
    if have_prices
        t_is(branch(:,MU_ANGMAX )-branch(:,MU_ANGMIN ), e, 4, [t 'branch ang diff mu']);
    else
        t_skip(1, [t 'branch ang diff mu: lambdas not computed']);
    end

    t = [t0 'w/ignored angle diff lims : '];
    mpopt1 = mpoption(mpopt, 'opf.ignore_angle_lim', 1);
    r = rundcopf(mpc, mpopt1);
    [bus, gen, branch, f, success] = deal(r.bus, r.gen, r.branch, r.f, r.success);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    if have_prices
        t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    else
        t_skip(1, [t 'bus lam: lambdas not computed']);
    end
    if have_prices
        t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    else
        t_skip(1, [t 'bus mu: lambdas not computed']);
    end
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    if have_prices
        t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    else
        t_skip(1, [t 'gen mu: lambdas not computed']);
    end
    t_is(branch(:,ibr_data  ), mpc.branch(:,ibr_data   ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    if have_prices
        t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
    else
        t_skip(1, [t 'branch mu: lambdas not computed']);
    end

    %%-----  run OPF with extra linear user constraints & costs  -----
    %% two new z variables
    %%      0 <= z1, P3 - P1 <= z1
    %%      0 <= z2, P3 - P2 <= z2
    %% with A and N sized for DC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1;1;2;2;2],[10;12;13;12;11;14],[-1;1;-1;1;-1;-1],2,14);
    mpc.u = [0; 0];
    mpc.l = [-Inf; -Inf];
    mpc.zl = [0; 0];

    mpc.N = sparse([1;2], [13;14], [1;1], 2, 14);   %% new z variables only
    mpc.fparm = ones(2,1) * [1 0 0 1];              %% w = r = z
    mpc.H = sparse(2,2);                            %% no quadratic term
    mpc.Cw = [1000;1];

    t = [t0 'w/extra constraints & costs 1 : '];
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(r.gen(1, PG), 116.15974, 5, [t 'Pg1 = 116.15974']);
    t_is(r.gen(3, PG), 116.15974, 5, [t 'Pg3 = 116.15974']);
    t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
    t_is(r.cost.usr, 0.3348, 4, [t 'user costs']);

    %% with A and N sized for AC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1;1;2;2;2],[19;21;25;21;20;26],[-1;1;-1;1;-1;-1],2,26);
    mpc.u = [0; 0];
    mpc.l = [-Inf; -Inf];
    mpc.zl = [0; 0];

    mpc.N = sparse([1;2], [25;26], [1;1], 2, 26);   %% new z variables only
    mpc.fparm = ones(2,1) * [1 0 0 1];              %% w = r = z
    mpc.H = sparse(2,2);                            %% no quadratic term
    mpc.Cw = [1000;1];

    t = [t0 'w/extra constraints & costs 2 : '];
    [r, success] = rundcopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(r.gen(1, PG), 116.15974, 5, [t 'Pg1 = 116.15974']);
    t_is(r.gen(3, PG), 116.15974, 5, [t 'Pg3 = 116.15974']);
    t_is(r.var.val.z, [0; 0.3348], 4, [t 'user vars']);
    t_is(r.cost.usr, 0.3348, 4, [t 'user costs']);

    t = [t0 'infeasible : '];
    %% with A and N sized for DC opf
    mpc = loadcase(casefile);
    mpc.A = sparse([1;1], [10;11], [1;1], 1, 14);   %% Pg1 + Pg2
    mpc.u = Inf;
    mpc.l = 600;
    %% Opt Tbx 5.1 has fatal bug here!
    try
        [r, success] = rundcopf(mpc, mpopt);
    catch
        success = 0;
    end
    t_ok(~success, [t 'no success']);

    %% OPF with all buses isolated
    t = [t0 'all buses isolated : '];
    mpc = loadcase(casefile);
    mpc.bus(:, BUS_TYPE) = NONE;
    try
        r = rundcopf(mpc, mpopt);
        t_is(r.success, 0, 12, [t 'success = 0']);
    catch
        t_ok(0, [t 'unexpected fatal error']);
    end
  end
else
    t_skip(num_tests, 'Optimization Toolbox not available');
end

if have_feature('octave')
    warning(s1.state, file_in_path_warn_id);
end

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [g1,g2] = sdlfmaXsdrbfKernGradient(sdlfmaKern, sdrbfKern, t1, ...
    t2, covGrad)

% SDLFMAXSDRBFKERNGRADIENT Gradients cross kernel between a SDLFM and SDRBF
% FORMAT
% DESC computes a cross gradient for a cross kernel between a switching
% dynamical LFM kernel (accel.) and a switching dynamical RBF kernel.
% ARG sdlfmvKern : the kernel structure associated with the SDLFM
% kernel (acceleration).
% ARG sdrbfKern : the kernel structure associated with the SDRBF kernel.
% ARG t1 : inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the first kernel, for
% ordering see lfmKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for
% ordering see lfmKernExtractParam.
%
% FORMAT
% DESC computes a cross gradient for a cross kernel between a switching
% dynamical LFM kernel (accel.) and a switching dynamical RBF kernel.
% ARG sdlfmvKern : the kernel structure associated with the SDLFM
% kernel (acceleration).
% ARG sdrbfKern : the kernel structure associated with the SDRBF kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% ARG covGrad : gradient of the objective function with respect to
% the elements of the cross kernel matrix.
% RETURN g1 : gradient of the parameters of the first kernel, for
% ordering see lfmKernExtractParam.
% RETURN g2 : gradient of the parameters of the second kernel, for
% ordering see lfmKernExtractParam.
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

if nargin == 4
    covGrad = t2;
    t2 = t1;
end

[g1, g2] = sdlfmXsdrbfKernGradient(sdlfmaKern, sdrbfKern, t1, t2, covGrad, 'Accel');

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function K = calWholeStiffnessMatrix(coord, unit_topology_table, materials,cal_type)
% 这个函数计算一维半带宽方法下的整体的劲度矩阵K
% 输入：
%     节点坐标表
%     单元拓扑表
%     材料情况
%       计算模式 cal_type
%            单元应力问题
%            单元应变问题
% 输出：
%     整体刚度矩阵 K
%     K 为cell
%     K{1，1}为转化为一维带状储存的K
%     K{1，2}为主对角元在一维矩阵中位置

%使用openKspeace（）函数计算带状矩阵K所需要的空间
    [K, K_info] = openKspeace(unit_topology_table);
    m           = size(unit_topology_table,1);
    t           = 1;
    for i = 1:m
        element_X = coord(unit_topology_table(i, :)',1);
        element_Y = coord(unit_topology_table(i, :)',2);
        matrixB   = calMatrixB(element_X, element_Y);
        matrixD   = calMatrixD(materials(i,1), materials(i,2), cal_type);
        element_k = calElementStiffnessMatrix(matrixB, matrixD, t,...
                                              calAera(element_X, element_Y));
        position  = zeros(1,6);
        for ii = 1:3
            position(ii*2-1) = unit_topology_table(i,ii)*2 - 1;
            position(ii*2)   = unit_topology_table(i,ii)*2;
        end
        for j = 1:6
            for jj = 1:6
                if position(j) >= position(jj)
                oned_x    = K_info(position(j)) - (position(j)-position(jj));
                K(oned_x) = K(oned_x) + element_k(j,jj);
                end
            end
        end
    end
    temp = K;
    clear K;
    K{1} = temp;
    K{2} = K_info;

author: Meelfy

repo: FEM

save_path: github-repos/MATLAB/Meelfy-FEM


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

%   Copyright (C) 2012 Jonathan Currie (I2C2)

% Print out possible values of properties.
if ((nargin == 0) && (nargout == 0))
    printfields();
    return
end

%Names and Defaults
Names = {'swarm_size';'vectorized';'mu';'nu';'iweight';'fweight';'delta';'idelta';'ddelta'};
Defaults = {42;0;0.5;0.5;0.9;0.4;Inf;2.0;0.5};

%Enter and check user args
try
    options = opticheckset(Names,Defaults,@checkfield,varargin{:});
catch ME
    throw(ME);
end

function checkfield(field,value)
%Check a field contains correct data type
switch lower(field)
    %Scalar double
    case {'mu','nu','iweight','fweight','delta','idelta','ddelta'}
        err = opticheckval.checkScalarDbl(value,field);
    %Scalar 0/1
    case 'vectorized'
        err = opticheckval.checkScalar01(value,field);
    %Non-zero scalar double
    case 'swarm_size'
        err = opticheckval.checkScalarIntGrtZ(value,field);

    otherwise
        err = sprintf('Unrecognized parameter name ''%s''.', field);
end
if(~isempty(err)), throw(err); end

function printfields()
%Print out fields with defaults
fprintf('        swarm_size: [ Swarm Size {42} ] \n');
fprintf('        vectorized: [ Objective function is vectorized {0}, 1 ] \n');
fprintf('                mu: [ Cognitial Parameter {0.5} ] \n');
fprintf('                nu: [ Social Parameter {0.5} ] \n');
fprintf('           iweight: [ Initial Weight {0.9} ] \n');
fprintf('           fweight: [ Final Weight {0.4} ] \n');
fprintf('             delta: [ Initial Delta {Inf} ] \n');
fprintf('            idelta: [ Increase Delta {2.0} ] \n');
fprintf('            ddelta: [ Decrease Delta {0.5} ] \n');
fprintf('\n');


text:
% function Css = synsq_adm(type, opt)
%
% Calculate the Synchrosqueezing admissibility constant, the term
% R_\psi in Eq. 3 of [1].  Note, here we multiply R_\psi by the
% inverse of log(2)/nv (found in Alg. 1 of that paper).
%
% 1. E. Brevdo, N.S. Fučkar, G. Thakur, and H-T. Wu, "The
% Synchrosqueezing algorithm: a robust analysis tool for signals
% with time-varying spectrum," 2011.
%
% Uses numerical integration.
%
% Input:
%   type: type of wavelet (see help wfiltfn)
%   opt: options structure (wavelet parameters, see help wfiltfn)
%
% Output:
%   Css: proportional to 2*int(conj(f(w))/w, w=0..inf)
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------
function Css = synsq_adm(type, opt)
    if nargin<2, opt=struct(); end
    switch type
      % case 'sombrero',
      %   if ~isfield(opt,'s'), s = 1; else s = opt.s; end
      %   Cpsi = (4/3)*s*sqrt(pi);
      % case 'shannon',
      %   Cpsi = log(2);
      otherwise
        psihfn = wfiltfn(type, opt);
        Css = quadgk(@(x) conj(psihfn(x))./x, 0, Inf);
    end

    % Normalization constant, due to logarithmic scaling in wavelet
    % transform
    Css = Css / (sqrt(2*pi)*2*log(2));
end

author: HeLiangHIT


text:
function ord = BF_LinkageOrdering(distMat,links)
% BF_LinkageOrdering attempts to use optimalleaforder for dendrogram orderings

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

%-------------------------------------------------------------------------------
% Find the size threshold:
if any(size(distMat)==1)
    % A vector of pairwise distances
    numItems = (1+sqrt(1+8*length(distMat)))/2;
else
    % A square matrix of pairwise distances
    numItems = length(distMat);
end

%-------------------------------------------------------------------------------
% Do either optimalleaforder, or normal dendrogram ordering, applying a threhold
% on the number of items for which it is efficient to apply the optimalleaforder
% algorithm. 2000 is the number

if numItems < 2000 % small enough to try optimalleaforder
    try
        ord = optimalleaforder(links,distMat);
        % [~,~,ord] = dendrogram(links,0,'r',ord);
        fprintf(1,'Using optimalleaforder reordering!\n');
    catch
        fprintf(1,'Using dendrogram reordering.\n');
        [~,~,ord] = dendrogram(links,0);
    end
else
    fprintf(1,'Too many objects to reorder using optimalleaforder, using dendrogram instead.\n');
    [~,~,ord] = dendrogram(links,0);
end

end


text:
function test_bug1760

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_multiplotER ft_multiplotTFR

cd(dccnpath('/home/common/matlab/fieldtrip/data/test'))
load bug1760.mat

figure
ft_multiplotER(cfg, freq);


text:
% plot_UWB_channel.m
clear, clf
no_output_files = 0;  % non-zero: avoids writing output files of continuous-time responses
Ts = 0.167;        % sampling time (nsec)
num_ch=100; % number of channel impulse responses to generate
randn('state',12);  % initialize state of function for repeatability
rand('state',12);   % initialize state of function for repeatability
cm = 1;  % channel model number from 1 to 4
% get channel model params based on this channel model number
[Lam,lam,Gam,gam,nlos,sdi,sdc,sdr] = UWB_parameters(cm);
fprintf(1,['Model Parameters\n' ' Lam= %.4f, lam= %.4f, Gam= %.4f, gam= %.4f\n NLOS flag= %d, std_shdw= %.4f, std_ln_1= %.4f, std_ln_2= %.4f\n'],...
 Lam,lam,Gam,gam,nlos,sdi,sdc,sdr);
% get a bunch of realizations (impulse responses)
[h_ct,t_ct,t0,np] = UWB_model_ct(Lam,lam,Gam,gam,num_ch,nlos,sdi,sdc,sdr);
% now reduce continuous-time result to a discrete-time result
[hN,N] = convert_UWB_ct(h_ct,t_ct,np,num_ch,Ts);
% if we wanted complex baseband model or to impose some filtering function,
% this would be a good place to do it
h = resample(hN,1,N);  % decimate the columns of hN by factor N
h = h*N;  % correct for 1/N scaling imposed by decimation
channel_energy = sum(abs(h).^2);  % channel energy
h_len = size(h,1);
t = [0:(h_len-1)]*Ts;  % for use in computing excess & RMS delays

for k=1:num_ch
  % determine excess delay and RMS delay
  sq_h = abs(h(:,k)).^2/channel_energy(k);
  t_norm = t - t0(k);  % remove the randomized arrival time of first cluster
  excess_delay(k) = t_norm*sq_h;
  rms_delay(k) = sqrt((t_norm-excess_delay(k)).^2*sq_h);
  % determine # of significant paths (paths within 10 dB from peak)
  threshold_dB = -10;   % dB
  temp_h = abs(h(:,k));
  temp_thresh = 10^(threshold_dB/20)*max(temp_h);
  num_sig_paths(k) = sum(temp_h>temp_thresh);
  % determine number of sig. paths (captures x % of energy in channel)
  x = 0.85;
  temp_sort = sort(temp_h.^2);  % sorted in ascending order of energy
  cum_energy = cumsum(temp_sort(end:-1:1));  % cumulative energy
  index_e = min(find(cum_energy >= x*cum_energy(end)));
  num_sig_e_paths(k) = index_e;
end
energy_mean = mean(10*log10(channel_energy));
energy_stddev = std(10*log10(channel_energy));
mean_excess_delay = mean(excess_delay);
mean_rms_delay = mean(rms_delay);
mean_sig_paths = mean(num_sig_paths);
mean_sig_e_paths = mean(num_sig_e_paths);

fprintf(1,'Model Characteristics\n');
fprintf(1,'  Mean delays: excess (tau_m) = %.1f ns, RMS (tau_rms) = %1.f\n', ...
    mean_excess_delay, mean_rms_delay);
fprintf(1,'  # paths: NP_10dB =  %.1f, NP_85%% = %.1f\n', ...
    mean_sig_paths, mean_sig_e_paths);
fprintf(1,'  Channel energy: mean = %.1f dB, std deviation = %.1f dB\n', ...
  energy_mean, energy_stddev);

subplot(421), plot(t,h), grid on
title('Impulse response realizations'), xlabel('Time (nS)')

subplot(422), plot([1:num_ch], excess_delay, 'b-', ...
  [1 num_ch], mean_excess_delay*[1 1], 'r--' );
grid on, title('Excess delay (nS)'), xlabel('Channel number')

subplot(423), plot([1:num_ch], rms_delay, 'b-', ...
  [1 num_ch], mean_rms_delay*[1 1], 'r--' );
grid on, title('RMS delay (nS)'), xlabel('Channel number')

subplot(424), plot([1:num_ch], num_sig_paths, 'b-', ...
  [1 num_ch], mean_sig_paths*[1 1], 'r--');
grid on, title('Number of significant paths within 10 dB of peak')
xlabel('Channel number')

subplot(427), plot([1:num_ch], num_sig_e_paths, 'b-', ...
  [1 num_ch], mean_sig_e_paths*[1 1], 'r--');
grid on, title('Number of significant paths capturing > 85% energy')
xlabel('Channel number')

temp_average_power = sum(h'.*(h)')/num_ch;
temp_average_power = temp_average_power/max(temp_average_power);
average_decay_profile_dB = 10*log10(temp_average_power);
subplot(425), plot(t,average_decay_profile_dB); grid on
axis([0 t(end) -60 0]), title('Average Power Decay Profile')
xlabel('Delay (nsec)'), ylabel('Average power (dB)')

subplot(426)
figh = plot([1:num_ch],10*log10(channel_energy),'b-', ...
  [1 num_ch], energy_mean*[1 1], 'g--', ...
  [1 num_ch], energy_mean+energy_stddev*[1 1], 'r:', ...
  [1 num_ch], energy_mean-energy_stddev*[1 1], 'r:');
xlabel('Channel number'), ylabel('dB'), title('Channel Energy');
legend(figh, 'Per-channel energy', 'Mean', '\pm Std. deviation', 0)

if no_output_files, return;  end
%%% save continuous-time (time,value) pairs to files
save_fn = sprintf('cm%d_imr', cm);
% A complete self-contained file for Matlab users
save([save_fn '.mat'], 't_ct', 'h_ct', 't0', 'np', 'num_ch', 'cm');
% Two comma-delimited text files for non-Matlab users:
% File #1: cmX_imr_np.csv lisTs the number of paths in each realization
dlmwrite([save_fn '_np.csv'], np, ',');  % number of paths
% File #2: cmX_imr.csv can open with Excel
%     n'th pair of columns contains the (time,value) pairs for the n'th realization
th_ct = zeros(size(t_ct,1),2*size(t_ct,2));
th_ct(:,1:2:end) = t_ct;  % odd columns are time
th_ct(:,2:2:end) = h_ct; % even columns are values
fid = fopen([save_fn '.csv'], 'w');
if fid < 0,
  error('unable to write .csv file for impulse response, file may be open in another application');
end
for k = 1:size(th_ct,1)
  fprintf(fid,'%.4f,%.6f,', th_ct(k,1:end));
  fprintf(fid,'\r\n'); % \r\n for Windoze end-of-line
end
fclose(fid);

author: LyricYang


text:
function bc = specific_bc(xbd,ybd)
%quad_bc   channel flow boundary condition for 4:1 quadrilateral domain
%   bc = specific_bc(xbd,ybd);
%   input
%          xbd          x boundary coordinate vector
%          ybd          y boundary coordinate vector
%
%   specifies streamfunction associated with Poiseuille flow
%   needs editing for general aspect ratios
%   IFISS function: DJS; 6 March 2005.
% Copyright (c) 2005 D.J. Silvester, H.C. Elman, A. Ramage
aspect=4; deform=4;
bc=0*xbd;
k=find(xbd==0); bc(k)=2*ybd(k).*ybd(k) -(4/3)*ybd(k).*ybd(k).*ybd(k);
slope=(1/aspect)*(1/deform-1);
k3=find(ybd==slope*xbd +1);  bc(k3)=2/3;
return

author: johannesgerer

repo: jburkardt-m


text:
% std_saveindvar - save independent variables of a STUDY
%
% Usage:
%   std_saveindvar(STUDY, indvars, filename);
%   std_saveindvar(STUDY, indvars, filename, session, group);
%
% Input:
%   STUDY - EEGLAB STUDY structure
%   indvars - [cell] cell array of independent variables
%   filename - [string] file name to save variables in
%   session  - [array of value] list of session, first save all the
%              variable for the first session in the list etc...
%   group    - [array of string] list of group, first save all the
%              variable for the first group in the list etc...
%
% Authors: A. Delorme, CERCO/CNRS and SCCN/UCSD, 2022-

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

function std_saveindvar(STUDY, indvars, filename, sessions, groups, interaction_sess, interaction_group)

if nargin < 1
    help std_saveindvar;
    return;
end

if nargin < 4
    sessions = [];
    interaction_sess = false;
end
if nargin < 5
    groups = {};
    interaction_group = false;
end

if length(sessions) > 1 && length(groups) > 1
    error('Cannot process groups and session - code untested');
end

fid = fopen(filename, 'w');
if fid == -1
    error('Cannot open file for writing');
end
for iGroup = 1:max(1, length(groups))
    for iSess = 1:max(1, length(sessions))
        for iSubj = 1:length(STUDY.subject)
            for iVar  = 1:length(indvars)
                indS = strmatch(STUDY.subject{iSubj}, {STUDY.datasetinfo.subject}, 'exact');

                % sessions
                if ~isempty(sessions)
                    indSess = find([STUDY.datasetinfo.session] == sessions(iSess));
                    indS = intersect(indS, indSess);
                end

                % groups
                if ~isempty(groups)
                    indGp = strmatch(groups{iGroup}, {STUDY.datasetinfo.group}, 'exact');
                    indS = intersect(indS, indGp);
                end

                if ~isempty(indS)
                    allVals = [ STUDY.datasetinfo(indS).(indvars{iVar})];
                    if isempty(allVals), allVals = NaN; end
                    if length(unique(allVals)) > 1
                        fprintf('Warning: more than 1 value for "%s" for subject %s, taking the average\n', indvars{iVar}, STUDY.subject{iSubj});
                        allVals = mean(allVals);
                    end
                else
                    allVals = NaN;
                    fprintf('Warning: no value for "%s" subject %s (session:%d,group:%d), using NaN\n', indvars{iVar}, STUDY.subject{iSubj}, iSess, iGroup);
                end

                % interactions
                interVal11 = compute_interaction(allVals, sessions, iSess , interaction_sess);
                interVal12 = compute_interaction(allVals, groups  , iGroup, interaction_group);
                vals = [  allVals(1) interVal11 interVal12 ];
                fprintf(fid, '%s', sprintf('%f\t',vals) );
            end
            fprintf(fid, '\n');
        end
    end
end
fclose(fid);

function interVal = compute_interaction(allVals, var, iVar, interaction_var)

interVal = [];
if interaction_var
    if length(var) > 2
        error('Cannot compute interaction with more than 2 sessions');
    end
    if iVar == 1
        interVal = -allVals;
    else
        interVal = allVals;
    end
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
classdef Step3 < adi.StepBase

    %   Copyright 2017 The MathWorks, Inc.

    methods
        function obj = Step3(varargin)
            obj@adi.StepBase(varargin{:});

            %       obj.Workflow.InstallDriverWorkflow = true;

            obj.Title.Text = 'ADI Board Support Package Installer: Step3';

            obj.InfoText1.Position = [20 340 400 30];
            obj.InfoText1.Text = 'Test1: Step3';

            obj.InfoText2.Position = [20 290 400 50];
            obj.InfoText2.Text = 'Test2: Step3';

            obj.HelpText.AboutSelection = 'About Selection: Step3';
            obj.HelpText.WhatToConsider = 'Test3: Step3';
        end

%         function out = getNextScreenID(obj)
%             % HSA infrastructure requires a "getNextScreenID" method implemented
%             % in the leaf class. Otherwise, it renders "Finish" button instead of
%             % "Next".
%             out = getNextScreenIDImpl(obj);
%         end

        function restoreScreen(obj)
            obj.enableScreen();
        end
    end

    methods (Access = protected)

        function out = getNextScreenIDImpl(~)
            out = 'adi.Step3';
        end

    end

end

author: analogdevicesinc


text:
function [scan3M,dose3mC,strMaskC,xyzGridC,strColorC] = ...
    getScanDoseStrVolumes(scanNum,doseNumV,structNamC,planC)
% function [scan3M,dose3mC,strMaskC] = ...
% getScanDoseStr(scanNum,doseNumV,structnamC,planC)
%
% This function returns volumetric matrices for scan, dose and structure
% masks on the original scan grid.
%
% Example usage
% cerrFileName = 'L:\Data\RTOG0617\CERR_files_tcia_registered_0617-489880_09-09-2000-50891\0617-489880_09-09-2000-50891.mat';
% scanNum = 1;
% doseNum = 1;
% structNamC = {'DL_HEART_MT','DL_AORTA','DL_LA','DL_LV','DL_RA',...
%       'DL_RV','DL_IVC','DL_SVC','DL_PA'};
%
% % Load planC from file
% planC = loadPlanC(cerrFileName,tempdir);
% planC = updatePlanFields(planC);
% planC = quality_assure_planC(cerrFileName,planC);
% % Get scan,dose, struct volumes
% [scan3M,dose3M,strMaskC] = getScanDoseStr(scanNum,doseNum,structnamC,planC);
%
% APA, 11/2/2020

indexS = planC{end};

% Extract scan grid
[xValsV, yValsV, zValsV] = getScanXYZVals(planC{indexS.scan}(1));
xyzGridC = {xValsV, yValsV, zValsV};

% Extract scan
scan3M = double(planC{indexS.scan}(scanNum).scanArray);
ctOffset = double(planC{indexS.scan}(scanNum).scanInfo(1).CTOffset);
scan3M = scan3M - ctOffset;

% Extract Dose on scan grid
scanType = 'normal';
for iDose = 1:length(doseNumV)
    doseNum = doseNumV(iDose);
    dose3mC{iDose} = getDoseOnCT(doseNum, scanNum, scanType, planC);
end

% Extract structure
numStructs = length(structNamC);
strC = {planC{indexS.structures}.structureName};
strMaskC = cell(1,numStructs);
strColorC = cell(1,numStructs);
for iStr = 1:numStructs
    strNum = getMatchingIndex(structNamC{iStr},strC,'exact');
    strMask3M = getStrMask(strNum,planC);
    strMaskC{iStr} = strMask3M;
    strColorC{iStr} = planC{indexS.structures}(strNum).structureColor;
end


text:
function F = spm_Fcdf(x,v,w)
% Cumulative Distribution Function (CDF) of F (Fisher-Snedecor) distribution
% FORMAT F = spm_Fcdf(x,df)
% FORMAT F = spm_Fcdf(x,v,w)
%
% x  - F-variate   (F has range [0,Inf) )
% df - Degrees of freedom, concatenated along last dimension
%      Eg. Scalar (or column vector) v & w. Then df=[v,w];
% v  - Shape parameter 1 /   numerator degrees of freedom (v>0)
% w  - Shape parameter 2 / denominator degrees of freedom (w>0)
% F  - CDF of F-distribution with [v,w] degrees of freedom at points x
%__________________________________________________________________________
%
% spm_Fcdf implements the Cumulative Distribution Function of the F-distribution.
%
% Definition:
%--------------------------------------------------------------------------
% The CDF F(x) of the F distribution with degrees of freedom v & w,
% defined for positive integer degrees of freedom v & w, is the
% probability that a realisation of an F random variable X has value
% less than x F(x)=Pr{X<x} for X~F(v,w). The F-distribution is defined
% for v>0 & w>0, and for x in [0,Inf) (See Evans et al., Ch16).
%
% Variate relationships: (Evans et al., Ch16 & 37)
%--------------------------------------------------------------------------
% The square of a Student's t variate with w degrees of freedom is
% distributed as an F-distribution with [1,w] degrees of freedom.
%
% For X an F-variate with v,w degrees of freedom, w/(w+v*X^2) has
% distributed related to a Beta random variable with shape parameters
% w/2 & v/2.
%
% Algorithm:
%--------------------------------------------------------------------------
% Using the relationship with the Beta distribution: The CDF of the
% F-distribution with v,w degrees of freedom is related to the
% incomplete beta function by:
%       Pr(X<x) = 1 - betainc(w/(w+v*x^2),w/2,v/2)
% See Abramowitz & Stegun, 26.6.2; Press et al., Sec6.4 for
% definitions of the incomplete beta function. The relationship is
% easily verified by substituting for w/(w+v*x^2) in the integral of the
% incomplete beta function.
%
% MATLAB's implementation of the incomplete beta function is used.
%
%
% References:
%--------------------------------------------------------------------------
% Evans M, Hastings N, Peacock B (1993)
%       "Statistical Distributions"
%        2nd Ed. Wiley, New York
%
% Abramowitz M, Stegun IA, (1964)
%       "Handbook of Mathematical Functions"
%        US Government Printing Office
%
% Press WH, Teukolsky SA, Vetterling AT, Flannery BP (1992)
%       "Numerical Recipes in C"
%        Cambridge
%
%__________________________________________________________________________
% Copyright (C) 1992-2013 Wellcome Trust Centre for Neuroimaging

% Andrew Holmes
% $Id: spm_Fcdf.m 5602 2013-08-12 13:35:52Z guillaume $

%-Format arguments, note & check sizes
%--------------------------------------------------------------------------
if nargin<2, error('Insufficient arguments'), end

%-Unpack degrees of freedom v & w from single df parameter (v)
if nargin<3
    vs = size(v);
    if prod(vs)==2
        %-DF is a 2-vector
        w = v(2); v = v(1);
    elseif vs(end)==2
        %-DF has last dimension 2 - unpack v & w
        nv = prod(vs);
        w  = reshape(v(nv/2+1:nv),vs(1:end-1));
        v  = reshape(v(1:nv/2)   ,vs(1:end-1));
    else
        error('Can''t unpack both df components from single argument')
    end
end

%-Check argument sizes
ad = [ndims(x);ndims(v);ndims(w)];
rd = max(ad);
as = [[size(x),ones(1,rd-ad(1))];...
      [size(v),ones(1,rd-ad(2))];...
      [size(w),ones(1,rd-ad(3))]];
rs = max(as);
xa = prod(as,2)>1;
if sum(xa)>1 && any(any(diff(as(xa,:)),1))
    error('non-scalar args must match in size'), end

%-Computation
%--------------------------------------------------------------------------
%-Initialise result to zeros
F = zeros(rs);

%-Only defined for strictly positive v & w. Return NaN if undefined.
md = ( ones(size(x))  &  v>0  &  w>0 );
if any(~md(:))
    F(~md) = NaN;
    warning('Returning NaN for out of range arguments');
end

%-Non-zero where defined and x>0
Q  = find( md  &  x>0 );
if isempty(Q), return, end
if xa(1), Qx=Q; else Qx=1; end
if xa(2), Qv=Q; else Qv=1; end
if xa(3), Qw=Q; else Qw=1; end

%-Compute
F(Q) = 1 - betainc(w(Qw)./(w(Qw) + v(Qv).*x(Qx)),w(Qw)/2,v(Qv)/2);

author: spm


text:
%DEMGLM1 Demonstrate simple classification using a generalized linear model.
%
%	Description
%	 The problem consists of a two dimensional input matrix DATA and a
%	vector of classifications T.  The data is  generated from two
%	Gaussian clusters, and a generalized linear model with logistic
%	output is trained using iterative reweighted least squares. A plot of
%	the data together with the 0.1, 0.5 and 0.9 contour lines of the
%	conditional probability is generated.
%
%	See also
%	DEMGLM2, GLM, GLMTRAIN
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Generate data from two classes in 2d
input_dim = 2;

% Fix seeds for reproducible results
randn('state', 42);
rand('state', 42);

ndata = 100;
% Generate mixture of two Gaussians in two dimensional space
mix = gmm(2, 2, 'spherical');
mix.priors = [0.4 0.6];              % Cluster priors
mix.centres = [2.0, 2.0; 0.0, 0.0];  % Cluster centres
mix.covars = [0.5, 1.0];

[data, label] = gmmsamp(mix, ndata);
targets = label - ones(ndata, 1);

% Plot the result

clc
disp('This demonstration illustrates the use of a generalized linear model')
disp('to classify data from two classes in a two-dimensional space. We')
disp('begin by generating and plotting the data.')
disp(' ')
disp('Press any key to continue.')
pause

fh1 = figure;
plot(data(label==1,1), data(label==1,2), 'bo');
hold on
axis([-4 5 -4 5])
set(gca, 'box', 'on')
plot(data(label==2,1), data(label==2,2), 'rx')
title('Data')

clc
disp('Now we fit a model consisting of a logistic sigmoid function of')
disp('a linear combination of the input variables.')
disp(' ')
disp('The model is trained using the IRLS algorithm for 5 iterations')
disp(' ')
disp('Press any key to continue.')
pause

net = glm(input_dim, 1, 'logistic');
options = foptions;
options(1) = 1;
options(14) = 5;
net = glmtrain(net, options, data, targets);

disp(' ')
disp('We now plot some density contours given by this model.')
disp('The contour labelled 0.5 is the decision boundary.')
disp(' ')
disp('Press any key to continue.')
pause
x = -4.0:0.2:5.0;
y = -4.0:0.2:5.0;
[X, Y] = meshgrid(x,y);
X = X(:);
Y = Y(:);
grid = [X Y];
Z = glmfwd(net, grid);
Z = reshape(Z, length(x), length(y));
v = [0.1 0.5 0.9];
[c, h] = contour(x, y, Z, v);
title('Generalized Linear Model')
set(h, 'linewidth', 3)
clabel(c, h);

clc
disp('Note that the contours of constant density are straight lines.')
disp(' ')
disp('Press any key to end.')
pause
close(fh1);
clear all;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function need = inits ( os, nos, eta )

%*****************************************************************************80
%
%% INITS finds the number of Chebyshev terms needed to achieve a given accuracy.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2009
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    FORTRAN90 version by John Burkardt.
%
%  Parameters:
%
%    Input, real OS(NOS), the array of coefficients.
%
%    Input, integer NOS, the number of coefficients.
%
%    Input, real ETA, the requested accuracy.
%    A typical value of ETA is ( EPSILON ( 1.0 ) ) / 10.0.
%
%    Output, integer NEED, the number of terms needed for
%    the given accuracy.
%
  if ( nos < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'INITS - Fatal error!\n' );
    fprintf ( 1, '  The number of coefficients is less than 1.\n' );
    fprintf ( 1, '  NOS = %d\n', nos );
    error ( 'INITS - Fatal error!' );
  end

  need = nos;

  err = 0.0;
  for i = nos : -1 : 1
    err = err + abs ( os(i) );
    if ( eta < err )
      need = i;
      return
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'INITS - Warning!\n' );
  fprintf ( 1, '  The requested accuracy, ETA = %e\n', eta );
  fprintf ( 1, '  is smaller than the Chebyshev coefficients\n' );
  fprintf ( 1, '  can guarantee.\n' );

  return
end


text:
function dist = distancePoints3d(p1, p2, varargin)
%DISTANCEPOINTS3D Compute euclidean distance between pairs of 3D Points
%
%   D = distancePoints3d(P1, P2) return distance between points P1 and
%   P2, given as [X Y Z].
%
%   If P1 and P2 are two arrays of points, result is a N1*N2 array
%   containing distance between each point of P1 and each point of P2.
%
%
%   D = distancePoints3d(P1, P2, NOR)
%   with NOR being 1, 2, or Inf, corresponfing to the norm used. Default is
%   2 (euclidean norm). 1 correspond to manhattan (or taxi driver) distance
%   and Inf to maximum difference in each coordinate.
%
%
%   See also:
%   points3d, minDistancePoints, distancePoints
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 18/02/2005.
%

%   HISTORY
%   21/02/2005: add different norms
%   28/08/2007: deprecate

norm = 2;
if length(varargin)==1
    norm = varargin{1};
end

% compute difference of coordinate for each pair of points
ptsDiff = bsxfun(@minus, p2, p1);

% Return dist based on the type of measurement requested
switch(norm)
    case 1
        dist = sum(abs(ptsDiff),2);
    case 2
        dist = vectorNorm3d(ptsDiff);
    case Inf
        dist = max(abs(ptsDiff), [], 2);
    otherwise
        dist = power(sum(power(ptsDiff, norm),2), 1/norm);
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
%% EBSD Simulation
%
%%
% There is an implementation by Rüdiger Kilian ...

author: mtex-toolbox

repo: mtex

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


text:
function v = theta(v)
%
% Gives the vector in the tangential plane in v in the direction theta
%

if nargin == 1
  v = vector3d(cos(v.rho).*cos(v.theta), sin(v.rho).*cos(v.theta), -sin(v.theta));
else
  v = zvector;
end

end


text:
function [y1 y2 y3] = hrate(a,vad)

y1 =  mean(vad(a==1)); %hit rate speech

y2 = 1-mean(vad(a==0)); %hit rate silence

y3 = sum(abs(a-vad))/length(a);

author: Sable


text:
function val = get( f, propName )
%GET GET method for BALLFUN class.
%   P = GET(F, PROP) returns the property P specified in the string PROP from
%   the BALLFUN F. Valid entries for the string PROP are:
%    'coeffs', 'isReal', 'domain'

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

% Get the properties.
switch ( propName )
    case 'coeffs'
        val = f.coeffs;
    case 'isReal'
        val = f.isReal;
    case 'domain'
        val = f.domain;
otherwise
        error('CHEBFUN:BALLFUN:get:propName', ...
            [propName,' is not a valid BALLFUN property.'])
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function test_ft_volumelookup

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_volumelookup ft_read_atlas atlas_lookup

atlasfilename = dccnpath('/home/common/matlab/fieldtrip/template/atlas/afni/TTatlas+tlrc.BRIK');
mrifilename   = dccnpath('/home/common/matlab/fieldtrip/external/spm8/templates/T1.nii');

mri = ft_read_mri(mrifilename);
mri.coordsys = 'mni';

cfg       = [];
cfg.atlas = atlasfilename;
cfg.roi   = {'Precentral Gyrus'};
mask1     = ft_volumelookup(cfg,mri);
cfg.roi   = {'Brodmann area 47'};
mask2     = ft_volumelookup(cfg,mri);
cfg.roi   = {'Precentral Gyrus';'Brodmann area 47'};
mask3     = ft_volumelookup(cfg,mri);

assert(isequal(sum(mask1(:)),8144));%8213));
assert(isequal(sum(mask2(:)),2096));%2093));
assert(isequal(sum(mask3(:)&(mask2(:)|mask1(:))),sum(mask3(:))));

% the following is just to check whether the functionality does not crash
% no idea who would use it anyway
cfg       = [];
cfg.roi   = [15 30 40;-40 -20 80];
cfg.sphere = [10 8];
mask4     = ft_volumelookup(cfg,mri);

mri.mask  = mask4;

cfg       = [];
cfg.atlas = atlasfilename;
cfg.maskparameter = 'mask';
mask5     = ft_volumelookup(cfg, mri);

atlasfilename = dccnpath('/home/common/matlab/fieldtrip/template/atlas/aal/ROI_MNI_V4.nii');
cfg       = [];
cfg.atlas = atlasfilename;
cfg.roi = 'Calcarine_R'; % right V1
mask6 = ft_volumelookup(cfg, mri);
assert(isequal(sum(mask6(:)),1861));

atlas_MNI = ft_read_atlas(atlasfilename);
atlas_MNI.coordsys = 'mni';
cfg = [];
cfg.roi = [52 -9 -45; 73 -37 -8];
cfg.atlas = atlasfilename;
cfg.output = 'single';
cfg.maxqueryrange = 29;
label_MNI = ft_volumelookup(cfg, atlas_MNI);
assert(size(label_MNI, 1) == 1);
cfg.output = 'label';
label_MNI = ft_volumelookup(cfg, atlas_MNI);
assert(size(label_MNI, 1) == 1);
cfg.output = 'multiple';
label_MNI = ft_volumelookup(cfg, atlas_MNI);
assert(size(label_MNI, 1) == 2);

author: fieldtrip

repo: fieldtrip


text:
function path = imGeodesicPath(img, source, target, varargin)
%IMGEODESICPATH Compute a geodesic path between two markers in an image
%
%   PATH = imGeodesicPath(MASK, SOURCE, TARGET)
%   Computea minimal geodesic-distance path between the two markers SOURCE
%   and TARGET. Both SOURCE and TARGET can be either a binary image the
%   same size as MASK, or a N-by-2 array containing a list of coordinates
%   (x first, y second).
%   The result PATH is a P-by-2 array containing coordinates of a polyline
%   with minimal geodesic length starting from one of the points specified
%   by SOURCE, and terminating at one of the points specified by TARGET.
%   Note that the result is not uniquely defined, and the returned solution
%   is one of the possible solutions.
%
%   PATH = imGeodesicPath(..., WEIGHTS)
%   Specify the weights to use for propagating chamfer distance. Default is
%   [3 4], as suggested by Borgefors.
%
%   PATH = imGeodesicPath(..., 'verbose', V)
%   Specify the verbosity. V can be either TRUE or FALSE.
%
%
%   Example
%     % read circle image, and create 2 markers
%     img = imread('circles.png');
%     imshow(img); hold on;
%     p1 = [130 130]; % (x1,y1)
%     p2 = [170 170]; % (x2,y2)
%     plot(p1(1), p1(2), 'bo');
%     plot(p2(1), p2(2), 'ro');
%     % Compute and display the path as a polyline
%     path = imGeodesicPath(img, p1, p2);
%     plot(path(:,1), path(:,2), 'color', 'm', 'linewidth', 2);
%
%   See also
%   imGeodesics, imGeodesicDistanceMap, imMaxGeodesicPath

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2011-04-04,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

%% Default values

% weights for computing geodesic lengths
ws = [3 4];

% no verbosity by default
verbose = 0;

%% process input arguments

% extract weights if present
if ~isempty(varargin)
    if isnumeric(varargin{1})
        ws = varargin{1};
        varargin(1) = [];
    end
end

% Extract options
while ~isempty(varargin)
    paramName = varargin{1};
    if strcmpi(paramName, 'verbose')
        verbose = varargin{2};
    else
        error(['Unkown option in imGeodesicLength: ' paramName]);
    end
    varargin(1:2) = [];
end

% forces input image to binary
img = img > 0;

%% Propagate distance from destination markers

% initialize marker image
if sum(size(target) == size(img)) == 2
    % initialize from binary mask
    markers = target > 0;

else
    % initialize from set of points
    markers = false(size(img));
    for i = 1:size(target, 1)
        markers(target(i,2), target(i,1)) = true;
    end
end

% compute distance map from marker
dist = imGeodesicDistanceMap(img, markers, ws, 'verbose', verbose);

%% Find position of first source point

% find position of closest point belonging to the source
if sum(size(source) == size(img)) == 2
    % compute in a binary mask
    [minDist, ind] = min(dist(source));
    [ys, xs] = ind2sub(size(source), ind);

else
    % initialize from set of points
    minDist = inf;
    for i = 1:size(source, 1)
        value = dist(source(i,2), source(i,1));
        if value < minDist
            xs = source(i,1);
            ys = source(i,2);
            minDist = value;
        end
    end
end

% check existence of path
if isempty(minDist)
    warning([mfilename ':NoPathFound'], ...
        'No path could be found between the two markers');
    path = [];
    return;
end

%% Create the path by returning to destination marker

% add a 1-pixel wide border around image, initialized with very high value
tmp = dist;
dist = zeros(size(dist)+2, class(dist));
dist(:) = max(tmp(:)) + 2;
dist(2:end-1, 2:end-1) = tmp;

% keep coordinates of geodesic extremity
x = xs + 1;
y = ys + 1;

% initialize path
path = [x y];

while true
    % check around current point
    neigh = dist(y-1:y+1, x-1:x+1);

    % look for the minimum
    [mini, ind] = min(neigh(:));

    % if minimum is the same as current value, minima is found
    if mini == dist(y, x)
        break;
    end

    % convert index to sub
    [iy, ix] = ind2sub([3 3], ind(1));

    % update coord
    x = x + ix - 2;
    y = y + iy - 2;

    % stores result
    path = [path; x y]; %#ok<AGROW>
end

% subtract 1, because of border
path = path - 1;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage