text:
function para = RemoveIOStream(para, streamIdxToBeRemoved)

nStreamRemove = length(streamIdxToBeRemoved);

para.nStream = para.nStream - nStreamRemove;

ArrayConfig = GetArrayConfig();

for i=1:length(ArrayConfig)
    if isfield(para, ArrayConfig{i}) && length(para.(ArrayConfig{i})) >= max(streamIdxToBeRemoved)
        para.(ArrayConfig{i})(streamIdxToBeRemoved) = [];
    end
end

if isfield(para, 'DataSyncSet')
    for i=1:length(para.DataSyncSet)
        if length(para.DataSyncSet{i})>1
            for j=1:length(streamIdxToBeRemoved)
                idx = para.DataSyncSet{i}==streamIdxToBeRemoved(j);
                if sum(idx)
                    para.DataSyncSet{i}(idx) = [];  % remove the stream from DataSyncSet if it is already removed
                end
            end
        end
        if length(para.DataSyncSet{i})<2    % if there is less than 2 streams, no need to sync
            para.DataSyncSet{i} = [];
        end
    end
end

if isfield(para, 'ApplyVADSet')
    for i=1:length(para.ApplyVADSet)
        if length(para.ApplyVADSet{i})>1
            for j=1:length(streamIdxToBeRemoved)
                idx = para.ApplyVADSet{i}==streamIdxToBeRemoved(j);
                if sum(idx)
                    para.ApplyVADSet{i}(idx) = [];  % remove the stream from DataSyncSet if it is already removed
                end
            end
        end
        if length(para.ApplyVADSet{i})<2    % if there is less than 2 streams, no need to sync
            para.ApplyVADSet{i} = [];
        end
    end
end

end

author: singaxiong

repo: SignalGraph


text:
%% mission_data.param4 contains some NaN data, which is invalid for timeserials.
%% so remove it as we are not using param4 yet.

[row, col] = size(Mission_Data.param4.Data);

for i = 1 : row * col
    if isnan(Mission_Data.param4.Data(i))
       Mission_Data.param4.Data(i) = single(0);
    end
end

author: Firmament-Autopilot

repo: FMT-Model


text:
s = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789';

hash = java.util.Hashtable;
C = {};
for i = 1:100
    key = randsample(s,20,true);
    value = randsample(s,40,true);
    hash.put(key, value);
    C{end+1} = {key, value};
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
% rigid registration of frames with offsets ds
function [dreg] = rigidMovie(data, ops1, dsall, yFOVs, xFOVs)

ix0 = 0;
Nbatch = 1000;
dreg = zeros(size(data), 'single');
while ix0<size(data,3)
    indxr = ix0 + (1:Nbatch);
    indxr(indxr>size(data,3)) = [];
    for l = 1:size(xFOVs,2)
        dreg(yFOVs(:,l), xFOVs(:,l), indxr)        = ...
            rigidRegFrames(data(yFOVs(:,l), xFOVs(:,l), indxr), ops1{1,l}, dsall(indxr,:,l));
    end
    ix0 = ix0 + Nbatch;
end

author: cortex-lab


text:
function [x, cost, info, options] = rlbfgs(problem, x0, options)
% Riemannian limited memory BFGS solver for smooth objective functions.
%
% function [x, cost, info, options] = rlbfgs(problem)
% function [x, cost, info, options] = rlbfgs(problem, x0)
% function [x, cost, info, options] = rlbfgs(problem, x0, options)
% function [x, cost, info, options] = rlbfgs(problem, [], options)
%
%
% This is a Riemannian limited memory BFGS solver (quasi-Newton method),
% which aims to minimize the cost function in the given problem structure.
% It requires access to the gradient of the cost function.
%
% Parameter options.memory can be used to specify the number of iterations
% the algorithm remembers and uses to approximate the inverse Hessian of
% the cost. Default value is 30.
% For unlimited memory, set options.memory = Inf.
%
%
% For a description of the algorithm and theorems offering convergence
% guarantees, see the references below.
%
% The initial iterate is x0 if it is provided. Otherwise, a random point on
% the manifold is picked. To specify options whilst not specifying an
% initial iterate, give x0 as [] (the empty matrix).
%
% The two outputs 'x' and 'cost' are the last reached point on the manifold
% and its cost.
%
% The output 'info' is a struct-array which contains information about the
% iterations:
%   iter (integer)
%       The iteration number. The initial guess is 0.
%   cost (double)
%       The corresponding cost value.
%   gradnorm (double)
%       The (Riemannian) norm of the gradient.
%   time (double)
%       The total elapsed time in seconds to reach the corresponding cost.
%   stepsize (double)
%       The size of the step from the previous to the new iterate.
%   accepted (Boolean)
%       true if step is accepted in the cautious update. 0 otherwise.
%   And possibly additional information logged by options.statsfun.
% For example, type [info.gradnorm] to obtain a vector of the successive
% gradient norms reached at each iteration.
%
% The options structure is used to overwrite the default values. All
% options have a default value and are hence optional. To force an option
% value, pass an options structure with a field options.optionname, where
% optionname is one of the following and the default value is indicated
% between parentheses:
%
%   tolgradnorm (1e-6)
%       The algorithm terminates if the norm of the gradient drops below
%       this. For well-scaled problems, a rule of thumb is that you can
%       expect to reduce the gradient norm by 8 orders of magnitude
%       (sqrt(eps)) compared to the gradient norm at a "typical" point (a
%       rough initial iterate for example). Further decrease is sometimes
%       possible, but inexact floating point arithmetic will eventually
%       limit the final accuracy. If tolgradnorm is set too low, the
%       algorithm may end up iterating forever (or at least until another
%       stopping criterion triggers).
%   maxiter (1000)
%       The algorithm terminates if maxiter iterations were executed.
%   maxtime (Inf)
%       The algorithm terminates if maxtime seconds elapsed.
%   minstepsize (1e-10)
%     The minimum norm of the tangent vector that points from the current
%     point to the next point. If the norm is less than minstepsize, the
%     program will terminate.
%   memory (30)
%     The number of previous iterations the program remembers. This is used
%     to approximate the inverse Hessian at the current point. Because of
%     difficulty of maintaining a representation of operators in terms of
%     coordinates, a recursive method is used. The number of steps in the
%     recursion is at most options.memory. This parameter can take any
%     integer value >= 0, or Inf, which is taken to be options.maxiter. If
%     options.maxiter has value Inf, then it will take value 10000 and a
%     warning will be displayed.
%   strict_inc_func (@(t) 1e-4*t)
%     The Cautious step needs a real function that has value 0 at t = 0,
%     and  is strictly increasing. See details in Wen Huang's paper
%     "A Riemannian BFGS Method without Differentiated Retraction for
%     Nonconvex Optimization Problems"
%   statsfun (none)
%       Function handle to a function that will be called after each
%       iteration to provide the opportunity to log additional statistics.
%       They will be returned in the info struct. See the generic Manopt
%       documentation about solvers for further information. statsfun is
%       called with the point x that was reached last.
%   stopfun (none)
%       Function handle to a function that will be called at each iteration
%       to provide the opportunity to specify additional stopping criteria.
%       See the generic Manopt documentation about solvers for further
%       information.
%   verbosity (2)
%       Integer number used to tune the amount of output the algorithm
%       generates during execution (mostly as text in the command window).
%       The higher, the more output. 0 means silent. 3 and above includes a
%       display of the options structure at the beginning of the execution.
%   debug (false)
%       Set to true to allow the algorithm to perform additional
%       computations for debugging purposes. If a debugging test fails, you
%       will be informed of it, usually via the command window. Be aware
%       that these additional computations appear in the algorithm timings
%       too, and may interfere with operations such as counting the number
%       of cost evaluations, etc. (the debug calls get storedb too).
%   storedepth (2)
%       Maximum number of different points x of the manifold for which a
%       store structure will be kept in memory in the storedb for caching.
%       If memory usage is an issue, you may try to lower this number.
%       Profiling may then help to investigate if a performance hit was
%       incurred as a result.
%   ls_initial_scale (@(gradnorm) 1/gradnorm)
%       A function handle that takes as input a real number (a gradient norm)
%       and outputs a real number for how to scale the initialization of
%       line-search.
%
%
% Please cite the Manopt paper as well as the research paper:
% @InBook{Huang2016,
%   title     = {A {R}iemannian {BFGS} Method for Nonconvex Optimization Problems},
%   author    = {Huang, W. and Absil, P.-A. and Gallivan, K.A.},
%   year      = {2016},
%   publisher = {Springer International Publishing},
%   editor    = {Karas{\"o}zen, B{\"u}lent and Manguo{\u{g}}lu, Murat and Tezer-Sezgin, M{\"u}nevver and G{\"o}ktepe, Serdar and U{\u{g}}ur, {\"O}m{\"u}r},
%   address   = {Cham},
%   booktitle = {Numerical Mathematics and Advanced Applications ENUMATH 2015},
%   pages     = {627--634},
%   doi       = {10.1007/978-3-319-39929-4_60}
% }
%
% We point out that, at the moment, this implementation of RLBFGS can be
% slower than the implementation in ROPTLIB by Wen Huang et al. referenced
% above. For the purpose of comparing to their work, please use their
% implementation.

% This file is part of Manopt: www.manopt.org.
% Original author: Changshuo Liu, July 19, 2017.
% Contributors: Nicolas Boumal
% Change log:
%
%   Nov. 27, 2017 (Wen Huang):
%       Changed the default strict_inc_func to @(t) 1e-4*t from @(t) t.
%
%   Jan. 18, 2018 (NB):
%       Corrected a bug related to the way the line search hint was defined
%       by default.
%
%   Aug. 2, 2018 (NB):
%       Using the new storedb.remove features to keep storedb lean, and
%       reduced the default value of storedepth from 30 to 2 as a result.
%
%   Aug. 2, 2022 (sfrcorne):
%       Added option ls_initial_scale, with default setting to ensure
%       the method is invariant to positive scaling of the cost function.

    % Verify that the problem description is sufficient for the solver.
    if ~canGetCost(problem)
        warning('manopt:getCost', ...
                'No cost provided. The algorithm will likely abort.');
    end
    if ~canGetGradient(problem) && ~canGetApproxGradient(problem)
        % Note: we do not give a warning if an approximate gradient is
        % explicitly given in the problem description, as in that case the user
        % seems to be aware of the issue.
        warning('manopt:getGradient:approx', ...
           ['No gradient provided. Using an FD approximation instead (slow).\n' ...
            'It may be necessary to increase options.tolgradnorm.\n' ...
            'To disable this warning: warning(''off'', ''manopt:getGradient:approx'')']);
        problem.approxgrad = approxgradientFD(problem);
    end

    % This solver uses linesearch_hint as a line search algorithm. By
    % default, try a step size of 1, so that if the BFGS approximation of
    % the Hessian (or inverse Hessian) is good, then the iteration is close
    % to a Newton step.
    if ~canGetLinesearch(problem)
        problem.linesearch = @(x, d) 1;
    end

    % Local defaults for the program
    localdefaults.minstepsize = 1e-10;
    localdefaults.maxiter = 1000;
    localdefaults.tolgradnorm = 1e-6;
    localdefaults.memory = 30;
    localdefaults.strict_inc_func = @(t) 1e-4*t;
    localdefaults.ls_max_steps = 25;
    localdefaults.ls_initial_scale = @(gradnorm) 1/gradnorm;
    localdefaults.storedepth = 2;

    % Merge global and local defaults, then merge w/ user options, if any.
    localdefaults = mergeOptions(getGlobalDefaults(), localdefaults);
    if ~exist('options', 'var') || isempty(options)
        options = struct();
    end
    options = mergeOptions(localdefaults, options);

    % To make sure memory in range [0, Inf)
    options.memory = max(options.memory, 0);
    if options.memory == Inf
        if isinf(options.maxiter)
            options.memory = 10000;
            warning('rlbfgs:memory', ['options.memory and options.maxiter' ...
              ' are both Inf; options.memory has been changed to 10000.']);
        else
            options.memory = options.maxiter;
        end
    end

    M = problem.M;

    timetic = tic();

    % Create a random starting point if no starting point is provided.
    if ~exist('x0', 'var')|| isempty(x0)
        xCur = M.rand();
    else
        xCur = x0;
    end

    % Create a store database and get a key for the current x
    storedb = StoreDB(options.storedepth);
    key = storedb.getNewKey();

    % __________Initialization of variables______________
    % Number of iterations since the last restart
    k = 0;
    % Total number of BFGS iterations
    iter = 0;

    % This cell stores step vectors which point from x_{t} to x_{t+1} for t
    % indexing the last iterations, capped at options.memory.
    % That is, it stores up to options.memory of the most recent step
    % vectors. However, the implementation below does not need step vectors
    % in their respective tangent spaces at x_{t}'s. Rather, it requires
    % them transported to the current point's tangent space by vector
    % transport. For details regarding the requirements on the the vector
    % transport, see the reference paper by Huang et al.
    % In this implementation, those step vectors are iteratively
    % transported to the current point's tangent space after every
    % iteration. Thus, at every iteration, vectors in sHistory are in the
    % current point's tangent space.
    sHistory = cell(1, options.memory);

    % This cell stores the differences for latest t's of the gradient at
    % x_{t+1} and the gradient at x_{t}, transported to x_{t+1}'s tangent
    % space. The memory is also capped at options.memory.
    yHistory = cell(1, options.memory);

    % rhoHistory{t} stores the reciprocal of the inner product between
    % sHistory{t} and yHistory{t}.
    rhoHistory = cell(1, options.memory);

    % Scaling of direction given by getDirection for acceptable step
    alpha = 1;

    % Norm of the step
    stepsize = 1;

    % Stores whether the step is accepted by the cautious update check.
    accepted = true;

    % Query the cost function and its gradient
    [xCurCost, xCurGradient] = getCostGrad(problem, xCur, storedb, key);

    xCurGradNorm = M.norm(xCur, xCurGradient);

    % Scaling of initial matrix, Barzilai-Borwein.
    scaleFactor = options.ls_initial_scale(xCurGradNorm);

    % Line-search statistics for recording in info.
    lsstats = [];

    % Flag to control restarting scheme to avoid infinite loops (see below)
    ultimatum = false;

    % Save stats in a struct array info, and preallocate.
    stats = savestats();
    info(1) = stats;
    info(min(10000, options.maxiter+1)).iter = [];

    if options.verbosity >= 2
        fprintf(' iter                   cost val            grad. norm           alpha\n');
    end

    % Main iteration
    while true

        % Display iteration information
        if options.verbosity >= 2
        fprintf('%5d    %+.16e        %.8e      %.4e\n', ...
                iter, xCurCost, xCurGradNorm, alpha);
        end

        % Start timing this iteration
        timetic = tic();

        % Run standard stopping criterion checks
        [stop, reason] = stoppingcriterion(problem, xCur, options, ...
                                           info, iter+1);

        % If none triggered, run specific stopping criterion check
        if ~stop
            if stats.stepsize < options.minstepsize
                % To avoid infinite loop and to push the search further
                % in case BFGS approximation of Hessian is off towards
                % the end, we erase the memory by setting k = 0;
                % In this way, it starts off like a steepest descent.
                % If even steepest descent does not work, then it is
                % hopeless and we will terminate.
                if ~ultimatum
                    if options.verbosity >= 2
                        fprintf(['stepsize is too small, restarting ' ...
                            'the bfgs procedure at the current point.\n']);
                    end
                    k = 0;
                    ultimatum = true;
                else
                    stop = true;
                    reason = sprintf(['Last stepsize smaller than '  ...
                        'minimum allowed; options.minstepsize = %g.'], ...
                        options.minstepsize);
                end
            else
                % We are not in trouble: lift the ultimatum if it was on.
                ultimatum = false;
            end
        end

        if stop
            if options.verbosity >= 1
                fprintf([reason '\n']);
            end
            break;
        end

        % Compute BFGS direction
        p = getDirection(M, xCur, xCurGradient, sHistory,...
                yHistory, rhoHistory, scaleFactor, min(k, options.memory));

        % Execute line-search
        [stepsize, xNext, newkey, lsstats] = ...
            linesearch_hint(problem, xCur, p, xCurCost, ...
                            M.inner(xCur, xCurGradient, p), ...
                            options, storedb, key);

        % Record the BFGS step-multiplier alpha which was effectively
        % selected. Toward convergence, we hope to see alpha = 1.
        alpha = stepsize/M.norm(xCur, p);
        step = M.lincomb(xCur, alpha, p);

        % Query cost and gradient at the candidate new point.
        [xNextCost, xNextGrad] = getCostGrad(problem, xNext, storedb, newkey);

        % Compute sk and yk
        sk = M.transp(xCur, xNext, step);
        yk = M.lincomb(xNext, 1, xNextGrad, ...
                             -1, M.transp(xCur, xNext, xCurGradient));

        % Computation of the BFGS step is invariant under scaling of sk and
        % yk by a common factor. For numerical reasons, we scale sk and yk
        % so that sk is a unit norm vector.
        norm_sk = M.norm(xNext, sk);
        sk = M.lincomb(xNext, 1/norm_sk, sk);
        yk = M.lincomb(xNext, 1/norm_sk, yk);

        inner_sk_yk = M.inner(xNext, sk, yk);
        inner_sk_sk = M.norm(xNext, sk)^2;    % ensures nonnegativity

        % If the cautious step is accepted (which is the intended
        % behavior), we record sk, yk and rhok and need to do some
        % housekeeping. If the cautious step is rejected, these are not
        % recorded. In all cases, xNext is the next iterate: the notion of
        % accept/reject here is limited to whether or not we keep track of
        % sk, yk, rhok to update the BFGS operator.
        cap = options.strict_inc_func(xCurGradNorm);
        if inner_sk_sk ~= 0 && (inner_sk_yk / inner_sk_sk) >= cap

            accepted = true;

            rhok = 1/inner_sk_yk;

            scaleFactor = inner_sk_yk / M.norm(xNext, yk)^2;

            % Time to store the vectors sk, yk and the scalar rhok.
            % Remember: we need to transport all vectors to the most
            % current tangent space.

            % If we are out of memory
            if k >= options.memory

                % sk and yk are saved from 1 to the end with the most
                % current recorded to the rightmost hand side of the cells
                % that are occupied. When memory is full, do a shift so
                % that the rightmost is earliest and replace it with the
                % most recent sk, yk.
                for  i = 2 : options.memory
                    sHistory{i} = M.transp(xCur, xNext, sHistory{i});
                    yHistory{i} = M.transp(xCur, xNext, yHistory{i});
                end
                if options.memory > 1
                    sHistory = sHistory([2:end, 1]);
                    yHistory = yHistory([2:end, 1]);
                    rhoHistory = rhoHistory([2:end 1]);
                end
                if options.memory > 0
                    sHistory{options.memory} = sk;
                    yHistory{options.memory} = yk;
                    rhoHistory{options.memory} = rhok;
                end

            % If we are not out of memory
            else

                for i = 1:k
                    sHistory{i} = M.transp(xCur, xNext, sHistory{i});
                    yHistory{i} = M.transp(xCur, xNext, yHistory{i});
                end
                sHistory{k+1} = sk;
                yHistory{k+1} = yk;
                rhoHistory{k+1} = rhok;

            end

            k = k + 1;

        % The cautious step is rejected: we do not store sk, yk, rhok but
        % we still need to transport stored vectors to the new tangent
        % space.
        else

            accepted = false;

            for  i = 1 : min(k, options.memory)
                sHistory{i} = M.transp(xCur, xNext, sHistory{i});
                yHistory{i} = M.transp(xCur, xNext, yHistory{i});
            end

        end

        % Update variables to new iterate.
        storedb.removefirstifdifferent(key, newkey);
        xCur = xNext;
        key = newkey;
        xCurGradient = xNextGrad;
        xCurGradNorm = M.norm(xNext, xNextGrad);
        xCurCost = xNextCost;

        % iter is the number of iterations we have accomplished.
        iter = iter + 1;

        % Make sure we don't use too much memory for the store database
        % (this is independent from the BFGS memory.)
        storedb.purge();

        % Log statistics for freshly executed iteration
        stats = savestats();
        info(iter+1) = stats;

    end

    % Housekeeping before we return
    info = info(1:iter+1);
    x = xCur;
    cost = xCurCost;

    if options.verbosity >= 1
        fprintf('Total time is %f [s] (excludes statsfun)\n', ...
                info(end).time);
    end

    % Routine in charge of collecting the current iteration stats
    function stats = savestats()
        stats.iter = iter;
        stats.cost = xCurCost;
        stats.gradnorm = xCurGradNorm;
        if iter == 0
            stats.stepsize = NaN;
            stats.time = toc(timetic);
            stats.accepted = NaN;
        else
            stats.stepsize = stepsize;
            stats.time = info(iter).time + toc(timetic);
            stats.accepted = accepted;
        end
        stats.linesearch = lsstats;
        stats = applyStatsfun(problem, xCur, storedb, key, options, stats);
    end

end

% BFGS step, see Wen's paper for details. This function takes in a tangent
% vector g, and applies an approximate inverse Hessian P to it to get Pg.
% Then, -Pg is returned.
%
% Theory requires the vector transport to be isometric and to satisfy the
% locking condition (see paper), but these properties do not seem to be
% crucial in practice. If your manifold provides M.isotransp, it may be
% good to do M.transp = M.isotransp; after loading M with a factory.
%
% This implementation operates in the tangent space of the most recent
% point since all vectors in sHistory and yHistory have been transported
% there.
function dir = getDirection(M, xCur, xCurGradient, sHistory, yHistory, ...
                            rhoHistory, scaleFactor, k)

    q = xCurGradient;

    inner_s_q = zeros(1, k);

    for i = k : -1 : 1
        inner_s_q(1, i) = rhoHistory{i} * M.inner(xCur, sHistory{i}, q);
        q = M.lincomb(xCur, 1, q, -inner_s_q(1, i), yHistory{i});
    end

    r = M.lincomb(xCur, scaleFactor, q);

    for i = 1 : k
         omega = rhoHistory{i} * M.inner(xCur, yHistory{i}, r);
         r = M.lincomb(xCur, 1, r, inner_s_q(1, i)-omega, sHistory{i});
    end

    dir = M.lincomb(xCur, -1, r);

end

author: NicolasBoumal

repo: manopt


text:
function [states] = tapas_sem_hier_init_states(data, model, inference)
%% Generates the data structure to store the samples.
%
% Input
%
%
% Output
%   data        -- Cell array with nx1 cells, where n = niter
%

% aponteeduardo@gmail.com
% copyright (C) 2016
%

states = cell(inference.niter, 1);

end

author: translationalneuromodeling

repo: tapas


text:
function [x,y,button] = ginput_zoom(n, style)
%  GINPUT_ZOOM Graphical input from mouse with Zoom In/Out function.
%     [X,Y] = GINPUT_ZOOM(N, STYLE) gets N points from the current axes and returns
%     the X- and Y-coordinates in length N vectors X and Y.  The cursor
%     can be positioned using a mouse (or by using the Arrow Keys on some
%     systems).  Data points are entered by pressing a mouse button
%     or any key on the keyboard except carriage return, which terminates
%     the input before N points are entered. The style of cursor can be
%     change by setting STYLE (please reference to "Specifying the Figure
%     Pointer" in the help, default: 'fullcross').
%
%     User can click and drag to Zoom In / Out
%     - Left-Top to Right-Down means Zoom In
%     - Right-Down to left-Top means Zoom Out
%
%     [X,Y] = GINPUT_ZOOM gathers an unlimited number (less then 100000)
%     of points until the return key is pressed.
%
%     [X,Y,BUTTON] = GINPUT_ZOOM(N) returns a third result, BUTTON, that
%     contains a vector of integers specifying which mouse button was
%     used (1,2,3 from left) or ASCII numbers if a key on the keyboard
%     was used.
%
%     Chang Yi-Chung 20091002

ox = xlim;     % get original x limit
oy = ylim;     % get original y limit
xlim('auto');  % zoom out to x max
ylim('auto');  % zoom out to y max
xmax = xlim;   % get max x
ymax = ylim;   % get max y
xlim(ox);      % restore x limit
ylim(oy);      % restore y limit

if exist('style') == 0
    style = 'fullcross';
end

set(gcf,'pointer', style);  % use fullcross, user can change this

if exist('n') == 0
    n = 100000;
end

for i = 1:n
    while 1
        k = waitforbuttonpress;
        point1 = get(gca,'CurrentPoint');    % button down detected
        finalRect = rbbox;                   % return figure units
        point2 = get(gca,'CurrentPoint');    % button up detected
        point1 = point1(1,1:2);              % extract x and y
        point2 = point2(1,1:2);
        if sum(point1 == point2) == 2        % same points
            if  k == 1                       % detects a key press
                x(i) = point1(1);
                y(i) = point1(2);
                button(i) = double(get(gcf,'CurrentCharacter'));   % get keyboad ascii number
                if n == 100000 & button(i) == 13                   % n is not assigned and press Enter(13) -> exit
                    i = n;
                    % compatible with ginput
                    x = x(1:end-1);
                    y = y(1:end-1);
                    button = button(1:end-1);
                end
            else                             % detects a mouse button press
                x(i) = point1(1);
                y(i) = point1(2);
                switch get(gcf,'SelectionType')
                    case 'normal'            % left button
                        button(i) = 1;
                    case 'extend'            % middle button
                        button(i) = 2;
                    case 'alt'               % right button
                        button(i) = 3;
                end
            end
            break;
        end
        offset = (point2-point1);
        if offset(1) > 0 & offset(2) < 0     % drag left-top to right-down -> zoom in
            xlim([point1(1) point2(1)]);
            ylim([point2(2) point1(2)]);
        end
        if offset(1) < 0 & offset(2) > 0     % drag right-down to left-top -> zoom out
            xlim(xmax);
            ylim(ymax);
        end
    end
    if i == n
        break;
    end
end

% compatible with ginput
x = x';
y = y';
button = button';

author: Sable

repo: mcbench-benchmarks


text:
function test007 ( )

%*****************************************************************************80
%
%% TEST007 tests TET_MESH_SEARCH_NAIVE and TET_MESH_SEARCH_DELAUNAY.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 August 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 5;
  tet_order = 4;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST007\n' );
  fprintf ( 1, '  TET_MESH_SEARCH_NAIVE uses a naive algorithm\n' );
  fprintf ( 1, '  to search a tetrahedral mesh for the tetrahedron\n' );
  fprintf ( 1, '  containing a point.\n' );
  fprintf ( 1, '  TET_MESH_SEARCH_DELAUNAY uses a faster algorithm\n' );
  fprintf ( 1, '  which is appropriate if the tet mesh is Delaunay.\n' );
%
%  Set up the example tetrahedron mesh.
%
  [ node_num, tet_num ] = tet_mesh_order4_example_size ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  This mesh has tetrahedron order %d\n', tet_order );
  fprintf ( 1, '  The number of tetrahedrons is   %d\n', tet_num );

  [ node_xyz, tet_node ] = tet_mesh_order4_example_set ( node_num, tet_num );
%
%  TET_NEIGHBOR is needed for the Delaunay search.
%
  tet_neighbor = tet_mesh_neighbor_tets ( tet_order, tet_num, tet_node );

  for test = 1 : test_num
%
%  Choose a tetrahedral index at random.
%
    [ tet1, seed ] = i4_uniform ( 1, tet_num, seed );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Point was chosen from tetrahedron    %8d\n', tet1 );
%
%  Choose a point at random from the tetrahedron.
%
    [ p, seed ] = tetrahedron_sample ( node_xyz(1:3,tet_node(1:4,tet1)), 1, seed );
%
%  Naive search.
%
    [ tet2, step_num ] = tet_mesh_search_naive ( node_num, node_xyz, ...
      tet_order, tet_num, tet_node, p );

    fprintf ( 1, '  Naive search ended in tetrahedron    %8d after %d steps\n', ...
      tet2, step_num );
%
%  Delaunay search.
%
    [ tet3, face, step_num ] = tet_mesh_search_delaunay ( node_num, ...
      node_xyz, tet_order, tet_num, tet_node, tet_neighbor, p );

    fprintf ( 1, ...
      '  Delaunay search ended in tetrahedron %8d after %d steps.\n', ...
      tet3, step_num );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function planC = createExpandedStructure2D(structNum, margin, planC)
% function createExpandedStructure2D(structNum, margin, planC)
%
% APA, 09/09/2014
% AI, 09/19/17 Updated to allow margin<0

if ~exist('planC','var')
    global planC
end

global stateS

indexS = planC{end};

scanNum = getStructureAssociatedScan(structNum,planC);

newStructNum = length(planC{indexS.structures}) + 1;

%newStructS = newCERRStructure(scanNum, planC, newStructNum);

rasterSegs = getRasterSegments(structNum, planC);

% Get Rastersegments +/- margin, thus creating an halo
if margin>0
halo = structMargin(rasterSegs, margin, scanNum, planC);
rasterSegsExpanded = structUnion(halo, rasterSegs, scanNum, planC); % Expand rastersegments by margin
else
halo = structMargin(rasterSegs, -margin, scanNum, planC);
rasterSegsExpanded = structDiff(rasterSegs, halo, scanNum, planC); % Shrink rastersegments by margin
end

% Get Contours from rasterSegs
contourS = rasterToPoly(rasterSegsExpanded, scanNum, planC);

% Generate Structure Name
strName = planC{indexS.structures}(structNum).structureName;
strName = [strName,' + 2D_',num2str(margin)];

% Create New structure
%Make an empty structure, assign name/contour.
newstr = newCERRStructure(scanNum, planC, newStructNum);
newstr.contour = contourS;
newstr.structureName = strName;
newstr.associatedScan = scanNum;
newstr.assocScanUID = planC{indexS.scan}(scanNum).scanUID;
numStructs = length(planC{indexS.structures});

%Append new structure to planC.
if ~isempty(planC{indexS.structures})
    planC{indexS.structures} = dissimilarInsert(planC{indexS.structures}, newstr, numStructs+1, []);
else
    planC{indexS.structures} = newstr;
end

%Update uniformized data.
planC = updateStructureMatrices(planC, numStructs+1);

% Refresh View
if ~isempty(stateS) && isfield(stateS,'handle') && isfield(stateS.handle,'CERRSliceViewer') && isnumeric(stateS.handle.CERRSliceViewer)
    stateS.structsChanged = 1;
    CERRRefresh
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function rtk=udion(rtk,tt,bl,ind)

global glc;
GAP_RESION=120;

for i=1:glc.MAXSAT
    j=rtk.ii+i;
    if rtk.x(j)~=0&&rtk.sat(i).outc(1)>GAP_RESION&&rtk.sat(i).outc(2)>GAP_RESION
        rtk.x(j)=0;
    end
end

for i=1:ind.ns
    j=rtk.ii+ind.sat(i);

    if rtk.x(j)==0
        rtk=initx(rtk,1e-6,(rtk.opt.std(2)*bl/1e4)^2,j);
    else
        el=rtk.sat(ind.sat(i)).azel(2);
        fact=cos(el);
        rtk.P(j,j)=rtk.P(j,j)+(rtk.opt.prn(2)*bl/1e4*fact)^2*abs(tt);
    end
end

return

repo: GINav


text:
function [X_norm, mu, sigma] = featureNormalize(X)
%FEATURENORMALIZE Normalizes the features in X
%   FEATURENORMALIZE(X) returns a normalized version of X where
%   the mean value of each feature is 0 and the standard deviation
%   is 1. This is often a good preprocessing step to do when
%   working with learning algorithms.

% You need to set these values correctly
X_norm = X;
mu = zeros(1, size(X, 2));
sigma = zeros(1, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: First, for each feature dimension, compute the mean
%               of the feature and subtract it from the dataset,
%               storing the mean value in mu. Next, compute the
%               standard deviation of each feature and divide
%               each feature by it's standard deviation, storing
%               the standard deviation in sigma.
%
%               Note that X is a matrix where each column is a
%               feature and each row is an example. You need
%               to perform the normalization separately for
%               each feature.
%
% Hint: You might find the 'mean' and 'std' functions useful.
%

mu = mean(X);
sigma = std(X);
X_norm = (X - repmat(mu, [size(X,1), 1]))./repmat(sigma, [size(X,1), 1]);

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

end

repo: machine-learning-coursera


text:
%% Rearranging CHARGUI to Make it Resizable and Dockable
% This script reverse-engineers |charGUI| and rebuilds it in a way that is
% both resizable and dockable.  It does so using |uipanel| and positioning
% with normalized units.

%% Create a New CHARGUI
% First, create a new instance of |charGUI|.
f = charGUI;

%% Cache the Major GUI Components
% Next, find the major GUI components: |axes|, |uipanel|, and |uicontrol|.
% Keeping references to these items will allow us to lay them out later.
% In the meantime, turn off their visiblity so that the steps taken below
% appear clearly in the GUI.
axs = findobj(f,'type','axes');
pans = findobj(f,'type','uipanel');
uis = findobj(f,'type','uicontrol');

set(axs,'visible','off');
set(pans,'visible','off');
set(uis,'visible','off');

%% Construct a New Layout Framework
% Use |uipanel| to create the framework for the new GUI layout.  Divide the
% space into three portions using normalized units.  The first, occupying
% the upper left-hand portion of the GUI, will contain the main axes.  The
% second, occupying the right-hand quarter of the GUI, will contain the
% every |uicontrol|.  The third, occupuying the bottom 30% of the GUI in
% the left-hand 75%, will contain the processing |axes|.
newpans = zeros(1,3);
for i = 1:3
    newpans(i) = uipanel('parent',f,'units','normalized');
end
set(newpans(1),'position',[0 .3 .75 .7]);
set(newpans(2),'position',[.75 0 .25 1]);
set(newpans(3),'position',[0 0 .75 .3]);

%% Populate the First Panel
% By putting the main set of |axes| into a |uipanel|, it can be laid out
% relative to the panel, as opposed to the whole GUI.  Use normalized
% units to make sure that the |axes| grow when the |uipanel| grows.
% Position it starting 10% up from the bottom and 10% right from the
% left-hand side of the GUI and have it occupy 80% of the horizontal and
% vertical |uipanel| space.
set(axs(1),'parent',newpans(1),'units','normalized',...
    'position',[.1 .1 .8 .8]);
set(get(newpans(1),'children'),'visible','on');

%% Populate the Second Panel
% Parent the |uicontrol| items to the second |uipanel|.  Maintain the
% use of characters units so that each |uicontrol| will not grow and shrink
% along with the |uipanel|.  Instead, each will maintain the same size and
% placement.
set(uis,'parent',newpans(2));
newuipos = get(uis(1),'position');
newuipos(1:2) = [2 1];
set(uis(1),'position',newuipos);
newuipos = get(uis(3),'position');
newuipos(1:2) = [2 14];
set(uis(3),'position',newuipos);
newuipos(1:2) = [2 18];
set(uis(4),'position',newuipos);
newuipos(1:2) = [2 22];
set(uis(5),'position',newuipos);
newuipos(1:2) = [2 26];
set(uis(6),'position',newuipos);
newuipos(1:2) = [2 30];
set(uis(2),'position',newuipos);
newuipos(1:2) = [2 34];
set(uis(7),'position',newuipos);
set(get(newpans(2),'children'),'visible','on');

%% Populate the Third Panel
% Each processing |axes| are already contained within a |uipanel|, which
% merely needs to be positioned within one containing |uipanel|.  Again,
% use normalized units so that these will grow and shrink with the
% containing |uipanel|.
%
% Also apply normalized units and new positions to the |axes| contained in
% the subpanels.  Doing so will fill in the GUI space and allow us, as the
% GUI grows, to see more detail in the processing.
set(pans,'parent',newpans(3),'units','normalized');
set(pans(4),'position',[.01 .025 .23 .95]);
set(pans(3),'position',[.26 .025 .23 .95]);
set(pans(2),'position',[.51 .025 .23 .95]);
set(pans(1),'position',[.76 .025 .23 .95]);
set(axs(2:5),'units','normalized','position',[.05 .05 .9 .9]);
set(get(newpans(3),'children'),'visible','on');

%% Touch Up the GUI and Make it Resizable
% Put the finishing touches on the GUI.  Apply the same background colors
% that exsited before to each new |uipanel| and remove the border.  The
% |figure| can now be made resizable.  Try docking it by setting its
% |windowstyle| property!
for i = 1:3
    set(newpans(i),'backgroundcolor',get(pans(1),'backgroundcolor'),...
        'bordertype','none');
end
set(f,'resize','on');


text:
function polynomial_scale_test ( )

%*****************************************************************************80
%
%% POLYNOMIAL_SCALE_TEST tests POLYNOMIAL_ScALE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 October 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYNOMIAL_SCALE_TEST\n' );
  fprintf ( 1, '  POLYNOMIAL_SCALE scales a polynomial by a multiplier S.\n' );

  m = 3;
  o = 6;
  c = [ 7.0, - 5.0, 9.0, 11.0, 0.0, - 13.0 ];
  e = [ 1, 2, 4, 5, 12, 33 ];

  fprintf ( 1, '\n' );
  title = '  P(X) = ';
  polynomial_print ( m, o, c, e, title );

  s = - 0.5;
  fprintf ( 1, '\n' );
  fprintf ( '  Apply scale factor S = %g\n', s );
  [ o, c, e ] = polynomial_scale ( s, m, o, c, e );

  fprintf ( 1, '\n' );
  title = '  S * P(X) = ';
  polynomial_print ( m, o, c, e, title );

  return
end


text:
function RHS = SliceMultiProd(M, X)
% function RHS = SliceMultiProd(M, X)
%
% PURPOSE: Multiple matrix product. This function performs the same
%          task as MULTIPROD but works on different shaping convention
%          for input/output matrices
%
% INPUT
%   M  : 3D array (m x n x p)
%   X  : 3D array (n x q x p)
% OUTPUT
%   RHS: 3D array (m x q x p)
%
% Compute p matrix products:
%   M(:,:,k) * X(:,:,k) = RHS(:,:,k) for all k=1,2,...,p
%
% NOTE -- Special call with X: common right matrix
%   Input argument X can be contracted as (n x q), and automatically
%   expanded by the function
%       M   is 3D array (m x n x p)
%       X   is 2D array (n x q)
%       RHS is 3D array (m x q x p)
%       M(:,:,k) * X = RHS(:,:,k) for all k=1,2,...,p
%
% See also: MultiProd, SliceMultiSolver
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% History: original 11-August-2009

X = permute(X, [1 3 2]);
RHS = MultiProd(M, X);
if size(M,3)>1
    RHS = permute(RHS, [1 3 2]);
end


text:
function move_vline(handle,DoneFcn)
%MOVE_VLINE implements horizontal movement of line.
%
%  Example:
%    plot(sin(0:0.1:pi))
%    h=vline(1);
%    move_vline(h)
%
%Note: This tools strictly requires MOVEX_TEXT, and isn't much good
%      without VLINE by Brandon Kuczenski, available at MATLAB Central.
%<http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=1039&objectType=file>

% Copyright 2004-2010 The MathWorks, Inc.

% This seems to lock the axes position
set(gcf,'Nextplot','Replace')
set(gcf,'DoubleBuffer','on')

h_ax=get(handle,'parent');
h_fig=get(h_ax,'parent');
setappdata(h_fig,'h_vline',handle)
if nargin<2, DoneFcn=[]; end
setappdata(h_fig,'DoneFcn',DoneFcn)
set(handle,'ButtonDownFcn',@DownFcn)

function DownFcn(hObject,eventdata,varargin)
set(gcf,'WindowButtonMotionFcn',@MoveFcn)
set(gcf,'WindowButtonUpFcn',@UpFcn)

function UpFcn(hObject,eventdata,varargin)
set(gcf,'WindowButtonMotionFcn',[])
DoneFcn=getappdata(hObject,'DoneFcn');
if isstr(DoneFcn)
  eval(DoneFcn)
elseif isa(DoneFcn,'function_handle')
  feval(DoneFcn)
end

function MoveFcn(hObject,eventdata,varargin)
h_vline=getappdata(hObject,'h_vline');
h_ax=get(h_vline,'parent');
cp = get(h_ax,'CurrentPoint');
xpos = cp(1);
x_range=get(h_ax,'xlim');
if xpos<x_range(1), xpos=x_range(1); end
if xpos>x_range(2), xpos=x_range(2); end
XData = get(h_vline,'XData');
XData(:)=xpos;
set(h_vline,'xdata',XData)
%update text
text_obj = findobj('Type','Text','Tag','cbar_text');
movex_text(text_obj,xpos)


text:
% Showing the effect of the POCS partial Fourier reconstruction.
%
% Mark a cell (code beginning with "%%") and
% press <Ctrl> + <Enter> to evaluate the cell content only.

%% simple example (2D, single coil, no phase errors)
%
N = 256;
iter = 30;
pf = 9/16;      % a typical large partial fourier factor
p = single(phantom(N));
kspRef = fftshift(fftshift(  fft(fft(  fftshift(fftshift( p, 1),2), [],1),[],2), 1),2);
kspRef = kspRef + 3 * complex( randn(N), randn(N) );    % add noise
kspPF  = kspRef;
kspPF(:,1:floor((1-pf)*N)) = 0;

imReco = pocs( kspPF, iter, true );

imRef    = ifftshift(ifftshift(  ifft(ifft(  ifftshift(ifftshift( kspRef, 1),2), [],1),[],2), 1),2);
imDirect = ifftshift(ifftshift(  ifft(ifft(  ifftshift(ifftshift(  kspPF, 1),2), [],1),[],2), 1),2);

figure,
imagesc(abs([imRef  imReco  imDirect]), [0 2*mean(abs(imRef(:)))])
axis image
title('reference (full dataset)   |   POCS reco   |    zero-padded reco')
colormap(gray(256))

%% multicoil example with phase error
%
clear  p  kspRef  kspPF  imReco  imRef  imDirect
N = 256;
iter = 30;
pf = 9/16;
% Set up some pseudo coilmaps:
cmap(1,:,:) = repmat(        ( 3./(1:N)  .^0.10) .* exp(1i*(1:N)  *pi/(N/3)) , N, 1 );
cmap(2,:,:) = repmat(        ( 1./(1:N).'.^0.20) .* exp(1i*(1:N).'*pi/(N/4)) , 1, N );
cmap(3,:,:) = repmat( fliplr(( 2./(1:N)  .^0.04) .* exp(1i*(1:N)  *pi/(N/8))), N, 1 );
cmap(4,:,:) = repmat( flipud(( 2./(1:N).'.^0.15) .* exp(1i*(1:N).'*pi/(N/10))),1, N );
cmap = single(cmap);

Nc = size(cmap,1);
p = reshape( single(phantom(N)), 1, N, N );
subs = { 1, 164:181, 187:193 };
p(subs{:}) = 0;
subs = { 1, 166:179, 189:191 };
p(subs{:}) = 1 * exp(1i * pi/2);
p = bsxfun( @times, cmap, p );
kspRef = fftshift(fftshift(  fft(fft(  fftshift(fftshift( p, 2),3), [],2),[],3), 2),3);
kspRef = kspRef + 10 * complex( randn(Nc,N,N), randn(Nc,N,N) );    % add noise
kspPF  = kspRef;
kspPF(:,:,1:floor((1-pf)*N)) = 0;

imReco = pocs( kspPF, iter, true );

imReco   = squeeze(sqrt(sum( abs(imReco).^2, 1)));
imRef    = ifftshift(ifftshift(  ifft(ifft(  ifftshift(ifftshift( kspRef, 2),3), [],2),[],3), 2),3);
imRef    = squeeze(sqrt(sum( abs(imRef).^2, 1)));
imDirect = ifftshift(ifftshift(  ifft(ifft(  ifftshift(ifftshift(  kspPF, 2),3), [],2),[],3), 2),3);
imDirect = squeeze(sqrt(sum( abs(imDirect).^2, 1)));

figure,
imagesc(abs([imRef  imReco  imDirect]),[0 2*mean(abs(imRef(:)))])
axis image
title('reference (full dataset)   |   POCS reco   |    zero-padded reco')
colormap(gray(256))


text:
classdef TethaedraSubMeshConnecComputer < handle

    properties (Access = public)
        nSubCellNodes
        nSubCellsByElem
    end

    properties (Access = private)
        connecCase
    end

    properties (Access = private)
        nElemInCase
        cellNodes
    end

    methods (Access = public)

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

        function compute(obj,nodes,icase)

            switch mode(size(nodes,2))
                case 7
                    obj.nSubCellsByElem = 4;
                case 8
                    obj.nSubCellsByElem = 6;
            end

            obj.nElemInCase = size(nodes,1);
            obj.cellNodes   = nodes;
            nodesC = zeros(obj.nSubCellsByElem,obj.nSubCellNodes,obj.nElemInCase);
            isActive = true(obj.nElemInCase,1);
            switch mode(size(nodes,2))
                case 7
                     [nodesT,nodesP] = obj.computeNodesTAndNodesP(icase);
                    for inode = 1:obj.nSubCellNodes
                        nodesC(1,inode,isActive) = nodesT(inode);
                        nodesC(2,inode,isActive) = nodesP(1,inode);
                        nodesC(3,inode,isActive) = nodesP(2,inode);
                        nodesC(4,inode,isActive) = nodesP(3,inode);
                    end
                case 8
                    nodesCC = obj.computeNodesFourCutNodes(icase);
                    for inode = 1:obj.nSubCellNodes
                        nodesC(1,inode,isActive) = nodesCC(1,inode);
                        nodesC(2,inode,isActive) = nodesCC(2,inode);
                        nodesC(3,inode,isActive) = nodesCC(3,inode);
                        nodesC(4,inode,isActive) = nodesCC(4,inode);
                        nodesC(5,inode,isActive) = nodesCC(5,inode);
                        nodesC(6,inode,isActive) = nodesCC(6,inode);
                    end

            end
            obj.connecCase = nodesC;
        end

        function nodePartition = partition(obj,nodes)
            connec = permute(obj.connecCase,[3 1 2]);
            nodePartition = obj.initNodePartition();
            nNodesSubCell = size(nodes,2);
            for isubcell = 1:obj.nSubCellsByElem
                for inode = 1:obj.nSubCellNodes
                    node = connec(:,isubcell,inode);
                    ind  = obj.computeAbsoluteMatrixIndex(node,nNodesSubCell);
                    nodePartition(inode,isubcell,:) = nodes(ind);
                end
            end
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.nSubCellNodes = 4;
        end

        function nodes = computeNodesFourCutNodes(obj,icase)

            switch icase
                case 3
                    nodes1 = obj.prismaTriangulation([1 5 6],[3 7 8]);
                    nodes2 = obj.prismaTriangulation([2 5 7],[4 6 8]);
                    nodes = [nodes1;nodes2];
                case 2
                    nodes1 = obj.prismaTriangulation([2 7 8],[1 5 6]);
                    nodes2 = obj.prismaTriangulation([5 3 7],[6 4 8]);
                    nodes = [nodes1;nodes2];
                case 1
                    nodes1 = obj.prismaTriangulation([4 7 8],[1 5 6]);
                    nodes2 = obj.prismaTriangulation([5 2 7],[6 3 8]);
                    nodes = [nodes1;nodes2];
            end

        end

        function [nodesT,nodesP] = computeNodesTAndNodesP(obj,icase)
            switch icase
                case 4
                    nodesT = [5 6 7 4];
                    nodesP = obj.prismaTriangulation([5 6 7],[1 2 3]);
                    %nodesP = [1 5 6 7;1 2 7 6;1 2 3 7];
                case 3
                    nodesP = [5 6 7 1 ;2 7 6 1;2 4 7 1 ];
                    nodesP = obj.prismaTriangulation([5 7 6],[1 4 2]);
                    nodesT = [5 7 6 3];

                    %nodesP = [1 5 7 6;1 4 6 7;1 4 2 6];
                    %nodesP  = [5 6 7 2;1 5 7 2;4 1 7 2];
                case 2
                    nodesT = [5 6 7 2];
%                    nodesP = [1 5 6 7;1 4 7 6;6 1 4 3];
                    nodesP = obj.prismaTriangulation([5 6 7],[1 3 4]);

                case 1
                    nodesT = [6 7 1 5];
%                    nodesP = [5 6 7 2;4 2 7 6;2 3 6 4];
                    nodesP = obj.prismaTriangulation([5 7 6],[2 4 3]);

            end

        end

       function nodes = prismaTriangulation2(obj,nodesA,nodesB)
             nodes =  [nodesA(1) nodesA(2) nodesA(3) nodesB(1);
                       nodesA(1) nodesA(2) nodesA(3) nodesB(2);
                       nodesA(1) nodesA(2) nodesA(3) nodesB(3)];
       end

        function nodes = prismaTriangulation(obj,nodesA,nodesB)
             nodes =  [nodesB(1) nodesA(1) nodesA(2) nodesA(3);
                       nodesB(1) nodesB(2) nodesA(3) nodesA(2);
                       nodesB(1) nodesB(2) nodesB(3) nodesA(3)];
        end

        function c = initNodePartition(obj)
            c = zeros(obj.nSubCellNodes,obj.nSubCellsByElem,obj.nElemInCase);
        end

        function ind = computeAbsoluteMatrixIndex(obj,colums,nNodesSubCell)
            nElem = obj.nElemInCase;
            rows = (1:nElem)';
            ind = sub2ind([nElem,nNodesSubCell],rows,colums);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function []=panel2hddisp(N,n,T,Units,endo,hd_estimates,stringdates1,decimaldates1,pref)

% plot the figure
if pref.plot
hd=figure('Tag','BEARresults');
set(hd,'Color',[0.9 0.9 0.9]);
set(hd,'name','historical decomposition');
% initiate the count
count=0;
% loop over units
for ii=1:N
   % loop over endogenous variables
   for jj=1:n
      % loop over shocks
      for kk=1:n+1
      % increment count
      count=count+1;
      % then plot
      subplot(N*n,n+1,count)
      temp=hd_estimates{jj,kk,ii};
      hold on
      Xpatch=[decimaldates1' fliplr(decimaldates1')];
      Ypatch=[temp(1,:) fliplr(temp(3,:))];
      HDpatch=patch(Xpatch,Ypatch,[0.7 0.78 1]);
      set(HDpatch,'facealpha',0.5);
      set(HDpatch,'edgecolor','none');
      plot(decimaldates1,temp(2,:)','Color',[0.4 0.4 1],'LineWidth',2);
      plot([decimaldates1(1,1),decimaldates1(end,1)],[0 0],'k--');
      hold off
      set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
         % top labels
         if (jj==1 && kk<=n)
         title([Units{ii,1} '\_' endo{kk,1}],'FontWeight','normal');
         elseif (jj==1 && kk==n+1)
         title('Exogenous','FontWeight','normal');
         end
         % side labels
         if kk==1
         ylabel([Units{ii,1} '\_' endo{jj,1}],'FontWeight','normal');
         end
      end
   end
end
% top supertitle
ax=axes('Units','Normal','Position',[.11 .075 .85 .88],'Visible','off');
set(get(ax,'Title'),'Visible','on')
title('Contribution of shock:','FontSize',11,'FontName','Times New Roman','FontWeight','normal');
% side supertitle
ylabel('Fluctuation of:','FontSize',12,'FontName','Times New Roman','FontWeight','normal');
set(get(ax,'Ylabel'),'Visible','on')
end

% save on Excel

% create the cell that will be saved on excel
hdcell={};
% build preliminary elements: space between the tables
vertspace=repmat({''},T+3,1);
horzspace=repmat({''},2,5*(n+1));
% loop over units
for ii=1:N
   % loop over endogenous variables
   for jj=1:n
   % initiate the cell of results
   endocell={};
      % loop over shocks
      for kk=1:n
      % create a header
      header=[{['fluctuation of ' Units{ii,1} '_' endo{jj,1} ' due to ' Units{ii,1} '_' endo{kk,1} ' shocks']} {''} {''} {''};{''} {''} {''} {''};{''} {'lower bound'} {'median'} {'upper bound'}];
      % complete the cell
      tempcell=[[header;stringdates1 num2cell((hd_estimates{jj,kk,ii}(1,:))') num2cell((hd_estimates{jj,kk,ii}(2,:))') num2cell((hd_estimates{jj,kk,ii}(3,:))')] vertspace];
      % concatenate to the previous parts of unitcell
      endocell=[endocell tempcell];
      end
   % exogenous contribution
   header=[{['fluctuation of ' Units{ii,1} '_' endo{jj,1} ' due to exogenous']} {''} {''} {''};{''} {''} {''} {''};{''} {'lower bound'} {'median'} {'upper bound'}];
   % complete the cell
   tempcell=[[header;stringdates1 num2cell((hd_estimates{jj,n+1,ii}(1,:))') num2cell((hd_estimates{jj,n+1,ii}(2,:))') num2cell((hd_estimates{jj,n+1,ii}(3,:))')] vertspace];
   % concatenate to the previous parts of unitcell
   endocell=[endocell tempcell];
   % concatenate to the previous parts of irfcell
   hdcell=[hdcell;horzspace;endocell];
   end
hdcell=[hdcell;horzspace];
end
% trim
hdcell=hdcell(3:end-2,1:end-1);
% write in excel
if pref.results==1
    bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),hdcell,'hist decomposition','B2');
end

author: european-central-bank

repo: BEAR-toolbox


text:
function calSave1 (~, ~, A, B)
    % save data to an interactively chosen file

    msg.infodisp('  ','Save Data');

    [file1,path1] = uiputfile(fullfile(ZmapGlobal.Data.Directories.output, '*.dat'), 'Filename ');
    if file1 && path1
        if ~iscolumn(A)
            A = A';
        end
        if ~iscolumn(B)
            B = B';
        end
        data = [A, B]';
        fid = fopen([path1 file1],'w') ;
        fprintf(fid, '%12.5f  %12.5f\n' , data);
        fclose(fid) ;
    else
        msg.dbdisp('  ','cancelled save');
    end

end

author: CelsoReyes


text:
function [weight]=nut_sLORETA(Lp,data,flags)
% weight=nut_sLORETA(Lp,data,flags)
% inputs for regularization contant:
% [1] data.Ryy = sample covariance, for data-dependent regularization
% [2] flags.gamma = user defined regularization constant, or 'auto' for
%     leadfield-based regularization

if nargin<3, flags.gamma=[]; end

L = reshape(Lp,size(Lp,1),size(Lp,2)*size(Lp,3));
G = L*L'; clear L

if isfield(flags,'snr')
    if ~isfield(data,'C') || isempty(data.C)
        data.C = eye(size(data.Ryy));
    end
    gamma = trace(G)/(trace(data.C)*flags.snr^2)
    InvG = inv(G+gamma*data.C);
elseif ~isfield(flags,'gamma') || isempty(flags.gamma)
    gamma = 1e0*max(eig(data.Ryy))
    InvG = inv(G+gamma*eye(size(G)));
elseif isnumeric(flags.gamma)
    gamma = flags.gamma   %* max(eig(data.Ryy))
    InvG = inv(G+gamma*eye(size(G)));
else
    % leadfield-based normalization
    doplot=false;
    x = [-20:20];
    gamma=10.^x; %.* max(eig(data.Ryy));
    numg = length(gamma);
    InvG = zeros([size(G) numg]);
    meanuptr = zeros(numg,1);
    warning('off','MATLAB:nearlySingularMatrix')
    for k=1:numg
        InvG(:,:,k) = inv(G+gamma(k)*eye(size(G)));
        uptr = abs(triu(InvG(:,:,k),1));
        meanuptr(k) = mean(uptr( find(uptr(:)) ));  % mean of matrix without diagonal
    end
    warning('on','MATLAB:nearlySingularMatrix')
    if doplot, figure; plot(x(1:end-2),-diff(diff(log10(meanuptr)))); end
    [dum,idx]=findpeaks(-diff(diff(log10(meanuptr))),'MINPEAKHEIGHT',.1);  % max of second derivation of InvG across gamma magnitudes
    idx = max(idx);
    fprintf('Optimal gamma: 1e%d\n',x(idx));
    InvG = InvG(:,:,idx);
end

Lp1 = squeeze(Lp(:,1,:));
w1 = zeros(size(Lp1));
if size(Lp,2)>1
    Lp2 = squeeze(Lp(:,2,:));
    w2 = zeros(size(Lp2));
end
if size(Lp,2)>2
    Lp3 = squeeze(Lp(:,3,:));
    w3 = zeros(size(Lp3));
end

for i=1:size(Lp,3)
    InvGLp = InvG*Lp1(:,i);
    J = inv(sqrt(Lp1(:,i)'*InvGLp));
    w1(:,i) = InvGLp * J;

    if size(Lp,2)>1
        InvGLp = InvG*Lp2(:,i);
        J = inv(sqrt(Lp2(:,i)'*InvGLp));
        w2(:,i) = InvGLp * J;
    end

    if size(Lp,2)>2
        InvGLp = InvG*Lp3(:,i);
        J = inv(sqrt(Lp3(:,i)'*InvGLp));
        w3(:,i) = InvGLp * J;
    end
end

weight(:,1,:) = w1;
if size(Lp,2)>1
    weight(:,2,:) = w2;
end
if size(Lp,2)>2
    weight(:,3,:) = w3;
end
disp('done');
% end

author: spm


text:
function field = SafeGetField(str, fieldname, default, warn)
%SAFEGETFIELD A safe wrapper to GETFIELD.
%
% FIELD = SAFEGETFIELD(STR, FIELDNAME) is a wrapper to GETFIELD that
% returns an empty matrix (with a warning) rather than an error if the
% field is not found.
%
% FIELD = SAFEGETFIELD(STR, FIELDNAME, DEFAULT) works as above but
% returns the value DEFAULT instead of [] when the field is not found.
%
% FIELD = SAFEGETFIELD(STR, FIELDNAME, DEFAULT, 0) works as above, but
% does not give a warning.
%
% Examples:
%
% x.a = 1;
% SafeGetField(x, 'b')
%
% See also: getfield

% $Author: rcotton $	$Date: 2010/10/01 13:58:35 $	$Revision: 1.1 $
% Copyright: Health and Safety Laboratory 2010

SetDefaultValue(3, 'default', []);
SetDefaultValue(4, 'warn', true);

if hasfield(str, fieldname)
   field = getfield(str, fieldname);    %#ok<GFLD>
else
   if warn
      warning('SafeGetField:NonExistentField', ...
         ['The field ''' fieldname ''' does not exist.']);
   end
   field = default;
end

end

function yn = hasfield(str, fieldname)
   yn = any(strcmp(fieldname, fieldnames(str)));
end


text:
function P = getFamily(L,families)
[K,n] = size(L);
nf = length(families);
C = zeros(K,nf);
for i=1:nf
    indf = families{i};
    C(indf,i) = 1;
end
f0 = C*sum(C,1)'.^-1/size(C,2);
P = zeros(nf,n);
for i=1:n
    logPm = L(:,i) + log(f0);
    Pm = exp(logPm-max(logPm));
    Pm = Pm./sum(Pm);
    P(:,i) = C'*Pm;
end

author: MBB-team

repo: VBA-toolbox


text:
% same as demo.m and demo_1p.m but using the MotionCorrection object. Scroll at the end
% to the 1p example

clear
gcp;

name = 'granule_love2.tif';         % two photon dataset
if ~exist(name,'file')  % download file if it doesn't exist in the directory
    url = 'https://www.dropbox.com/s/mjmtwn4pdgydkny/granule_love2.tif.zip?dl=1';
    filename = 'granule_love2.tif.zip';
    fprintf('downloading the file...');
    outfilename = websave(filename,url);
    fprintf('...done. Now unzipping...')
    unzip(filename);
    fprintf('done.');
end

%% rigid motion correction

MC_rigid = MotionCorrection(name);
options_rigid = NoRMCorreSetParms('d1',MC_rigid.dims(1),'d2',MC_rigid.dims(2),'bin_width',50,'max_shift',15,'us_fac',50,'init_batch',200);
MC_rigid.motionCorrectSerial(options_rigid);  % can also try parallel

%% pw-rigid motion correction (in parallel)

MC_nonrigid = MotionCorrection(name);
options_nonrigid = NoRMCorreSetParms('d1',MC_nonrigid.dims(1),'d2',MC_nonrigid.dims(2),'grid_size',[32,32],'mot_uf',4,'bin_width',50,'max_shift',15,'max_dev',3,'us_fac',50,'init_batch',200);
MC_nonrigid.motionCorrectSerial(options_nonrigid);

%% compute metrics

MC_rigid.correlationMean(10);
MC_rigid.crispness(10);
MC_nonrigid.correlationMean(10);
MC_nonrigid.crispness(10);

%% do some plotting
nnY = quantile(MC_rigid.meanY(:),0.005);
mmY = quantile(MC_rigid.meanY(:),0.995);
T = MC_rigid.T;
cY = MC_rigid.corrY;
cM1 = MC_rigid.corrM;
cM2 = MC_nonrigid.corrM;
figure;
    ax1 = subplot(2,3,1); imagesc(MC_rigid.meanY,[nnY,mmY]);  axis equal; axis tight; axis off; title('mean raw data','fontsize',14,'fontweight','bold')
    ax2 = subplot(2,3,2); imagesc(MC_rigid.meanM,[nnY,mmY]);  axis equal; axis tight; axis off; title('mean rigid corrected','fontsize',14,'fontweight','bold')
    ax3 = subplot(2,3,3); imagesc(MC_nonrigid.meanM,[nnY,mmY]); axis equal; axis tight; axis off; title('mean non-rigid corrected','fontsize',14,'fontweight','bold')
    subplot(2,3,4); plot(1:T,cY,1:T,cM1,1:T,cM2); legend('raw data','rigid','non-rigid'); title('correlation coefficients','fontsize',14,'fontweight','bold')
    subplot(2,3,5); scatter(cY,cM1); hold on; plot([0.9*min(cY),1.05*max(cM1)],[0.9*min(cY),1.05*max(cM1)],'--r'); axis square;
        xlabel('raw data','fontsize',14,'fontweight','bold'); ylabel('rigid corrected','fontsize',14,'fontweight','bold');
    subplot(2,3,6); scatter(cM1,cM2); hold on; plot([0.9*min(cY),1.05*max(cM1)],[0.9*min(cY),1.05*max(cM1)],'--r'); axis square;
        xlabel('rigid corrected','fontsize',14,'fontweight','bold'); ylabel('non-rigid corrected','fontsize',14,'fontweight','bold');
    linkaxes([ax1,ax2,ax3],'xy')

%% plot shifts

patch_id = 1:size(MC_nonrigid.shifts_x,2);
str = strtrim(cellstr(int2str(patch_id.')));
str = cellfun(@(x) ['patch # ',x],str,'un',0);

figure;
    ax1 = subplot(311); plot(1:T,cY,1:T,cM1,1:T,cM2); legend('raw data','rigid','non-rigid'); title('correlation coefficients','fontsize',14,'fontweight','bold')
            set(gca,'Xtick',[])
    ax2 = subplot(312); plot(MC_nonrigid.shifts_x); hold on; plot(MC_rigid.shifts_x,'--k','linewidth',2); title('displacements along x','fontsize',14,'fontweight','bold')
            set(gca,'Xtick',[])
    ax3 = subplot(313); plot(MC_nonrigid.shifts_y); hold on; plot(MC_rigid.shifts_y,'--k','linewidth',2); title('displacements along y','fontsize',14,'fontweight','bold')
            xlabel('timestep','fontsize',14,'fontweight','bold')
    linkaxes([ax1,ax2,ax3],'x'); legend(str)

%%
%%%%%% 1p example following demo_1p.m  %%%%%%%%

clear;
gcp;

name = '/Users/epnevmatikakis/Documents/Ca_datasets/Miniscope/msCam13.avi';
MC_1p = MotionCorrection(name);
MC_1p.HPF(7,17);                                        % high pass filter
options_1p = NoRMCorreSetParms('d1',MC_1p.dims(1),'d2',MC_1p.dims(2),'bin_width',50,'max_shift',20,'iter',1,'correct_bidir',false);
MC_1p.motionCorrectParallel(options_1p);                % correct the high pass filtered file
MC_1p.M = MC_1p.applyShifts(MC_1p.file_orig);           % apply shifts to the original file

%% plot mean images
figure;subplot(121); imagesc(mean(MC_1p.file_orig,3)); title('Mean image (raw)'); axis tight; axis equal;
       subplot(122); imagesc(mean(MC_1p.M,3)); title('Mean image (corrected)'); axis tight; axis equal;

author: flatironinstitute

repo: NoRMCorre


text:
function alt = triaxial(rsc)

% geodetic altitude relative
% to a triaxial ellipsoid

% input

%  rsc = geocentric position vector (kilometers)

% output

%  alt = geodetic altitude (kilometers)

% Orbital Mechanics with Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global flat cx cy cz vcon ucon caxis

% semi-axes

aaxis = 6378.138;
baxis = 6367.0;
caxis = aaxis * (1 - flat);

ucon = aaxis * aaxis - caxis * caxis;

vcon = baxis * baxis - caxis * caxis;

cx = (aaxis * caxis * rsc(1)) ^ 2;
cy = (baxis * caxis * rsc(2)) ^ 2;
cz = caxis ^ 2 * rsc(3);

g = sqrt((rsc(1) / aaxis) ^ 2 + (rsc(2) / baxis) ^ 2 + ...
         (rsc(3) / caxis) ^ 2);

zp0 = rsc(3) / g;

% bracketing interval

zp1 = zp0 - 10;

zp2 = zp0 + 10;

rtol = 1.0e-8;

zp = brent('trifunc', zp1, zp2, rtol);

xp = (aaxis ^ 2 * rsc(1) * zp) / (caxis ^ 2 * rsc(3) + ucon * zp);

yp = (baxis ^ 2 * rsc(2) * zp) / (caxis ^ 2 * rsc(3) + vcon * zp);

dx = rsc(1) - xp;
dy = rsc(2) - yp;
dz = rsc(3) - zp;

alt = sqrt(dx ^ 2 + dy ^ 2 + dz ^ 2);


text:
function F = uminus( F )
%-   Unary minus of a SPHEREFUNV
%   -F returns the SPHEREFUNV negated componentwise.
%
%   UMINUS(F) is called by the syntax '-F'.

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

% Empty check:
if (isempty(F) )
    return
end

% Take uminus of each component:
F.components = cellfun(@uminus, F.components, 'UniformOutput', false);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [yn,en,S] = FXLMSadapt(un,dn,S)

% FXLMSadapt        Filtered-x LMS (FXLMS) Algorithm
%
%                   Perform over the entire length of the input sequence.
%                   The history of output, square error and coefficients of FIR
%                   filters are passed out to extenal
% Arguments:
% un                Input signal
% dn                Desired signal
% S                 Adptive filter parameters as defined in LMSinit.m
% yn                History of output signal
% en                History of error signal
%
% by Lee, Gan, and Kuo, 2008
% Subband Adaptive Filtering: Theory and Implementation
% Publisher: John Wiley and Sons, Ltd

M = length(S.coeffs);                     % Length of FIR filter
mu = S.step;                              % Step size of the NLMS algorithm
leak = S.leakage;                         % Leaky factor for the leaky LMS algorithm
AdaptStart = S.AdaptStart;
w = S.coeffs;                             % Coefficients of FIR filter
u = zeros(M,1);
fu = zeros(M,1);

ITER = length(un);                        % Length of input sequence
yn = zeros(1,ITER);                       % Initialize output sequence to zero
en = zeros(1,ITER);                       % Initialize error sequence to zero
%filt_un = zeros(1,ITER);                 % Initialize filtered input to zero
%est_filt_un = zeros(1,ITER);             % Initialize filtered input to zero (est)
filt_un = filter(S.sec_num,S.sec_den,un); % Filtered input
est_filt_un = filter(S.estsec,1,un);      % Filtered input (est)

if isfield(S,'unknownsys')
    b = S.unknownsys;
    norm_b = norm(b);
    eml = zeros(1,ITER);
    ComputeEML = 1;
else
    ComputeEML = 0;
end

for n = 1:ITER
   u = [filt_un(n); u(1:end-1)];          % Input signal vector (through the actual
                                          %   secondary path)
   fu = [est_filt_un(n); fu(1:end-1)];    % Filtered input signal vector (through the
                                          %   estimated secondary path)
   yn(n) = w'*u;                          % Compute output using inner product of w and u
   en(n) = dn(n)-yn(n);                   % Error computation
   if ComputeEML == 1;
        eml(n) = norm(b-w)/norm_b;        % System error norm (normalized)
   end
   if n >= AdaptStart
        w = (1-mu*leak)*w + (mu*en(n))*fu;% LMS algorithm in leaky mode
        S.iter = S.iter + 1;
   end
end

S.coeffs = w;                             % Coefficient values at the final iteration
if ComputeEML == 1;
    S.eml = eml;
end

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
function [x,zo]=overlapadd(f,win,inc)
%OVERLAPADD join overlapping frames together X=(F,WIN,INC)
%
% Usage for frequency-domain processing:
%       S=...;                              % input signal
%       OV=2;                               % overlap factor of 2 (4 is also often used)
%       INC=20;                             % set frame increment in samples
%       NW=INC*OV;                          % DFT window length
%       W=sqrt(hamming(NW,'periodic'));     % omit sqrt if OV=4
%       W=W/sqrt(sum(W(1:INC:NW).^2));      % normalize window
%       F=rfft(enframe(S,W,INC),NW,2);      % do STFT: one row per time frame, +ve frequencies only
%       ... process frames ...
%       X=overlapadd(irfft(F,NW,2),W,INC);  % reconstitute the time waveform (omit "X=" to plot waveform)
%
% Inputs:  F(NR,NW) contains the frames to be added together, one
%                   frame per row.
%          WIN(NW)  contains a window function to multiply each frame.
%                   WIN may be omitted to use a default rectangular window
%                   If processing the input in chunks, WIN should be replaced by
%                   ZI on the second and subsequent calls where ZI is the saved
%                   output state from the previous call.
%          INC      gives the time increment (in samples) between
%                   succesive frames [default = NW].
%
% Outputs: X(N,1) is the output signal. The number of output samples is N=NW+(NR-1)*INC.
%          ZO     Contains the saved state to allow a long signal
%                 to be processed in chunks. In this case X will contain only N=NR*INC
%                 output samples.
%

%	   Copyright (C) Mike Brookes 2009
%      Version: $Id: overlapadd.m 2470 2012-11-02 15:27:24Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nr,nf]=size(f);            % number of frames and frame length
if nargin<2
    win=nf;                 % default increment
end
if isstruct(win)
    w=win.w;
    if ~numel(w) && length(w)~=nf
        error('window length does not match frames size');
    end
    inc=win.inc;
    xx=win.xx;
else
    if nargin<3
        inc=nf;
    end
    if numel(win)==1 && win==fix(win) && nargin<3       % win has been omitted
        inc=win;
        w=[];
    else
        w=win(:).';
        if length(w)~=nf
            error('window length does not match frames size');
        end
        if all(w==1)
            w=[];
        end
    end
    xx=[];      % partial output from previous call is null
end
nb=ceil(nf/inc);        % number of overlap buffers
no=nf+(nr-1)*inc;       % buffer length
z=zeros(no,nb);                      % space for overlapped output speech
if numel(w)
    z(repmat(1:nf,nr,1)+repmat((0:nr-1)'*inc+rem((0:nr-1)',nb)*no,1,nf))=f.*repmat(w,nr,1);
else
    z(repmat(1:nf,nr,1)+repmat((0:nr-1)'*inc+rem((0:nr-1)',nb)*no,1,nf))=f;
end
x=sum(z,2);
if ~isempty(xx)
    x(1:length(xx))=x(1:length(xx))+xx;     % add on leftovers from previous call
end
if nargout>1            % check if we want to preserve the state
    mo=inc*nr;          % completed output samples
    if no<mo
        x(mo,1)=0;
        zo.xx=[];
    else
        zo.xx=x(mo+1:end);
        zo.w=w;
        zo.inc=inc;
        x=x(1:mo);
    end
elseif ~nargout
    if isempty(xx)
        k1=nf-inc;  % dubious samples at start
    else
        k1=0;
    end
    k2=nf-inc;      % dubious samples at end
    plot(1+(0:nr-1)*inc,x(1+(0:nr-1)*inc),'>r',nf+(0:nr-1)*inc,x(nf+(0:nr-1)*inc),'<r', ...
        1:k1+1,x(1:k1+1),':b',k1+1:no-k2,x(k1+1:end-k2),'-b',no-k2:no,x(no-k2:no),':b');
    xlabel('Sample Number');
    title(sprintf('%d frames of %d samples with %.0f%% overlap = %d samples',nr,nf,100*(1-inc/nf),no));
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function H=comp_warpedfreqresponse(wintype,fc,bw,fs,L,freqtoscale,scaletofreq,varargin)
%COMP_WARPEDFREQRESPONSE  Transfer function of warped filter
%   Usage: H=comp_warpedfreqresponse(wintype,fc,bw,fs,L,freqtoscale);
%          H=comp_warpedfreqresponse(wintype,fc,bw,fs,L,freqtoscale,normtype);
%
%   Input parameters:
%      wintype     : Type of window (from firwin)
%      fc          : Centre frequency, in scale units.
%      bw          : Bandwith, in scale units.
%      fs          : Sampling frequency in Hz.
%      L           : Transform length (in samples).
%      freqtoscale : Function to convert Hz into scale units.
%      scaletofreq : Function to convert scale units into Hz.
%      normtype    : Normalization flag to pass to |setnorm|.

definput.import={'setnorm'};
definput.flags.symmetry = {'nonsymmetric','symmetric'};
[flags,kv]=ltfatarghelper({},definput,varargin);

fcwasnegative = fc < 0;

if fcwasnegative && flags.do_symmetric
   fc = -fc;
end

fcscale = freqtoscale(fc);

if ~flags.do_symmetric
   % Compute the values in Aud of the channel frequencies of an FFT of
   % length L.
   bins_lo   = freqtoscale(modcent(fs*(0:L-1)/L,fs)).';
else
   bins_lo   = freqtoscale(fs*(0:L-1)/L).';
end

% This one is necessary to represent the highest frequency filters, which
% overlap into the negative frequencies.
nyquest2  = 2*freqtoscale(fs/2);
bins_hi   = nyquest2+bins_lo;

% firwin makes a window of width 1 centered around 0 on the scale, so we rescale the
% bins in order to pass the correct width to firwin and subtract fc
bins_lo=(bins_lo-fcscale)/bw;
bins_hi=(bins_hi-fcscale)/bw;

pos_lo=comp_warpedfoff(fc,bw,fs,L,freqtoscale,scaletofreq,flags.do_symmetric);
% The "floor" below often cuts away a non-zero sample, but it makes
% the support stay below the limit needed for the painless case. Same
% deal 4 lines below.
pos_hi=floor(scaletofreq(fcscale+.5*bw)/fs*L);

if pos_hi>L/2
    % Filter is high pass and spilling into the negative frequencies
    pos_hi=floor(scaletofreq(fcscale+.5*bw-nyquest2)/fs*L);
end;

win_lo=firwin(wintype,bins_lo);
win_hi=firwin(wintype,bins_hi);

H=win_lo+win_hi;
H(isnan(H)) = 0;

H=setnorm(H,flags.norm);

H=circshift(H,-pos_lo);
upidx=modcent(pos_hi-pos_lo,L);

% ------ Testing ---------------
if 0
    bb=circshift(bins_lo,-pos_lo);
    if bb(1)<-0.5
        % Adjust bin_lo
        error('Could do better here.');
    end;
    if (bb(upidx+1)<0.5) && (bb(upidx+1)>0)
        disp('Chopped non-zero sample.');
        bb(upidx+1)
    end;
end;

H=H(1:upidx);

if fcwasnegative && flags.do_symmetric
   H = H(end:-1:1);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function best_acc = FDDL_top(dataset, N_train, k, lambda1, lambda2)
% function FDDL_top(dataset, N_train, k, lambda1, lambda2)
    %% Dependencies
    addpath('utils');
    addpath('LRSDL_FDDL');
    addpath('ODL');
    %% test mode
    if nargin == 0
        dataset = 'myARgender';
        N_train = 350;
        k = 25;
        dataset = 'myYaleB';
        N_train = 10;
        k = 8;
        k0 = 5;
        lambda1 = 0.001;
        lambda2 = 0.05;
    end
    %% get data
    t = getTimeStr();
    [dataset, Y_train, Y_test, label_train, label_test] = train_test_split(...
        dataset, N_train);
    %% main
    [acc, rt] = FDDL_wrapper(Y_train, label_train, Y_test , label_test, ...
                            k, lambda1, lambda2)
    fprintf('rt = %5.1f\n', rt);
    %% output filename
    if ~exist('results', 'dir')
        mkdir('results');
    end
    if ~exist(fullfile('results', 'FDDL'), 'dir')
        mkdir('results', 'FDDL');
    end
    fn = fullfile('results', 'FDDL', strcat(dataset, ...
        '_N_', num2str(N_train), '_k_', num2str(k),'_l1_', ...
        num2str(lambda1), '_l2_', num2str(lambda2),'_', t, '.mat'));
    disp(fn);
    save(fn, 'acc', 'rt');
    best_acc = max(acc);
end

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function CSPFilter = spatial_filtering( x_flt, oldSMC, numPatterns, verbose )

if verbose == 1
    fprintf( '\tSpatial Filtering...\n\t\t' );
end

CSPFilter = cell( 1, oldSMC.numBands );
for i=1:oldSMC.numBands
    if verbose == 1
        if mod(i, 5) == 0
            fprintf( '%d', i );
        else
            fprintf( '.' );
        end
        if mod(i, 100) == 0
            fprintf( '\n' );
        end
    end

    D1 = x_flt{1, i};
    D2 = x_flt{2, i};
    [W, D] = myTrainCSP( D1, D2 );
    CSPFilter{i}.W = W( :, [1:numPatterns, end-numPatterns+1:end] );
    Dd = diag(D);
    CSPFilter{i}.D = Dd([1:numPatterns, end-numPatterns+1:end]);
end
fprintf( '\n' );

author: PatternRecognition

repo: OpenBMI


text:
function b = eq(G1,G2)
% implements G1 == G2

b = all(G1.points == G2.points);

author: mtex-toolbox

repo: mtex

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


text:
function stopnow = mystopfun(problem, x, info, last)
    stopnow = (last >= 3 && info(last-2).cost - info(last).cost < 1e-3);
end

author: guohuayan


text:
clear;
clc;

root1 = './Train_argu/shadow/';
image_path1 = dir(fullfile(root1,'*.jpg'));
root3 = './Train_argu/shadow_free/';
image_path3 = dir(fullfile(root3,'*.jpg'));

save_root1 = './Train_argu/shadow_lab/';
save_root3 = './Train_argu/shadow_free_lab/';

colorTransform = makecform('srgb2lab');

for i=1: length(image_path1)
    i

    image1 = imread([root1 image_path1(i).name]);
    image3 = imread([root3 image_path3(i).name]);

    image1_lab = applycform(image1, colorTransform);
    image3_lab = applycform(image3, colorTransform);

    imwrite(image1_lab,[save_root1  image_path1(i).name]);
    imwrite(image3_lab,[save_root3  image_path1(i).name]);
end

author: xw-hu

repo: DSC

save_path: github-repos/MATLAB/xw-hu-DSC


text:
function [gg h] = collapsegroup(g)
% [gg h] = collapsegroup(g)
% g: group values. col vec of ints, or col cellstr.
% gg: col vec of ints taking values in the interval [1,Ngroups].
% h: if g is a col int vec, h is an int col vec. if g is a col cellstr, h
% is a col cellstr. h is an Ngroups-by-1 vector.
%
% gg(i) is the new grouping value for g(i). h(gg(i)) is equal to g(i).
%
% This is similar to stats/grp2idx, but does something I prefer with nans.

assert(isnumeric(g)||iscellstr(g));
g = g(:);

h = unique(g);

% nans are not uniqued; if there is a nan, include one in h
if isnumeric(h) && any(isnan(h))
    h(isnan(h),:) = [];
    h(end+1,1) = nan;
end

[~,gg] = ismember(g,h);

% if there are nans in g/h, gg will have 0 values since a nan in g is not
% considered to match up with a nan in h.
if isnumeric(h) && any(isnan(h))
    gg(gg==0) = numel(h);
end

author: kristinbranson

repo: JAABA


text:
%this is an example, demonstrating a Radial Basis Function method of vector
%acoustic tomography, using acoustic time-of-flight tomography as an
%example.
%This involves sending a sonic signal across a measurement area and timing
%how long it takes to travel.  A number of transducers are situated around
%the measurement area to achieve this.  Since the absolute sonic speed is
%effected by temperature and wind speed, so is the time of flight.  It is
%possible to reconstruct the temperature and wind velocity from the
%collection of time-of-flight data.
%For more information, please see ref [1] or http://blog.nutaksas.com
%
%References
%Wiens, Travis "Sensing of Turbulent Flows Using Real-Time Acoustic
%Tomography." X1Xth Biennial Conference of the NewZealand Acoustical
%Society, 2008.

%    Copyright Travis Wiens 2008
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%    If you would like to request that this software be licensed under a less
%    restrictive license (i.e. for commercial closed-source use) please
%    contact Travis at travis.mlfx@nutaksas.com

N_xd=16;%number of transducers
sig_dt=1e-6;%(s) std of error to add to dt
N_r=25;%number of radial basis functions

recalculate=true;%set this to false to skip recalculating
rbf_data_fname='gerris_rbf_data.mat';
if recalculate

    fprintf('Calculating forward problem...')
    gerris_fname='simsave-5.0.gfs';
    %this data was calculated using the Gerris Flow Solver
    %(http:gfs.sourceforge.net), from the Karman vortex street example

    %since Gerris data is not regular, we'll fit a RBF to the temperature
    %and velocity data, then use that to calculate the flight times
    data=load(gerris_fname);%load [x y z u v]

    N_rbf_f=300;%number of rbf centers for forward problem
    X_scale=10;%rescale data
    Sc_G=[3 3]*X_scale;%scale factor for gerris
    OS_G=[1 0];%(m) offset for gerris data

    X_limits=[-2 2;-2 2]*X_scale;%limits for training
    Xc_fu=[X_limits(1,1)+diff(X_limits(1,:))*rand(1,N_rbf_f);X_limits(2,1)+diff(X_limits(2,:))*rand(1,N_rbf_f)]';%model rbf centres
    Xc_fv=[X_limits(1,1)+diff(X_limits(1,:))*rand(1,N_rbf_f);X_limits(2,1)+diff(X_limits(2,:))*rand(1,N_rbf_f)]';
    Xc_fT=[X_limits(1,1)+diff(X_limits(1,:))*rand(1,N_rbf_f);X_limits(2,1)+diff(X_limits(2,:))*rand(1,N_rbf_f)]';

    x_G=(data(:,1)-OS_G(1))*Sc_G(1);%x for Gerris data
    y_G=(data(:,2)-OS_G(2))*Sc_G(2);%x for Gerris data
    U_idx=find((x_G>X_limits(1,1))&(x_G<X_limits(1,2))&(y_G>X_limits(2,1))&(y_G<X_limits(2,2)));
    x_G=x_G(U_idx);
    y_G=y_G(U_idx);

    k_fu=1*ones(N_rbf_f,1)/X_scale^2;%prescalar
    k_fv=1*ones(N_rbf_f,1)/X_scale^2;
    k_fT=4*ones(N_rbf_f,1)/X_scale^2;
    k_fu(1)=0;%allow for bias
    k_fv(1)=0;
    k_fT(1)=0;

    u0=mean(data(U_idx,4));%(m/s) base velocity
    v0=mean(data(U_idx,5));%(m/s) base velocity
    T=293+data(U_idx,6);%(K)  temperature
    T0=mean(T);%(K) base temperature
    [W_fu]=train_rbf([x_G y_G],data(U_idx,4)-u0,Xc_fu,k_fu);
    [W_fv]=train_rbf([x_G y_G],data(U_idx,5)-v0,Xc_fv,k_fv);
    [W_fT]=train_rbf([x_G y_G],T-T0,Xc_fT,k_fT);

    %set up transducers in a square
    N_dim=2;%number of dimensions
    X_xd=zeros(N_xd,N_dim);%XD positions
    phi=linspace(0,2*pi,N_xd+1);%angle of XD
    phi=phi+0.01*randn(size(phi));%move the xducers a bit
    r_xd=1*X_scale;%radius of circle inscribed in square
    X_xd_c=mean(X_limits,2)';%centre of circle
    for i=1:N_xd;
        if abs(tan(phi(i)))<1%intersects with verical lines
            X_xd(i,:)=X_xd_c+r_xd*[sign(cos(phi(i))) tan(phi(i))*sign(cos(phi(i)))];%set XD on a square
        else
            X_xd(i,:)=X_xd_c+r_xd*[ sign(sin(phi(i)))*tan(phi(i)).^-1 sign(sin(phi(i)))];%set XD on a square
        end
    end

    %now calculate pairs of XD to make paths
    X_xd0=[];
    X_xd1=[];
    idx=[];
    count=1;
    for i=1:N_xd
        for j=1:N_xd
            if i~=j%ignore path from XD to itself
                X_xd0(count,:)=X_xd(i,:);
                X_xd1(count,:)=X_xd(j,:);
                idx(count,:)=[i j];
                count=count+1;
            end
        end
    end
    N_paths=count-1;%number of paths

    X_limits_plot=(repmat(X_xd_c,2,1)+[-r_xd -r_xd;+r_xd +r_xd])';%limits to plot
    N_plot1=[25 25];%number of points to plot [x y]
    x_plot1=linspace(X_limits_plot(1,1),X_limits_plot(1,2),N_plot1(1))';
    y_plot1=linspace(X_limits_plot(2,1),X_limits_plot(2,2),N_plot1(2))';
    [x_mesh1 y_mesh1]=meshgrid(x_plot1,y_plot1);
    X_plot1=[reshape(x_mesh1,[],1) reshape(y_mesh1,[],1)];

    N_plot2=[100 100];%%higher resolution
    x_plot2=linspace(X_limits_plot(1,1),X_limits_plot(1,2),N_plot2(1))';
    y_plot2=linspace(X_limits_plot(2,1),X_limits_plot(2,2),N_plot2(2))';
    [x_mesh2 y_mesh2]=meshgrid(x_plot2,y_plot2);
    X_plot2=[reshape(x_mesh2,[],1) reshape(y_mesh2,[],1)];

    [u_rbf]=sim_rbf(Xc_fu,X_plot1,W_fu,k_fu);%simulate rbf fields
    [v_rbf]=sim_rbf(Xc_fv,X_plot1,W_fv,k_fv);
    [T_rbf]=sim_rbf(Xc_fT,X_plot1,W_fT,k_fT);

    [u_rbf2]=sim_rbf(Xc_fu,X_plot2,W_fu,k_fu);
    [v_rbf2]=sim_rbf(Xc_fv,X_plot2,W_fv,k_fv);
    [T_rbf2]=sim_rbf(Xc_fT,X_plot2,W_fT,k_fT);

    %calculate flight times analytically from RBF networks
    c0=340;%(m/s) sonic speed

    u_int=rbfn_integral(Xc_fu,X_xd0,X_xd1,W_fu,k_fu);%line integrals across rbfn
    v_int=rbfn_integral(Xc_fv,X_xd0,X_xd1,W_fv,k_fv);
    T_int=rbfn_integral(Xc_fT,X_xd0,X_xd1,W_fT,k_fT);
    theta=atan2(X_xd1(:,2)-X_xd0(:,2),X_xd1(:,1)-X_xd0(:,1));%(rad) angle of lines
    l=sqrt(sum((X_xd1-X_xd0).^2,2));%(m) path length
    dt0=l./c0.*(1-(u0*cos(theta)+v0*sin(theta))/c0)-...
            1/c0*(T_int/(2*T0)+u_int.*cos(theta)/c0+v_int.*sin(theta)/c0);%(s) time of flight

    save(rbf_data_fname,'dt0','l','theta','c0','T0','u_rbf','v_rbf','T_rbf','u_rbf2','v_rbf2','T_rbf2','X_plot1',...
        'X_plot2','x_plot1','x_plot2','y_plot1','y_plot2','N_plot1','N_plot2','N_paths','N_dim','X_xd','X_xd0',...
        'X_xd1','r_xd','u0','v0','X_xd_c','X_limits','X_scale','X_limits_plot')
    fprintf('done\n')
else
    load(rbf_data_fname)
end

%%%%%%%%%%%%%
%%%%%%%%%%%
%inversion
%%%%%%%%%%
%%%%%%%%%%%%

fprintf('Initializing inversion...')

basisfunction='gaussian';
k_scale=1/X_scale^2;%prescaler for basis function
alpha_scale=1.0e-6;%parameter for regularization alpha

dt=dt0+sig_dt*randn(size(dt0));%add noise to time of flight measurements

Xc_i=[X_limits(1,1)+diff(X_limits(1,:))*rand(1,N_r);X_limits(2,1)+diff(X_limits(2,:))*rand(1,N_r)]';%inversion rbf centres
k_i=k_scale*ones(N_r,1);
Omegapsi=Omegapsi_rbf(Xc_i,X_xd0,X_xd1,k_i,c0,basisfunction);
OmegaT=rbf_integral(Xc_i,X_xd0,X_xd1,k_i,basisfunction )/(c0*2*T0);
Omega=[OmegaT Omegapsi];
alpha_reg=alpha_scale*norm(Omega'*Omega);%regulatization parameter
A_rbf=pinv(Omega'*Omega+alpha_reg*eye(size(Omega,2)))*Omega';%solve!

%calculate matrices to get fields from weights
[zeta_u1 zeta_v1]=zeta_uv(Xc_i,X_plot1,k_i,basisfunction);
[zeta_u2 zeta_v2]=zeta_uv(Xc_i,X_plot2,k_i,basisfunction);
W_T_tmp=zeros(N_r,1);%placeholder
[tmp zeta_T1]=sim_rbf(Xc_i,X_plot1,W_T_tmp,k_i,basisfunction);
[tmp zeta_T2]=sim_rbf(Xc_i,X_plot2,W_T_tmp,k_i,basisfunction);

fprintf('done\n')

fprintf('Solving for weights...')
d_rbf=-dt+l/c0.*(1-(u0*cos(theta)+v0*sin(theta))/c0);%data vector
W=A_rbf*d_rbf;%weights
W_T=W(1:N_r);
W_U=W((N_r+1):end);
fprintf('done\n')

fprintf('Calculating fields...')
%calculate fields
u_hat1=zeta_u1*W_U;
v_hat1=zeta_v1*W_U;
u_hat2=zeta_u2*W_U;
v_hat2=zeta_v2*W_U;
T_hat1=zeta_T1*W_T;
T_hat2=zeta_T2*W_T;
fprintf('done\n')

EU_rms=sqrt(mean((u_rbf2-u_hat2).^2+(v_rbf2-v_hat2).^2));
ET_rms=sqrt(mean((T_rbf2-T_hat2).^2));

fprintf('EU_rms=%f m/s ET_rms=%f K\n',EU_rms,ET_rms)
fprintf('Cond(Omega^T*Omega)=%0.2e, Cond(Omega^T*Omega+alpha*eye)=%0.2e\n',cond(Omega'*Omega),cond(Omega'*Omega+alpha_reg*eye(size(Omega,2))))
fprintf('%d data, %d parameters, ratio=%f\n',size(Omega,1),size(Omega,2),size(Omega,2)/size(Omega,1))

figure(1)
clim=[min([T_rbf2; T_hat2]) max([T_rbf2; T_hat2])];
clf
subplot(2,1,1)
plot(nan,nan)
hold on
imagesc(x_plot2,y_plot2,reshape(T_rbf2,N_plot2(1),[]),clim)
plot([X_xd0(:,1) X_xd1(:,1)]',[X_xd0(:,2) X_xd1(:,2)]','wx:')
hold off
h=colorbar;
ylabel(h,'dT (K)')
axis tight
xlabel('x (m)')
ylabel('y (m)')
title('Target')

subplot(2,1,2)
plot(nan,nan)
hold on
imagesc(x_plot2,y_plot2,reshape(T_hat2,N_plot2(1),[]),clim)
plot([X_xd0(:,1) X_xd1(:,1)]',[X_xd0(:,2) X_xd1(:,2)]','wx:')
hold off
h=colorbar;
ylabel(h,'dT_{hat} (K)')
axis tight
xlabel('x (m)')
ylabel('y (m)')
title('Reconstruction')

figure(2)
plot(nan,nan)
hold on
imagesc(x_plot2,y_plot2,reshape(sqrt((T_hat2-T_rbf2).^2),N_plot2(1),[]))
%plot([X_xd0(:,1) X_xd1(:,1)]',[X_xd0(:,2) X_xd1(:,2)]','kx-')
hold off
h=colorbar;
ylabel(h,'RMSE (K)')
axis tight
xlabel('x (m)')
ylabel('y (m)')

figure(3)
plot(nan,nan)
hold on
imagesc(x_plot2,y_plot2,reshape(sqrt((u_hat2-u_rbf2).^2+(v_hat2-v_rbf2).^2),N_plot2(1),[]))
h=colorbar;
h1=quiver(x_plot1,y_plot1,reshape(u_rbf,N_plot1(1),[]),reshape(v_rbf,N_plot1(1),[]),'k');
h2=quiver(x_plot1,y_plot1,reshape(u_hat1,N_plot1(1),[]),reshape(v_hat1,N_plot1(1),[]));
hold off
axis tight
xlabel('x (m)')
ylabel('y (m)')
ylabel(h,'E (m/s)')
legend([h1;h2],'Target','Reconstruction')

author: Sable

repo: mcbench-benchmarks


text:
function [v] = var(x, normalizeflag, dim, flag)

% [V] = VAR(X, NORMALIZEFLAG, DIM, FLAG) computes the variance,
% across all cells in x along the dimension dim. Normalizeflag = 1 normalizes
% by N, normalizeflag = [] or 0 normalizes by N-1.
%
% X should be an linear cell-array of matrices for which the size in at
% least one of the dimensions should be the same for all cells. If flag==1, the mean will
% be subtracted first (default behavior, but to save time on already demeaned data, it
% can be set to 0).

if nargin<2 || isempty(normalizeflag)
  normalizeflag = 0;
end

D = checkinput(x);
if nargin<3 || isempty(dim),
  dim = find(D,1,'last');
else
  % check whether the requested dim can be used
  if ~D(dim), error('data can not be concatenated across dimension %d',dim); end
end

if nargin<4,
  flag = 1;
end

if flag,
  m    = mean(x, dim);
  x    = cellvecadd(x, -m);
end

nx   = max(size(x));
nsmp = size2(x, dim, 'cell');
ssmp = cellfun(@sumsq,   x, repmat({dim},1,nx), 'UniformOutput', 0);

if normalizeflag
  N = sum(nsmp);
else
  N = sum(nsmp)-1;
end

v = sum(cell2mat(ssmp), dim)./N;

function [s] = sumsq(x, dim)

s = sum(x.^2, dim);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function test_poincareball()

    k = 3;
    n = 10;
    manifold = poincareballfactory(k, n);

    problem.M = manifold;

    P = manifold.rand();
    problem.cost = @(X) .5*norm(X-P, 'fro').^2;
    problem.egrad = @(X) X-P;
    problem.ehess = @(X, U) U;

    % Check numerically whether gradient and Ressian are correct
    checkgradient(problem);
    drawnow;
    pause;
    checkhessian(problem);
    drawnow;
    pause;

    % Initialization
    X0 = [];

    % Options (not mandatory)
    options.maxiter = inf;
    options.maxinner = 30;
    options.maxtime = 120;
    options.tolgradnorm = 1e-10;

    % Pick an algorithm to solve the problem
    Xopt = trustregions(problem, X0, options);
    % [Xopt, costopt, info] = steepestdescent(problem, X0, options);
    % [Xopt, costopt, info] = conjugategradient(problem, X0, options);

    % Curiously enough, the optimizer sometimes doesn't find the solution
    % P. When this happens, the erroneous columns of Xopt tend to have norm
    % very close to 1, which is also where numerics become tricky.
    Xopt - P
    sum(Xopt.^2)
    sum(P.^2)

    evs = hessianspectrum(problem, Xopt);
    evs = real(evs);
    stairs(sort(evs));
    title(['Eigenvalues of the Hessian of the cost function ' ...
           'at the solution']);
    fprintf('Hessian condition number at solution: %g\n', ...
            max(abs(evs))/min(abs(evs)));
end

author: NicolasBoumal

repo: manopt


text:
function save_gc_hdc_trajectory(gcFileName, hdcFileName, gcData, hdcData)
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

    fileID = fopen(gcFileName,'w');
    len = length(gcData(:,1));
    for ind = 1 : len
        fprintf(fileID,'%d %f %f %f \r\n',ind,gcData(ind,:, :));
    end
    fclose(fileID);

    fileID = fopen(hdcFileName,'w');
    len = length(hdcData(:,1));
    for ind = 1 : len
        fprintf(fileID,'%d %f %f \r\n',ind,hdcData(ind, :));
    end
    fclose(fileID);

end

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function bundle  =  readBDL(filename)
% READBDL Read mesh(es) from a .bdl bundle file
%
% bundle  =  readBDL(filename)
%
% Input:
%   filename  path to .bdl file
% Outputs:
%   bundle  struct containing mesh data

  function mesh = readBDLMeshChunk(fp)
  % read a mesh chunk from a file stream
  %
  % Usage:
  %    mesh = readMeshChunk( fp)
  %
  % mesh fields:
  % V  (required)  vertices
  % F              faces
  % col            colors
  % baseColor      a single color for the whole mesh
  % UV             uv coords

      magicID   = fread(fp, 1, 'uint32');
          chunkSize = fread(fp, 1,  'int64');
          nV        = fread(fp, 1, 'uint32');
          nF        = fread(fp, 1, 'uint32');
          nVperFace = fread(fp, 1, 'uint32');
          nCperPos  = fread(fp, 1, 'uint32');
          nCperUV   = fread(fp, 1, 'uint32');
          nCperCol  = fread(fp, 1, 'uint32');
          nCperNor  = fread(fp, 1, 'uint32');
          strip     = fread(fp, 1, 'uchar');
          mesh.baseColor = fread(fp, 4, 'float32');

          mesh.V   = fread(fp, [nCperPos,nV],'float32');
      mesh.V = mesh.V';
          mesh.UV  = fread(fp, [nCperUV,nV], 'float32');
      mesh.UV = mesh.V';
          mesh.col = fread(fp, [nCperCol,nV],'float32');
      mesh.col = mesh.col';
          if nCperNor ~= 0 && nCperNor ~= 3
            error('nCperNor has to be 0 or 3')
          end
          if nCperNor == 3
            mesh.nor = fread(fp, [nCperNor,nV],'float32');
        mesh.nor = mesh.nor';
          end
          mesh.F = fread(fp,  [nVperFace,nF],'uint32');
          mesh.F = mesh.F+1;
  end

  fp = fopen(filename,'r');
  if fp == -1
     error('could not open file')
  end

  chunk = 1;
  while ~feof(fp)
     magicID = fread(fp,1,'uint32');
     if feof(fp)
        break
     end
     bck =  -sizeof('uint32');
     fseek(fp, bck, 'cof');
     if magicID == 1
    bundle.meshes{chunk} = readBDLMeshChunk(fp);
    chunk = chunk+1;
     else
        chunkSize = fread(fp, 1, 'uint64');
        fseek(fp, chunkSize, 'cof');
     end
  end

author: alecjacobson

repo: gptoolbox


text:
function fx1 = p00_fx1 ( prob, x )

%*****************************************************************************80
%
%% P00_FX1: first derivative of a function specified by problem number.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the number of the problem.
%
%    Input, real X, the point at which F is to be evaluated.
%
%    Output, real FX1, the first derivative of the function at X.
%
  if ( prob == 1 )
    fx1 = p01_fx1 ( x );
  elseif ( prob == 2 )
    fx1 = p02_fx1 ( x );
  elseif ( prob == 3 )
    fx1 = p03_fx1 ( x );
  elseif ( prob == 4 )
    fx1 = p04_fx1 ( x );
  elseif ( prob == 5 )
    fx1 = p05_fx1 ( x );
  elseif ( prob == 6 )
    fx1 = p06_fx1 ( x );
  elseif ( prob == 7 )
    fx1 = p07_fx1 ( x );
  elseif ( prob == 8 )
    fx1 = p08_fx1 ( x );
  elseif ( prob == 9 )
    fx1 = p09_fx1 ( x );
  elseif ( prob == 10 )
    fx1 = p10_fx1 ( x );
  elseif ( prob == 11 )
    fx1 = p11_fx1 ( x );
  elseif ( prob == 12 )
    fx1 = p12_fx1 ( x );
  elseif ( prob == 13 )
    fx1 = p13_fx1 ( x );
  elseif ( prob == 14 )
    fx1 = p14_fx1 ( x );
  elseif ( prob == 15 )
    fx1 = p15_fx1 ( x );
  elseif ( prob == 16 )
    fx1 = p16_fx1 ( x );
  elseif ( prob == 17 )
    fx1 = p17_fx1 ( x );
  elseif ( prob == 18 )
    fx1 = p18_fx1 ( x );
  elseif ( prob == 19 )
    fx1 = p19_fx1 ( x );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_FX1 - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem number = %d\n', prob );
    error ( 'P00_FX1 - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [E,J]=SynthMeasWatsonSHCylSingleRadTortIsoVIsoDot_GPD_B0(x, protocol, fibredir, roots)
% Substrate: Impermeable cylinders with one radius in a homogeneous background.
% Orientation distribution: Watson's distribution with SH approximation
% Signal approximation: Gaussian phase distribution.
% Notes: This version estimates the hindered diffusivity from the free diffusivity
% and packing density using Szafer et al's tortuosity model for randomly
% packed cylinders.
% This version includes an isotropic diffusion compartment with its own
% diffusivity.
% This version includes a stationary water compartment.
% Includes a free parameter for the measurement at b=0.
%
% [E,J]=SynthMeasWatsonSHCylSingleRadTortIsoVIsoDot_GPD_B0(x, protocol, fibredir, roots)
% returns the measurements E according to the model and the Jacobian J of the
% measurements with respect to the parameters.  The Jacobian does not
% include derivates with respect to the fibre direction.
%
% x is the list of model parameters in SI units:
% x(1) is the volume fraction of the intracellular space.
% x(2) is the free diffusivity of the material inside and outside the cylinders.
% x(3) is the radius of the cylinders.
% x(4) is the concentration parameter of the Watson's distribution.
% x(5) is the volume fraction of the isotropic compartment.
% x(6) is the diffusivity of the isotropic compartment.
% x(7) is the volume fraction of the isotropic restriction.
% x(8) is the measurement at b=0.
%
% protocol is the object containing the acquisition protocol.
%
% fibredir is a unit vector along the symmetry axis of the Watson's
% distribution.  It must be in Cartesian coordinates [x y z]' with size [3 1].
%
% roots contains solutions to the Bessel function equation from function
% BesselJ_RootsCyl.
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

S0 = x(8);

% Call the other function to get normalized measurements.
if(nargout == 1)
    Enorm=SynthMeasWatsonSHCylSingleRadTortIsoVIsoDot_GPD(x, protocol, fibredir, roots);
else
   [Enorm,Jnorm]=SynthMeasWatsonSHCylSingleRadTortIsoVIsoDot_GPD(x, protocol, fibredir, roots);
end

E = Enorm*S0;

if(nargout>1)
    J = Jnorm*S0;
    J(:,8) = Enorm;
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function ccl_sparse_test ( )

%*****************************************************************************80
%
%% CCL_SPARSE_TEST uses CCL to build a sparse grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Local parameters:
%
%    Local, integer D, the spatial dimension.
%
%    Local, integer MAXK, the maximum level to check.
%
  d = 10;
  maxk = 7;
  func = 'prod( exp(-(x/2).^2/2)/2/sqrt(2*pi), 2)';
  trueval = fu_integral ( d );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CCL_SPARSE_TEST:\n' );
  fprintf ( 1, '  CCL sparse grid:\n' );
  fprintf ( 1, '  Clenshaw-Curtis Linear sparse grid.\n' );
  fprintf ( 1, '  Exact integral is %g\n', trueval );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   D  Level   Nodes    SG error    MC error\n' );
  fprintf ( 1, '\n' );

  for k = 1 : maxk
%
%  Compute the sparse grid estimate.
%
    [ x, w ] = nwspgr ( 'ccl', d, k );
    fx = eval ( func );
    SGappr = w' * fx;
    SGerror = sqrt ( ( SGappr - trueval ).^2 ) / trueval;
%
%  Average 1000 Monte Carlo estimates.
%
    numnodes = length ( w );
    sim = zeros ( 1000, 1 );
    for r = 1 : 1000
      x = rand ( numnodes, d );
      fx = eval ( func );
      sim(r) = mean ( fx );
    end
    simerror = sqrt ( mean ( ( sim - trueval ).^2 ) ) / trueval;

    fprintf( '  %2d     %2d  %6d  %10.5g  %10.5g\n', ...
      d, k, numnodes, SGerror, simerror )

  end

  return
end