text:
function [x, infos] = sagmu_nmf(V, rank, in_options)
% Stochastic averaging gradient multiplicative update for non-negative matrix factorization (SAGMU-NMF) algorithm.
%
% Inputs:
%       matrix      V
%       rank        rank
%       options     options
% Output:
%       w           solution of w
%       infos       information
%
%
% Reference:
%       H. Kasai,
%       "Accelerated stochastic multiplicative update with gradient averaging for nonnegative matrix factorizations,"
%       EUSIPCO, 2018.
%
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai on March 22, 2017
%
%       Feb. 26, 2018 (Hiroyuki Kasai): Fixed algorithm.
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    m = size(V, 1);
    n = size(V, 2);

    % set local options
    local_options.fast_calc             = true;
    local_options.permute_on            = true;
    local_options.sub_mode              = 'STD';
    local_options.accel                 = false;
    local_options.ls                    = false;
    local_options.h_repeat              = 1;
    local_options.rep_mode              = 'fix';
    local_options.stepsize_ratio        = 1; % stepsize ratio
    local_options.robust                = false;
    local_options.lambda                = 1;

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

    number_of_batches = floor(n/options.batch_size);

    if ~isfield(options, 'accel')
        options.accel = false;
        if options.ls
            options.sub_mode = 'LS';
        else
            options.sub_mode = 'STD';
        end
    else
        if options.accel
            options.sub_mode = 'ACC';
        else
            if options.ls
                options.sub_mode = 'LS';
            else
                options.sub_mode = 'STD';
            end
        end
    end

    if options.accel
        if ~isfield(options, 'h_repeat')
            options.h_repeat = 1;
        else
        end

        if ~isfield(options, 'rep_mode')
            options.rep_mode = 'fix';
        else
        end
    else
        options.h_repeat = 1;
    end

    if options.h_repeat == 1
        if options.ls
            options.sub_mode = 'LS';
        else
            options.sub_mode = 'STD';
        end
    end

    if strcmp(options.rep_mode, 'adaptive')
        rhoh = 1+(m+m*rank)/(1*(rank+1));
        alpha = 2;
        delta = 0.01;
    end

    if ~isfield(options, 'x_init')
        Wt  = rand(m, rank);
        H   = rand(rank, n);
        R   = rand(m, n);
    else
        Wt  = options.x_init.W;
        H   = options.x_init.H;
        R   = options.x_init.R;
    end

    % permute samples
    if options.permute_on
        perm_idx = randperm(n);
    else
        perm_idx = 1:n;
    end
    V = V(:, perm_idx);
    H = H(:, perm_idx);
    R = R(:, perm_idx);

    if options.robust
        mode = 'R-SAGMU-NMF';
    else
        mode = 'SAGMU-NMF';
        R = zeros(m, n);
    end

    % initialize
    method_name = sprintf('%s (%s)', mode, options.sub_mode);
    epoch = 0;
    l = zeros(m, options.batch_size) + options.lambda;
    grad_calc_count = 0;

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

    % prepare arrays for vht and Whht
    vht = cell(number_of_batches,1);
    Whht = cell(number_of_batches,1);

    % store vht and Whht
    cnt = 0;
    for t=1: options.batch_size : n - 1
        cnt = cnt + 1;
        vt = V(:,t:t+options.batch_size-1);
        ht = H(:,t:t+options.batch_size-1);
        rt = R(:,t:t+options.batch_size-1);
        vht{cnt} = vt * ht';
        Whht{cnt} = (Wt * ht + rt) * ht';
    end

    % prepare Delta_minus and Delta_plus
    if options.fast_calc
        Delta_minus = zeros(m, rank);
        Delta_plus = zeros(m, rank);
    end

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, Wt, H, R, options, [], epoch, grad_calc_count, 0);

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

    % set start time
    start_time = tic();

    % main outer loop
    while true

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

        cnt = 0;
        % main inner loop
        for t = 1 : options.batch_size : n - 1
            cnt = cnt + 1;

            % retrieve vt and ht
            vt = V(:,t:t+options.batch_size-1);
            ht = H(:,t:t+options.batch_size-1);

            if ~options.robust

                % uddate ht
                Wtv = Wt' * vt;
                WtW = Wt' * Wt;
                if strcmp(options.sub_mode, 'ACC')
                    if strcmp(options.rep_mode, 'adaptive')
                        gamma = 1;
                        eps0 = 1;
                        j = 1;
                        rhoh_alpha = rhoh*alpha;
                        %while j <= floor(1+rhoh*alpha) &&  gamma >= delta*eps0
                        ht0 = ht;
                        while j <= rhoh_alpha && gamma >= delta*eps0
                            ht = ht .* (Wtv) ./ (WtW * ht);
                            ht = ht + (ht<eps) .* eps;
                            if j == 1
                                eps0 = norm(ht0-ht);
                            end
                            gamma = norm(ht0-ht);
                            j = j+1;
                        end
                    else
                        for iii=1:options.h_repeat
                            ht = ht .* (Wtv) ./ (WtW * ht);
                            ht = ht + (ht<eps) .* eps;
                        end
                        %ht = ht + (ht<eps) .* eps;
                    end
                elseif strcmp(options.sub_mode, 'LS')
                    ht = calc_nls_nmf(vt, Wt, 1e-16);
                    ht = ht + (ht<eps) * eps;
                else
                    ht = ht .* (Wtv) ./ (WtW * ht);
                    ht = ht + (ht<eps) .* eps;
                end
            else

                rt = R(:, t:t+options.batch_size-1);

                % uddate ht
                Wtv = Wt' * vt;
                %WtW = Wt' * Wt;
                if strcmp(options.sub_mode, 'ACC')
                    if strcmp(options.rep_mode, 'adaptive')
                        gamma = 1;
                        eps0 = 1;
                        j = 1;
                        rhoh_alpha = rhoh*alpha;
                        %while j <= floor(1+rhoh*alpha) &&  gamma >= delta*eps0
                        ht0 = ht;
                        while j <= rhoh_alpha && gamma >= delta*eps0
    %                         ht = ht .* (Wtv) ./ (Wt' * (Wt * ht + rt));
    %                         ht = ht + (ht<eps) .* eps;
    %                         rt = rt .* vt ./ (Wt * ht + rt + l);
                            Wh_r = Wt * ht + rt;
                            ht = ht .* (Wtv) ./ (Wt' * Wh_r);
                            ht = ht + (ht<eps) .* eps;
                            rt = rt .* vt ./ (Wh_r + l);
                            if j == 1
                                eps0 = norm(ht0-ht);
                            end
                            gamma = norm(ht0-ht);
                            j = j+1;
                        end
                    else
                        for iii=1:options.h_repeat
    %                         ht = ht .* (Wtv) ./ (Wt' * (Wt * ht + rt));
    %                         ht = ht + (ht<eps) .* eps;
    %                         rt = rt .* vt ./ (Wt * ht + rt + l);
                            Wh_r = Wt * ht + rt;
                            ht = ht .* (Wtv) ./ (Wt' * Wh_r);
                            ht = ht + (ht<eps) .* eps;
                            rt = rt .* vt ./ (Wh_r + l);
                        end
                    end
                else
    %                 ht = ht .* Wtv ./ (Wt' * (Wt * ht + rt));
    %                 ht = ht + (ht<eps) .* eps;
    %                 rt = rt .* vt ./ (Wt * ht + rt + l);
                    Wh_r = Wt * ht + rt;
                    ht = ht .* (Wtv) ./ (Wt' * Wh_r);
                    ht = ht + (ht<eps) .* eps;
                    rt = rt .* vt ./ (Wh_r + l);
                end

            end

            % store vht and Whht
            Whht_org = Whht{cnt};
            vht_org = vht{cnt};

            % update
            vht{cnt} = vt * ht';
            if ~options.robust
                Whht{cnt} = Wt * (ht * ht');
            else
                Whht{cnt} = (Wt * ht + rt) * ht';
            end

            if epoch > 0
                if options.fast_calc
                    delta_minus = Delta_minus/n + (vht{cnt} + Whht_org)/options.batch_size;
                    delta_plus = Delta_plus/n + (Whht{cnt} + vht_org)/options.batch_size;

                    % update Delta_minus and Delta_plus
                    Delta_minus = Delta_minus + (vht{cnt} - vht_org);
                    Delta_plus = Delta_plus + (Whht{cnt} - Whht_org);
                else
                    delta_minus = zeros(m, rank);
                    delta_plus = zeros(m, rank);

                    for jj=1:number_of_batches
                        delta_minus = delta_minus + vht{jj};
                        delta_plus = delta_plus + Whht{jj};
                    end

                    delta_minus = delta_minus/n;
                    delta_plus = delta_plus/n;

                    delta_minus = delta_minus + (vht{cnt} + Whht_org)/options.batch_size;
                    delta_plus = delta_plus + (Whht{cnt} + vht_org)/options.batch_size;
                end

            else
                delta_minus = vht{cnt}/options.batch_size;
                delta_plus = Whht{cnt}/options.batch_size;

                if options.fast_calc
                    % update Delta_minus and Delta_plus
                    Delta_minus = Delta_minus + vht{cnt};
                    Delta_plus = Delta_plus + Whht{cnt};
                end
            end

            % update W
            if options.stepsize_ratio == 1
                Wt = Wt .* (delta_minus ./ delta_plus);
            else
                Wt = (1-ratio)* Wt + ratio * Wt .* (delta_minus ./ delta_plus);
            end

            Wt = Wt + (Wt<eps) .* eps;

            % store new h
            H(:,t:t+options.batch_size-1) = ht;

            % update R
            if options.robust
                R(:,t:t+options.batch_size-1) = rt;
            end

            grad_calc_count = grad_calc_count + m * options.batch_size;
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, Wt, H, R, options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = Wt;
    x.H(:,perm_idx) = H;
    x.R = R;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
% Copyright 2013 The MathWorks, Inc.
function y=playControls()

persistent state;
persistent hpause;
persistent hexit; %#ok<PUSE>
persistent isPaused;

if isempty(state)
    initializeUI;
end

% deal with the state
if state == 2 % exit
    y = 0; % terminate the loop
    state = [];
else
    % block the return until pause is turned off
    while(isPaused)
        drawnow; % keep refreshing the UI
        pause(0.01); % don't hog the CPU
    end

    y = 1;
end

    function exit_callback(~, ~)
        state = 2; % 2 signifies exit
        isPaused = false;
        close gcf;
    end

    function pause_callback(obj, ~)
        % Handle button string toggling
        if get(hpause,'value')
            isPaused = true;
            set(obj, 'string', 'PLAY');
        else
            isPaused = false;
            set(obj, 'string', 'PAUSE');
        end
    end

    function initializeUI
        state = 0;

        scrsz = get(0,'ScreenSize');
        % [left, bottom, width, height]
        %figure('Position',[1 scrsz(4)-50 80 50]);
        h = figure('Position',[scrsz(3)*.25 scrsz(4)*.75 80 50], 'Name', 'Loop controls', ...
            'NumberTitle','off','MenuBar','none');

        hpause = uicontrol('unit','pixel','style','togglebutton','string','PAUSE',...
            'position',[10 10 50 25],'callback',@pause_callback);
        hexit = uicontrol('unit','pixel','style','pushbutton','string','EXIT',...
            'position',[60 10 50 25],'callback',@exit_callback);

        set(h,'HandleVisibility','callback');

        isPaused = false;
    end

end

author: Sable

repo: mcbench-benchmarks


text:
% Digital Video Stabilization and Rolling Shutter Correction using Gyroscopes
% Copyright (C) 2011 Alexandre Karpenko
%
% 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
% 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/>.

function [ff td ts wd err] = solve_projective(gyro, frame_time, W, video_file)

gyro(:,[2 1]) = gyro(:,1:2);    % x & y are swapped
gyro(:,3) = -gyro(:,3);         % rotations about z
dft = diff(frame_time);
dT = diff(gyro(:,4));
theta = ((gyro(1:end-1,1:3) + gyro(2:end,1:3)) / 2) .* dT(:,[1 1 1]);
theta = [0 0 0; cumsum(theta, 1)];

ff = -1572;
td = -0.1564;
ts = 0;

alpha_f = 1e-12;
alpha_td = 1e-12;
alpha_ts = 1e-12;
alpha_wd = 1e-10;
wd = -mean(gyro(:,1:3),1);
%d = [0 0 0];        % gyro drift

% read corresponding movie
xyloObj = mmreader(['data/' video_file '.mov']);
num_frames = xyloObj.NumberOfFrames;
vh = xyloObj.Height;
vw = xyloObj.Width;

frame1 = read(xyloObj, 200);
frame2 = read(xyloObj, 201);

for iter = 1:5000
    for f=100:num_frames-1
        dE_df = 0;
        dE_dtd = 0;
        dE_dts = 0;

        err = 0;
        denom = 0;

        x = W(f).x1;
        y = W(f).x2;
        for i = 1:size(x,1)
            K = [1 0 -vw/2; 0 1 -vh/2; 0 0 ff];

            X = K*[x(i,:)'; 1];
            Y = K*[y(i,:)'; 1];

            tx = frame_time(f) - td - ts * X(2) / vh;
            ty = frame_time(f+1) - td - ts * Y(2) / vh;

            w = interp1(gyro(:,4), gyro(:,1:3), [tx ty], 'linear', 'extrap');
            dw = diff(w);
            dws = (w(2,:) * Y(2) - w(1,:) * X(2)) / vh;
            dth = diff(interp1(gyro(:,4), theta(:,1:3), [tx ty], 'linear', 'extrap'));

            Rx = [1            0           0;
                  0 cos(dth(1)) -sin(dth(1));
                  0 sin(dth(1))  cos(dth(1))];

            Ry = [ cos(dth(2)) 0 sin(dth(2));
                   0           1           0;
                  -sin(dth(2)) 0 cos(dth(2))];

            Rz = [cos(dth(3)) -sin(dth(3)) 0;
                  sin(dth(3))  cos(dth(3)) 0;
                  0           0            1];

            dRx =[0             0           0;
                  0 -sin(dth(1)) -cos(dth(1));
                  0 cos(dth(1))  -sin(dth(1))];

            dRy =[-sin(dth(2)) 0  cos(dth(2));
                   0           0            0;
                  -cos(dth(2)) 0 -sin(dth(2))];

            dRz =[-sin(dth(3)) -cos(dth(3)) 0;
                   cos(dth(3)) -sin(dth(3)) 0;
                   0           0            0];

            R = Rz * Ry * Rx;
            RY = R * Y;

            XXY_XYX = X'*X*RY' - X'*RY*X';
            dE_dts = dE_dts + 2 * (XXY_XYX) * (Rz * Ry * dRx * dws(1) + Rz * dRy * Rx * dws(2) + dRz * Ry * Rx * dws(3)) * Y;
            dE_dtd = dE_dtd + 2 * (XXY_XYX) * (Rz * Ry * dRx * dw(1) + Rz * dRy * Rx * dw(2) + dRz * Ry * Rx * dw(3)) * Y;
            dE_df = dE_df + 2*(RY'*RY)*X(3)*ff + 2*(X'*X)*(RY'*(R*[0;0;f])) - 2*(X'*RY)*(ff*RY(3) + X'*(R*[0;0;f]));
            %{
            if f == 200 && mod(iter-1,20) == 0
                figure(1); clf;
                h1 = imshow(frame2); hold on;
                h2 = imshow(imtransform(frame1, maketform('affine', M'), 'XData', [1 vid_width], 'YData', [1 vid_height]));
                set(h2, 'AlphaData', 0.6);
                hold off;
                pause(0.01);
            end
            %}
            err = err + (X'*X)*(Y'*Y) - (X'*RY)^2;
            denom = denom + 1;
        end

        err = err / denom;

        ff = ff + alpha_f * dE_df;
        %td = td - alpha_td * dE_dtd;
        %ts = ts - alpha_ts * dE_dts;

        fprintf('f: %f, td: %f, ts: %f, d: (%f, %f, %f), error: %f\n', ff, td, ts, 0, 0, 0, err);
    end
end

author: alex-golts

repo: Video-Stabilization


text:
function [c, frec, info] = multidgtrealmp(f,dicts,varargin)
%MULTIDGTREALMP  Matching Pursuit Decomposition with Multi-Gabor Dictionary
%   Usage:  c = multidgtrealmp(f,dicts)
%           c = multidgtrealmp(f,dicts,errdb,maxit)
%           [c,frec,info] = multidgtrealmp(...)
%
%   Input parameters:
%       f        : Input signal
%       dicts    : Dictionaries. Format {g1,a1,M1,g2,a2,M2,...}
%       errdb    : Target normalized approximation error in dB
%       maxit    : Maximum number of iterations.
%   Output parameters:
%       c        : Sparse representation
%       frec     : Reconstructed signal
%       info     : Struct with additional output paramerets
%
%   `multidgtrealmp(f,{g1,a1,M1,g2,a2,M2,...,gW,aW,MW})` returns sparse
%   representation of a signal in `W` Gabor dictionaries using the
%   fast matching pursuit algorithm. `gw` is a Gabor window defined
%   as in |dgt| and |dgtreal|, `aw` is a hop factor, `Mw` is the number of
%   frequency channels. All `aw` and `Mw` must be divisible by `min(a1,...,aW)`.
%   The algorithm will try to reach -40 dB relative approximation error
%   in at most `numel(f)` iterations.
%   The function returns a cell-array with elements storing coefficients
%   for individual Gabor systems such that they can be directly used in
%   |idgtreal|.
%
%   `multidgtrealmp(f,dicts,errdb,maxit)` tries to reach normalized
%   approximation error *errdb* dB in at most *maxit* iterations.
%
%   `[c,frec,info] = multidgtrealmp(...)` in addition returns the
%   aproximation *frec* and a struct `info` with the following fields:
%
%     .iter     Number of iterations done.
%
%     .atoms    Number of atoms selected.
%
%     .relres   Final approximation error. If resets are enabled, this is a
%               vector additionally containing the approximation error at
%               every reset.
%
%     .g        Cell array of numeric windows used in the multi-dictionary
%
%     .a        Array of hop factors for indivitual dictionaries
%
%     .M        Array of numbers of channels for individual dictionaries
%
%     .synthetize  Anonymous function which can be used to synthetize from
%                  the (modified) coefficients as
%                  `frec = sum(info.synthetize(c),dim)`
%                  where `dim=2` if the input *f* was a column vector and
%                  `dim=1` if it was a row vector.
%
%   The normalized approximation error is computed as
%   `errdb=20*log10(norm(f-frec)/norm(f))`.
%
%   The function takes the following optional parameters at the end of
%   the line of input arguments:
%
%     'kenrnthr',kt     Kernel threshold. Must be in range ]0,1].
%                       Default is 1e-4.
%
%     'timeinv'         Use the time invariant phase convention. The
%                       default is 'freqinv'.
%
%     'pedanticsearch'  Be pedantic about the energy of pairs of
%                       conjugated atoms in the selection step.
%                       Disabled by default.
%
%     'algorithm',alg   Algorithm to use. Available:
%                       'mp'(default),'selfprojmp','cyclicmp'
%
%     'reset'           Reset the decomposition until stopping criteria
%                       are met. The reset means a complete synthesis and
%                       re-analysis. Available:
%                       'noreset' (default), 'reset', 'quickreset'
%                       Option 'quickreset' applies a heuristic to stop the
%                       algorithm more quickly when the targer error is
%                       reached.
%
%   Reset conditions can be controlled by the following key-value pairs:
%
%     'resetit',it      Reset the decomposition every `it` iteration.
%                       Must be a positive integer.
%
%     'reseterrdb',err  Reset the decomposition when the error estimate
%                       in dB decreases below `err`. Must be negative.
%                       Default is `10*log10(1.1*kv.kernthr)`.
%
%   The computational routine is only available in C. Use |ltfatmex| to
%   to compile it.
%
%   Algorithms
%   ----------
%
%   By default, the function uses the fast MP using approximate update
%   in the coefficient domain as described in::
%
%   "Z. Prusa, N. Holighaus, P. Balazs: Fast Matching Pursuit with Multi-Gabor Dictionaries"
%
%   The kernel threshold limits the minimum approximation error which
%   can be reached. For example the default threshold 1e-4 in general
%   allows achieving at least -40 dB.
%
%   Note that this algorithm internally computes the error as
%   `errdb = 10*log10(err(k)/norm(f)^2)`, with `err(0) = norm(f)^2` and
%   `err(k) = err(k-1) - abs(cselected(k))^2`. Here, `cselected(k)` is the
%   Gabor coefficient selected at the k-th MP step. This formula is exact
%   only if the kernel threshold equals 0. Otherwise it serves as a cheap
%   estimate. The exact error is only computed when this function is run
%   with resets enabled and only in between resets.
%
%   Examples
%   --------
%
%   The following example shows the decomposition in 3 dictionaries and
%   plots contributions from the individual dictionaries and the residual.:::
%
%       [f,fs] = gspi;
%       [c, frec, info] = multidgtrealmp(f,...
%       {'blackman',128,512,'blackman',512,2048,'blackman',2048,8192});
%       frecd = info.synthetize(c);
%       figure(1);
%       xvals = (0:numel(f)-1)/fs;
%       subplot(4,1,1); plot(xvals,frecd(:,1));ylim([-0.5,0.5]);
%       subplot(4,1,2); plot(xvals,frecd(:,2));ylim([-0.5,0.5]);
%       subplot(4,1,3); plot(xvals,frecd(:,3));ylim([-0.5,0.5]);
%       subplot(4,1,4); plot(xvals,f-frec);ylim([-0.5,0.5]); xlabel('Time (s)');
%
%   See also: dgtreal idgtreal
%
%   References: ltfatnote052 mazh93 stchr10 rerose17

%AUTHOR: Zdenek Prusa

thismfile = upper(mfilename);
complainif_notenoughargs(nargin,2,thismfile);
% Define initial value for flags and key/value pairs.
definput.keyvals.errdb=-40;
definput.keyvals.maxit=[];
definput.keyvals.resetit=0;
definput.keyvals.reseterrdb = 0;
%definput.keyvals.iterstep=[];
definput.keyvals.kernthr = 1e-4;
%definput.flags.print={'quiet','print'};
definput.flags.algversion={'fast','slow'};
definput.flags.algorithm={'mp','selfprojmp','cyclicmp'};
definput.flags.search={'plainsearch','pedanticsearch'};
definput.flags.phaseconv={'freqinv','timeinv'};
definput.flags.reset={'noreset','reset','quickreset'};

[flags,kv]=ltfatarghelper({'errdb','maxit'},definput,varargin);

if exist('comp_multidgtrealmp','file') ~= 3 && flags.do_fast
    error(['%s: MEX/OCT file is missing. Either compile the MEX/OCT ',...
           'interfaces or re-run the function with ''slow'''], thismfile);
end

if flags.do_slow
    error('%s: ''slow'' is not supported yet.',thismfile)
end

if flags.do_reset || flags.do_quickreset % Check reseterrdb or resetit
   if kv.reseterrdb == 0
        if kv.resetit == 0 % Set both to conservative default values depending on kernthr
            kv.reseterrdb = 10*log10(1.1*kv.kernthr);
            kv.resetit = round(1000/sqrt(kv.kernthr));
        else % If only resetit is given, ignore reseterrdb
            kv.reseterrdb = -flintmax;
        end
   elseif kv.resetit == 0 % If only reseterrdb is given, ignore resetit
        kv.resetit = flintmax;
   end
   if kv.reseterrdb >= 0
       error('%s: Reset error tolerance must be negative.',thismfile);
   end
   if kv.resetit < 0 || kv.resetit ~= round(kv.resetit)
       error('%s: Number of iterations before reset must be a positive integer.',thismfile);
   end
end

%% ----- step 1 : Verify f and determine its length -------
% Change f to correct shape.
[f,~,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,[],[],upper(mfilename));

if W>1
    error('%s: Input signal can be single channel only.',upper(mfilename));
end

if kv.errdb > 0
    error('%s: Target error must be lower than 0 dB.',upper(mfilename));
end

if ~(kv.kernthr > 0 && kv.kernthr <= 1)
    error('%s: Kernel threshold must be in range ]0,1].',upper(mfilename));
end

if ~iscell(dicts), error('%s: dicts must be cell',thismfile); end
if rem(numel(dicts),3) ~= 0 || ~all(cellfun(@(x)isscalar(x), dicts([2:3:end,3:3:end])))
    error('%s: bad format of dicts. Check {g1,a1,M1,g2,a2,M2,...,gW,aW,MW}',...
        thismfile);
end

dictno = numel(dicts)/3;
gin = dicts(1:3:end);
a = cell2mat(dicts(2:3:end));
M = cell2mat(dicts(3:3:end));

if any(rem(M,a) ~= 0) || any(M./a<2)
    error(['%s: Only integer oversampling greater than 1 is allowed ',...
           'i.e. M/a must be an integer>=2.'],...
    upper(mfilename));
end

if dictno > 1
    asort  = sort(a);
    Msort  = sort(M);
    if any(rem(asort(2:end),asort(1:end-1)) ~= 0)
        error('%s: all au and av must be divisible by min(au,av)',thismfile);
    end
    if any(rem(Msort(2:end),Msort(1:end-1)) ~= 0)
        error('%s: all Mu and Mv must be divisible by min(Mu,Mv)',thismfile);
    end
end

info.a = a;
info.M = M;
info.iter = 0;
info.relres = [];
fnorm = norm(f);

L = filterbanklength(Ls,[a(:);M(:)]);
if isempty(kv.maxit), kv.maxit = L; end

info.g = cell(dictno,1);
for dIdx = 1:dictno
    info.g{dIdx} = setnorm(gabwin(gin{dIdx},a(dIdx),M(dIdx),L),'2');
end

for dIdx = 1:dictno
    condnum = gabframebounds(info.g{dIdx},a(dIdx),M(dIdx));
    if condnum > 1e3
        error('%s: Dictionary %d is badly conditioned.',dIdx,upper(mfilename));
    end
end

% Initial residuum
fpad = postpad(f,L);

% Determine reset/no reset mode and run MP decomposition
if flags.do_noreset
    [c,info.atoms,info.iter,info.status] = ...
        comp_multidgtrealmp(fpad,info.g,a,M,flags.do_timeinv,...
                            kv.kernthr,kv.errdb,kv.maxit,kv.maxit,...
                            flags.do_pedanticsearch, flags.algorithm );
else
    c = cell(1,dictno);
    for dIdx = 1:dictno
        c{dIdx} = zeros(floor(M(dIdx)/2) + 1, L/a(dIdx),class(f));
    end

    remainingit = kv.maxit;
    currerrdb = kv.errdb;

    info.resets = -1;
    reseterrdb = kv.reseterrdb;

    while remainingit > 0 && currerrdb < 0
        info.resets = info.resets + 1;

        % Adjust number of iterations
        currit = min([remainingit,kv.resetit]);
        if flags.do_quickreset % Adjust reset error condition
            reseterrdb = max([1.2*currerrdb,kv.reseterrdb]);
        end

        [ctmp,info.atoms,info.iter,info.status] = ...
            comp_multidgtrealmp(fpad,info.g,a,M,flags.do_timeinv,...
                                kv.kernthr,reseterrdb,currit,currit,...
                                flags.do_pedanticsearch, flags.algorithm );

        c = cellfun(@(cEl,ctmpEl) cEl + ctmpEl, c(:)',ctmp(:)','UniformOutput', 0);

        remainingit = remainingit - info.iter;

        %The following means that some other stopping criterion has been reached
        if info.status ~=0 && info.status ~=3 && currit ~= info.iter, break; end

        if remainingit > 0
            % Recompute residuum
            frec = sum(cell2mat(cellfun(@(cEl,gEl,aEl,MEl) idgtreal(cEl,gEl,aEl,MEl,flags.phaseconv),...
                   c(:)',info.g(:)',num2cell(a(:))',num2cell(M(:))','UniformOutput',0)),2);

            fpad = postpad(f,L) - frec;
            relres = norm(fpad(:)) / fnorm;
            % Since we are changing the residual, we must adjust the target error
            currerrdb = kv.errdb - 20*log10(relres);
            if currerrdb < 0
                info.relres(info.resets+1) = relres;
            end
            if info.resets >= 1 && info.relres(end) > info.relres(end-1)
                info.status = 6;
                break;
            end
        end
    end

    %Fix the info
    info.atoms = sum(cellfun(@(cEl) numel(find( abs(cEl) > 0 )), c));
    info.iter = kv.maxit - remainingit;
end

if nargout>1
  permutedsize2 = permutedsize; permutedsize2(2) = dictno;
  info.synthetize = @(c) ...
      assert_sigreshape_post(...
      postpad(cell2mat(cellfun(@(cEl,gEl,aEl,MEl) idgtreal(cEl,gEl,aEl,MEl,flags.phaseconv),...
      c(:)',info.g(:)',num2cell(a(:))',num2cell(M(:))','UniformOutput',0)),Ls),...
      dim,permutedsize2,order);
  dim2 = 2;
  if dim == 2, dim2 = 1; end
  frec = sum(info.synthetize(c),dim2);
  if fnorm == 0
    info.relres(end+1) = 0;
  else
    info.relres(end+1) = norm(frec(:)-f(:))/fnorm;
  end

  % Fix the status
  if info.relres(end) <= 10^(kv.errdb/20)
      info.status = 0;
  elseif info.iter == kv.maxit
      info.status = 2;
  end

  status_str = {...
  'Target error reached',...
  'Maximum number of atoms reached',...
  'Maximum number of iterations reached',...
  'Stalled (abs. norm. error estimate became negative)',...
  'Selected coefficient tolerance reached',...
  'All zeros',...
  'Stalled (Residual has increased since last reset. Try to reduce resetit.)'...
  };
  info.message = status_str{1 + info.status};
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
classdef Conv3D < dagnn.Filter
  properties
    size = [0 0 0 0]
    hasBias = true
    opts = {'cuDNN'}
  end
  methods
    function outputs = forward(obj, inputs, params)
      if ~obj.hasBias, params{2} = [] ; end

      sz = size(inputs{1}); if numel(sz) < 4, sz(4) = 1; end

      nFrames = sz(4) / obj.net.meta.curBatchSize ;

      inputs{1} = reshape(inputs{1}, sz(1), sz(2), sz(3),  nFrames, sz(4) / nFrames ) ;
      inputs{1} = permute(inputs{1}, [1 2 4 3 5]);
      sz_in = size(inputs{1});

      outputs{1} = mex_conv3d(...
        inputs{1}, params{1}, params{2}, ...
        'pad', obj.pad, ...
        'stride', obj.stride) ;

      outputs{1} = permute(outputs{1}, [1 2 4 3 5]);
      sz_out = size(outputs{1});
      if numel(sz_out) == 4, sz_out(5) = 1; end % fixes the case of single batch
      nFrames = sz_out(4); % reset nframes
      outputs{1} = reshape(outputs{1}, sz_out(1), sz_out(2), sz_out(3),  nFrames*sz_out(5) ) ;

    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      if ~obj.hasBias, params{2} = [] ; end
      sz_in = size(inputs{1}); sz_out = size(derOutputs{1});
      if numel(sz_in) == 3, sz_in(4) = 1; end
      if numel(sz_out) == 3, sz_out(4) = 1; end
      nFramesIn = sz_in(4) / obj.net.meta.curBatchSize ;
      nFramesOut = sz_out(4) / obj.net.meta.curBatchSize ;

      inputs{1} = reshape(inputs{1}, sz_in(1), sz_in(2), sz_in(3),  nFramesIn, sz_in(4) / nFramesIn ) ;
      inputs{1} = permute(inputs{1}, [1 2 4 3 5]);
      derOutputs{1} = reshape(derOutputs{1}, sz_out(1), sz_out(2), sz_out(3),  nFramesOut, sz_out(4) / nFramesOut ) ;
      derOutputs{1} = permute(derOutputs{1}, [1 2 4 3 5]);
      [derInputs{1}, derParams{1}, derParams{2}] = mex_conv3d(...
        inputs{1}, params{1}, params{2}, derOutputs{1}, ...
        'pad', obj.pad, ...
        'stride', obj.stride) ;

      derInputs{1} = permute(derInputs{1}, [1 2 4 3 5]);
      sz_out = size(derInputs{1});
      if numel(sz_out) == 4, sz_out(5) = 1; end % fixes the case of single batch
      nFrames = sz_out(4); % reset nframes

      derInputs{1} = reshape(derInputs{1}, sz_out(1), sz_out(2), sz_out(3),  nFrames*sz_out(5) ) ;

    end

    function kernelSize = getKernelSize(obj)
      kernelSize = obj.size(1:2) ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      outputSizes = getOutputSizes@dagnn.Filter(obj, inputSizes) ;
      outputSizes{1}(3) = obj.size(4) ;
    end

    function params = initParams(obj)
      sc = sqrt(2 / prod(obj.size(1:3))) ;
      params{1} = randn(obj.size,'single') * sc ;
      params{2} = zeros(obj.size(4),1,'single') * sc ;
    end

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

    function rfs = getReceptiveFields(obj)
      ks = obj.getKernelSize() ;
      y1 = 1 - obj.pad(1) ;
      y2 = 1 - obj.pad(1) + ks(1) - 1 ;
      x1 = 1 - obj.pad(2) ;
      x2 = 1 - obj.pad(2) + ks(2) - 1 ;
      h = y2 - y1 + 1 ;
      w = x2 - x1 + 1 ;
      rfs.size = [h, w] ;
      rfs.stride = obj.stride(1:2) ;
      rfs.offset = [y1+y2, x1+x2]/2 ;
    end

  end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function mmap=filterjsonmmap(mmap, patterns, isinclude)
%
% mmap=filterjsonmmap(mmap, patterns, isinclude)
%
% filter JSON mmap keys based on inclusive or exclusive string patterns
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
% initially created on 2022/02/13
%
% input:
%      mmap: memory-map returned by loadjson/loadbj of the same data
%            important: mmap must be produced from the same file/string,
%            otherwise calling this function may cause data corruption
%      patterns: a string or a cell array of strings, each string will
%            be tested to match the JSONPath keys in mmap
%      isinclude: 1 (default) to include all mmap entries that match at
%            least one of the patterns, and 0 - exclude those that match
%
% output:
%      mmap: a filtered JSON mmap
%
% examples:
%      str='{"arr":[[1,2],"a",{"c":2}],"obj":{"k":"test"}}';
%      [dat, mmap]=loadjson(str);
%      savejson('',mmap)
%      newmmap=filterjsonmmap(mmap,{'arr.[1]', 'obj.k'});
%      savejson('',newmmap)
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%
    if(nargin<3)
        isinclude=1;
    end
    if(nargin>1 && ~isempty(patterns))
        keylist=[mmap{:}];
        keylist=keylist(1:2:end);
        if(~iscell(patterns))
            patterns={patterns};
        end
        mask=zeros(1,length(keylist));
        for i=1:length(patterns)
            mask=mask+cellfun(@length, strfind(keylist,patterns{i}));
        end
        if(isinclude)
            mmap=mmap(mask>0);
        else
            mmap(mask>0)=[];
        end
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [cl, imgdat] = extract_from_rois(imgs_to_extract_from, mask_image, varargin)
% Generic function for extracting image data from a mask or atlas image,
% and returning the data and averages within regions specified by the user.
%
% :Usage:
% ::
%
%    [cl, imgdat] = extract_from_rois(imgs_to_extract_from, mask_image, varargin)
%
%  Regions to average over can be either regions of contiguous voxels
%  bounded by voxels with values of 0 or NaN, which are considered non-data
%  values, or regions defined by unique integer codes in the mask image
%  (i.e., for atlas images with unique codes for each defined region.)
%
%  Mask/Atlas image does NOT have to be in the same space as the images to
%  extract from.  It will be remapped/resliced.
%
%  extracted data is returned in single data format.
%
% :Inputs:
%
%   **char array** of strings containing 4D image file names (data extracted from these)
%
%   **mask_image** to extract from
%
% :Optional Inputs:
%
%   **average_over:**
%        Default = 'unique_mask_values' to average over unique integer codes in the mask image
%        (i.e., for atlas images with unique codes for each defined region)
%
%        OPT = 'contiguous_regions' to average over contiguous voxels
%        bounded by voxels of 0 or NaN (non-data values)
%
%
% :Examples:
% ::
%
%    imgs_to_extract_from = filenames('w*.nii','char');
%    mask_image = which('anat_lbpa_thal.img');
%    [cl, imgdat] = extract_from_rois(imgs_to_extract_from, mask_image);
%
% : Programmer Notes:
% ::
%   Marianne corrected a mistaken variable name from dat to regiondat in
%   line 175. Error would go unnoticed if region dat was not empty (which
%   is normally the case). June 2016
%

    if isempty(mask_image)
        mask_image = which('anat_lbpa_thal.img');
    end

    average_over = 'unique_mask_values'; %'contiguous_regions'  or 'unique_mask_values';

    for varg = 1:length(varargin)
        if ischar(varargin{varg})
            switch varargin{varg}

                % reserved keywords
                case 'contiguous_regions', average_over = 'contiguous_regions';
                case 'unique_mask_values', average_over = 'unique_mask_values';
            end
        end
    end

    space_defining_image = deblank(imgs_to_extract_from(1, :));

    % Note: we need to write out the image here only because we need its
    % anatomical boundaries (in-mask areas), contig voxels, etc. for later

    tmpname = num2str(round(10*rand(1, 5))); tmpname(tmpname == ' ') = []; tmpname = ['tmp_mask_' tmpname '.img'];

    maskData = scn_map_image(mask_image, space_defining_image, 'write', tmpname);
    maskData = maskData(:);

    volInfo = iimg_read_img(tmpname, 2);
    maskData = maskData(volInfo.wh_inmask);

    % This needs to be made Windows-compatible. kludgy now.
    eval(['!rm ' tmpname])
    eval(['!rm ' tmpname(1:end-4) '.hdr'])

    % Now we have the mask data and volInfo structure, and we can extract

    %%
    clear imgdat
    switch spm('Ver')

        case {'SPM12', 'SPM8', 'SPM5'}

            imgdat = iimg_get_data(volInfo, imgs_to_extract_from, 'single', 'verbose', 'noexpand');

        case {'SPM2', 'SPM99'}
            % legacy, for old SPM

            imgdat = iimg_get_data(volInfo, imgs_to_extract_from, 'single', 'verbose');

        otherwise
            error('Unknown version of SPM! Update code, check path, etc.');
    end

    %% Now get averages by cluster, if requested

    switch average_over

        % Define integer codes for sets of voxels to average over.

        case 'unique_mask_values'
            maskData = round(maskData);
            u = unique(maskData)'; u(u == 0) = [];
            nregions = length(u);
            fprintf('Averaging over unique mask values, assuming integer-valued mask: %3.0f regions\n', nregions);

        case 'contiguous_regions'
            u = unique(volInfo.cluster); u(u == 0) = [];
            maskData = volInfo.cluster;

        case 'none'

        otherwise
            error('Illegal value for average_over.  See help for this function.');
    end

    %% Now get the average activity in each region and define a "cluster"

    nregions = length(u);

    cl = struct('title', 'title', 'threshold', NaN, 'Z', NaN, 'voxSize', abs(diag(volInfo.mat(1:3, 1:3)))', ...
        'XYZ', [0 0 0]', 'XYZmm', [0 0 0]', 'from_label', NaN, 'M', volInfo.mat, 'dim', volInfo.dim);

    cl(1:nregions) = cl;

    for i = 1:nregions
        imgvec = maskData == u(i);

        regiondat = imgdat(:, imgvec);

        if ~isempty(regiondat)
            regionmean = double(nanmean(regiondat')');

        else
            regionmean = NaN .* zeros(size(regiondat, 1), 1);
            disp('WARNING: Empty region dat, check your masks');
        end

        cl(i).average_data = regionmean;
        cl(i).XYZ = volInfo.xyzlist(imgvec,:)';
        cl(i).XYZmm = voxel2mm(cl(i).XYZ,volInfo.mat);

        cl(i).Z = ones(1,size(cl(i).XYZ,2));

    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function  gB = reorder(gB,varargin)
% nicely reorder grain boundaries

% compute tours through the adjecency matrix
[idV,idE] = EulerTours(gB.A_V);

% remove breaks
id = isnan(idE);
idE(id) = [];
idV(id) = [];

% reorder edges
F = gB.F(idE,:);

% and also all other properties
gB.ebsdId = gB.ebsdId(idE,:);
gB.grainId = gB.grainId(idE,:);
gB.phaseId = gB.phaseId(idE,:);
gB.misrotation = gB.misrotation(idE);

% which needs to be flipped?
flip = F(:,1) ~= idV(:);

% flip vertices
F(flip,:) = fliplr(F(flip,:));
gB.F = F;

author: mtex-toolbox

repo: mtex

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


text:
function res = snr(hat, star)

    sigma = std(reshape(star,size(star,1)*size(star,2),1));
    res = psnr(hat, star, sigma);

end


text:
function res = pickconditions(this, label, rejectbad)
% Method for returning indices of trials of a certain trial type.
% If input argument rejectbad == 0, the function will also return trial
% indices which are set to bad (i.e. rejected). If bad is omitted,
% the default is not to include rejected trials.
% FORMAT res = pickconditions(this)
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: pickconditions.m 3146 2009-05-26 09:54:23Z vladimir $

if nargin<3
    rejectbad = 1;
end

c = conditions(this);

if isa(label, 'char')
    label = {label};
end

res = find(ismember(c, label));

if rejectbad && ~isempty(res)
    res = res(~reject(this, res));
end


text:
function [oscilatior_corr_ifgs,oscilatior_corr_velocity] = env_oscilator_corr(envisat_flag,forced_sm_flag)
% This function perform oscialtor drift correction for envisat interferograms based on
% Petar Marinkovic presentation at ESA Living Planet 2013 in Edinburgh.
% Approximation formula:
% dR/year = c/2 * (slantRangeTime_FAR - slantRangeTime_NEAR) * corrPerYear
% 'Apparent displacement' correction:
% dR/year ~ (7.8m * 5000)*3.87e-7 ~ 0.01482m
% Correction uses the pixel information in range and in adition for
% interferograms the temporal baseline information
%
% OPTIONAL INPUT:
% envisat_flag                  'y' when envisat or 'n'. When empty, a
%                               search is performed in the master.res file
%
% OUTPUTS:
% oscilatior_corr_velocity      Correction in mm for the velocity
% oscilatior_corr_ifgs          Correction for individual interferograms in rad
%
% Correction is defined such that:
% Corrected interferogram/velocity = original interferogram/velocity - correction
%
% P. Marinkovic Envisat oscialtor drift correction coded by David Bekaert -- University of Leeds 2014
%
% cite as:
% P. Marinkovic (PPO.labs) and Y. Larsen (NORUT)
% Consequences of Long-Term ASAR Local Oscillator Frequency Decay - an Empirical Study of 10 Years of Data
% ESA Living Planet Symposium (2013)
%
% Modifications
% 04/2014       DB      Put non-envisat fix to zeros
% 04/2014       DB      Allow forced SM oscialtor drift computation
% 06/2014       DB      Fixed error for SM computation and added extra envisat check.
% 06/2014       DB      Fix in case not envisat and forced SM.
% 11/2014       DB      Fix to make windows and linux compatible
% 03/2015       DB      Clean script output
% 01/2017       DB      Bug fix for non-envisat SM case. n_ifg was not n_image

if nargin<1 || isempty(envisat_flag)
  % checking if this is envisat or not
  platform=getparm('platform');
  if isempty(platform)
    if exist('master.res','file')==2
        master_file = 'master.res';
    elseif exist('../master.res','file')==2
        master_file = '../master.res';
    else
        master_file = [];
    end

    if ~isempty(master_file)

        % make windows and linux compatible
        a = fileread(master_file);
        ix = strfind(a,'ASAR');
       if ~isempty(ix)
           platform='ENVISAT';
       end
    else
        fprintf('Could not check if this is Envisat \n')
    end
  end
  if strcmpi(platform,'ENVISAT')
      envisat_flag = 'y';
      fprintf('This is Envisat, oscilator drift is being removed... \n')
  else
       envisat_flag = 'n';
  end

end

small_baseline_flag = getparm('small_baseline_flag');
if nargin<2
    forced_sm_flag=0;
end

if forced_sm_flag==1
    small_baseline_flag='n';
end

if strcmp(envisat_flag,'y')
    load psver
    % use the ps2.mat data
    ps = load(['ps' num2str(psver) '.mat']);
    lambda = getparm('lambda');

    % velocity map correction:
    envisat_resolution = 7.8;                   % ground range resolution [m]
    Oscilator_drift_corr_year = 3.87e-7;         % drift correction in range [1/year]

    % velocity correction in mm
    oscilatior_corr_velocity = (envisat_resolution*ps.ij(:,3))*Oscilator_drift_corr_year*1000;

    % interferogram
    if strcmp(small_baseline_flag,'y')
        n_ifg = ps.n_ifg;
        delta_year = (ps.ifgday(:,2)-ps.ifgday(:,1))./365.25;
    else
        n_ifg = ps.n_image;
        delta_year = (ps.day-ps.master_day)./365.25;
    end
    oscilatior_corr_ifgs = -4.*pi./lambda.*repmat(oscilatior_corr_velocity,1,n_ifg)./1000.*repmat(delta_year',ps.n_ps,1);

else
    load psver
    ps = load(['ps' num2str(psver) '.mat']);

    if strcmp(small_baseline_flag,'y')
        oscilatior_corr_ifgs = zeros([ps.n_ps ps.n_ifg]);
    else
        n_ifg = ps.n_image;     % bug fix DB
        oscilatior_corr_ifgs = zeros([ps.n_ps n_ifg]);
    end
    oscilatior_corr_velocity = zeros([ps.n_ps 1]);

end

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function bnet = mk_map_hhmm(varargin)

% p is the prob of a successful move (defines the reliability of motors)
p = 1;
obs_model = 'unique';

for i=1:2:length(varargin)
  switch varargin{i},
   case 'p', p = varargin{i+1};
   case 'obs_model', obs_model = varargin{i+1};
  end
end

q = 1-p;
unique_obs = strcmp(obs_model, 'unique');

% assign numbers to the nodes in topological order
U = 1; A = 2; C = 3; F = 4;
if unique_obs
  onodes = 5;
else
  N = 5; E = 6; S = 7; W = 8; % north, east, south, west
  onodes = [N E S W];
end

% create graph structure

ss = 4 + length(onodes); % slice size
intra = zeros(ss,ss);
intra(U,F)=1;
intra(A,[C F onodes])=1;
intra(C,[F onodes])=1;

inter = zeros(ss,ss);
inter(U,[A C])=1;
inter(A,[A C])=1;
inter(F,[A C])=1;
inter(C,C)=1;

% node sizes
ns = zeros(1,ss);
ns(U) = 2; % left/right
ns(A) = 2;
ns(C) = 3;
ns(F) = 2;
if unique_obs
  ns(onodes) = 5; % we will assign each state a unique symbol
else
  ns(onodes) = 2;
end
l = 1; r = 2; % left/right
L = 1; R = 2;

% Make the DBN
bnet = mk_dbn(intra, inter, ns, 'observed', onodes);
eclass = bnet.equiv_class;

% Define CPDs for slice 1
% We clamp all the CPDs that are not tied,
% since we cannot learn them from a single sequence.

% uniform probs over actions (the input could be chosen from a policy)
bnet.CPD{eclass(U,1)} = tabular_CPD(bnet, U, 'CPT', mk_stochastic(ones(ns(U),1)), ...
				    'adjustable', 0);

% uniform probs over starting abstract state
bnet.CPD{eclass(A,1)} = tabular_CPD(bnet, A, 'CPT', mk_stochastic(ones(ns(A),1)), ...
				    'adjustable', 0);

% Uniform probs over starting concrete state, modulo the fact
% that corridor 2 is only of length 2.
CPT = zeros(ns(A), ns(C)); % CPT(i,j) = P(C starts in j | A=i)
CPT(1, :) = [1/3 1/3 1/3];
CPT(2, :) = [1/2 1/2 0];
bnet.CPD{eclass(C,1)} = tabular_CPD(bnet, C, 'CPT', CPT, 'adjustable', 0);

% Termination probs
CPT = zeros(ns(U), ns(A), ns(C), ns(F));
CPT(r,1,1,:) = [1 0];
CPT(r,1,2,:) = [1 0];
CPT(r,1,3,:) = [q p];
CPT(r,2,1,:) = [1 0];
CPT(r,2,2,:) = [q p];
CPT(l,1,1,:) = [q p];
CPT(l,1,2,:) = [1 0];
CPT(l,1,3,:) = [1 0];
CPT(l,2,1,:) = [q p];
CPT(l,2,2,:) = [1 0];

bnet.CPD{eclass(F,1)} = tabular_CPD(bnet, F, 'CPT', CPT);

% Observation model
if unique_obs
  CPT = zeros(ns(A), ns(C), 5);
  CPT(1,1,1)=1;  % Theo state 4
  CPT(1,2,2)=1;  % Theo state 5
  CPT(1,3,3)=1; % Theo state 6
  CPT(2,1,4)=1; % Theo state 9
  CPT(2,2,5)=1; % Theo state 10
  %CPT(2,3,:) undefined
  O = onodes(1);
  bnet.CPD{eclass(O,1)} = tabular_CPD(bnet, O, 'CPT', CPT);
else
  % north/east/south/west can see wall (1) or opening (2)
  CPT = zeros(ns(A), ns(C), 2);
  CPT(:,:,1) = q;
  CPT(:,:,2) = p;
  bnet.CPD{eclass(W,1)} = tabular_CPD(bnet, W, 'CPT', CPT);
  bnet.CPD{eclass(E,1)} = tabular_CPD(bnet, E, 'CPT', CPT);
  CPT = zeros(ns(A), ns(C), 2);
  CPT(:,:,1) = p;
  CPT(:,:,2) = q;
  bnet.CPD{eclass(S,1)} = tabular_CPD(bnet, S, 'CPT', CPT);
  bnet.CPD{eclass(N,1)} = tabular_CPD(bnet, N, 'CPT', CPT);
end

% Define the CPDs for slice 2

% Abstract

% Since the top level never resets, the starting distribution is irrelevant:
% A2 will be determined by sampling from transmat(A1,:).
% But the code requires we specify it anyway; we make it all 0s, a dummy value.
startprob = zeros(ns(U), ns(A));

transmat = zeros(ns(U), ns(A), ns(A));
transmat(R,1,:) = [q p];
transmat(R,2,:) = [0 1];
transmat(L,1,:) = [1 0];
transmat(L,2,:) = [p q];

% Qps are the parents we condition the parameters on, in this case just
% the past action.
bnet.CPD{eclass(A,2)} = hhmm2Q_CPD(bnet, A+ss, 'Fbelow', F, ...
				  'startprob', startprob, 'transprob', transmat);

% Concrete

transmat = zeros(ns(C), ns(U), ns(A), ns(C));
transmat(1,r,1,:) = [q p 0.0];
transmat(2,r,1,:) = [0.0 q p];
transmat(3,r,1,:) = [0.0 0.0 1.0];
transmat(1,r,2,:) = [q p 0.0];
transmat(2,r,2,:) = [0.0 1.0 0.0];
%
transmat(1,l,1,:) = [1.0 0.0 0.0];
transmat(2,l,1,:) = [p q 0.0];
transmat(3,l,1,:) = [0.0 p q];
transmat(1,l,2,:) = [1.0 0.0 0.0];
transmat(2,l,2,:) = [p q 0.0];

% Add a new dimension for A(t-1), by copying old vals,
% so the matrix is the same size as startprob

transmat = reshape(transmat, [ns(C) ns(U) ns(A) 1 ns(C)]);
transmat = repmat(transmat, [1 1 1 ns(A) 1]);

% startprob(C(t-1), U(t-1), A(t-1), A(t), C(t))
startprob = zeros(ns(C), ns(U), ns(A), ns(A), ns(C));
startprob(1,L,1,1,:) = [1.0 0.0 0.0];
startprob(3,R,1,2,:) = [1.0 0.0 0.0];
startprob(3,R,1,1,:) = [0.0 0.0 1.0];
%
startprob(1,L,2,1,:) = [0.0 0.0 010];
startprob(2,L,2,1,:) = [1.0 0.0 0.0];
startprob(2,R,2,2,:) = [0.0 1.0 0.0];

% want transmat(U,A,C,At,Ct), ie. in topo order
transmat = permute(transmat, [2 3 1 4 5]);
startprob  = permute(startprob, [2 3 1 4 5]);
bnet.CPD{eclass(C,2)} = hhmm2Q_CPD(bnet, C+ss, 'Fself', F, ...
				  'startprob', startprob, 'transprob', transmat);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function pass = test_sample( )
% Test diskfun sample() command

tol = 100*chebfunpref().cheb2Prefs.chebfun2eps;

% Function to test
f = diskfun(@(x,y) sin(pi*x.*y));

% Ensure the matrix of sampled values is correct.
[m,n] = length(f);
[nn,mm] = size(sample(f));
pass(1) = (m == mm) && (n == nn);

% Sample on fixed grids of various sizes to make sure the right size output
% is given.
m = 120;
n = 121;
[nn,mm] = size(sample(f, m, n));
pass(2) = (m == mm) && (n == nn);

m = 121;
n = 120;
[nn,mm] = size(sample(f, m, n));
pass(3) = (m == mm) && (n == nn);

% Check samples are correct.
% m even and n odd
m = 30;
n = 2*20+1;
cp = chebpts(n);
[t,r] = meshgrid(trigpts(m, [-pi, pi]), cp((n+1)/2:end));
F = f(t,r, 'polar');
G = sample(f, m, (n+1)/2);
pass(4) = norm(F(:) - G(:), inf) < tol;
[U, D, V] = sample(f, m, (n+1)/2);
G = U * D * V.';
pass(5) = norm(F(:) - G(:), inf) < tol;

% m odd and n odd
m = 31;
n = 2*20+1;
cp = chebpts(n);
[t,r] = meshgrid(trigpts(m, [-pi, pi]), cp((n+1)/2:end));
F = f(t,r, 'polar');
G = sample(f, m, (n+1)/2);
pass(6) = norm(F(:) - G(:), inf) < tol;
[U, D, V] = sample(f, m, (n+1)/2);
G = U * D * V.';
pass(7) = norm(F(:) - G(:), inf) < tol;

% m odd and n even
m = 31;
n = 2*20-1;
cp = chebpts(n);
[t,r] = meshgrid(trigpts(m, [-pi, pi]), cp((n+1)/2:end));
F = f(t,r, 'polar');
G = sample(f, m, (n+1)/2);
pass(8) = norm(F(:) - G(:), inf) < tol;
[U, D, V] = sample(f, m, (n+1)/2);
G = U * D * V.';
pass(9) = norm(F(:) - G(:), inf) < tol;

% m even and n even
m = 30;
n = 2*20-1;
cp = chebpts(n);
[t,r] = meshgrid(trigpts(m, [-pi, pi]), cp((n+1)/2:end));
F = f(t,r, 'polar');
G = sample(f, m, (n+1)/2);
pass(10) = norm(F(:) - G(:), inf) < tol;
[U, D, V] = sample(f, m, (n+1)/2);
G = U * D * V.';
pass(11) = norm(F(:) - G(:), inf) < tol;

% Sample should return all ones for the function 1.
f = diskfun(@(x,y) 1 + 0*x);
F = sample(f, 128, 128);
pass(12) = norm(F(:) - 1, inf) < tol;

% Check that errors are caught
try
    F = sample(f, 0, 20);
    pass(13) = false;
catch ME
    pass(13) = strcmp(ME.identifier, 'CHEBFUN:DISKFUN:sample:inputs');
end

try
    F = sample(f, 20, 0);
    pass(14) = false;
catch ME
    pass(14) = strcmp(ME.identifier, 'CHEBFUN:DISKFUN:sample:inputs');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function node = meshsmoothing(node,elem,step,rho,method,fixedNode)
%% MESHSMOOTHING improves the geometric mesh quality
%
% node = meshsmoothing(node,elem) improves shape regularity of
% triangles while keeping the density of vertices.
%
% node = meshsmoothing(node,elem,m) performs m steps mesh smoothing.
% The default setting is m=3.
%
% node = meshsmoothing(node,elem,m,rho) accept a non-uniform density
% given by elementwise function rho. The default choice rho = 1/|t|. The
% quasi-uniform grids corresponds to rho=1. The density function rho can be
% given by a user specified function or a posteriori error estimator in the
% setting of adaptive finite element method.
%
% node = meshsmoothing(node,elem,m,rho,method,fixedNode) accept the
% extra input arguments method and fixedNode. The string |method| is to
% specify the smoothing method: either 'CPT' or 'ODT'. The fixedNode array
% will fix certain nodes, e.g., vertices on a interface. The boundary nodes
% of the mesh will be automatically considered as fixed.
%
% The function meshsmoothing will keep the topology of the input mesh,
% i.e., the node index and connectivity of nodes are unchanged. In
% contrast, the edgeswap function will keep the node but change elem.
%
% The algorithm implemented is the simplified version of ODT smoothing; see
% ODTmesh.
%
% Example:
%     load airfoilperturbmesh
%     meshquality(node,elem);
%     node = meshsmoothing(node,elem);
%     meshquality(node,elem);
%
% See also  optmesh, bdsmoothing, edgeswap, rmisopoint
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('step','var'), step = 3; end
if ~exist('method','var'), method = 'ODT'; end

%% Smooth boundary nodes
% [node,elem,bdNode] = bdsmoothing(node,elem);

%% Find boundary and interior ndoes
[bdNode,bdEdge,isBdNode,isBdElem] = findboundary(elem); %#ok<*ASGLU>
N = size(node,1);  NT = size(elem,1);

%% Compute mesh quality associated to elements and nodes
qt = 1 - simpqual(node,elem);
qtp = sparse([1:NT,1:NT,1:NT], elem(1:NT,:), [qt, qt, qt], NT, N);
qp = max(qtp);
qp = qp' + 0.01*rand(N,1);  % break the equal case
qp(bdNode) = 0;             % do not include bdNode and fixed nodes
if exist('fixedNode','var'), qp(fixedNode) = 0; end

%% Decompose nodes into independent sets
t2p = sparse([1:NT,1:NT,1:NT], elem(1:NT,:), 1, NT, N);
A = t2p'*t2p;  % connectness of nodes
nodeSet = coloring(A,5,qp);

%% Move nodes one by one
for m = 1:step
    for k = 1:5
        % find elements containing nodeSet
        movingNode = nodeSet{k};
        if isempty(movingNode)
            break;
        end
        [i,j] = find(t2p(:,movingNode)); %#ok<*NASGU>
        idx = false(NT,1);
        idx(i) = true;
        idx = find(idx);
        % compute centers and areas
        switch upper(method)
          case 'CPT'
               center = (node(elem(idx,1),:)+node(elem(idx,2),:)+node(elem(idx,3),:))/3;
          case 'ODT'
               center = circumcenter(node,elem(idx,:));
               % modification neary boundary elements: using barycenter to replace circumcenter
               center(isBdElem(idx),:) = (node(elem(idx(isBdElem(idx)),1),:) ...
                                        + node(elem(idx(isBdElem(idx)),2),:) ...
                                        + node(elem(idx(isBdElem(idx)),3),:))/3;
        end
        % compute the weight
        if exist('rho','var') && ~isempty(rho)
            weight = simplexvolume(node,elem(idx,:));
            if (rho~=1)
                weight = rho(idx).*weight;
            end
        end
        % update to new averaging center
        oldpi = node(movingNode,:);     % keep the location before moving
        if exist('rho','var') && ~isempty(rho) % weighted average
            newnode(:,1) = accumarray([elem(idx,1);elem(idx,2);elem(idx,3)],...
                                repmat(weight.*center(:,1),3,1),[N,1]);
            newnode(:,2) = accumarray([elem(idx,1);elem(idx,2);elem(idx,3)],...
                                repmat(weight.*center(:,2),3,1),[N,1]);
            valence = accumarray([elem(idx,1);elem(idx,2);elem(idx,3)], ...
                                 [weight; weight; weight],[N 1]);
        else % simple average of centers
            newnode(:,1) = accumarray([elem(idx,1);elem(idx,2);elem(idx,3)],...
                                      [center(:,1);center(:,1);center(:,1)],[N,1]);
            newnode(:,2) = accumarray([elem(idx,1);elem(idx,2);elem(idx,3)],...
                                      [center(:,2);center(:,2);center(:,2)],[N,1]);
            valence = accumarray([elem(idx,1);elem(idx,2);elem(idx,3)],...
                                  ones(3*length(idx),1),[N 1]);
        end
        node(movingNode,:) = newnode(movingNode,:)./[valence(movingNode) valence(movingNode)];
        % check if the moving is valid
        ve2 = node(elem(idx,1),:)-node(elem(idx,3),:);
        ve3 = node(elem(idx,2),:)-node(elem(idx,1),:);
        area = 0.5*(-ve3(:,1).*ve2(:,2)+ve3(:,2).*ve2(:,1));
        invalidElem = idx(area<0);
        if any(invalidElem)
            isNotmoving = false(N,1);
            isNotmoving(elem(invalidElem,:)) = true;
            pidx = isNotmoving(movingNode);
            node(movingNode(pidx),:) = oldpi(pidx,:);
        end
    end
end

repo: ifem


text:
% ------------------------------------------------------------------------
%  Copyright (C)
%  Torr Vision Group (TVG)
%  University of Oxford - UK
%
%  Qizhu Li <liqizhu@robots.ox.ac.uk>
%  August 2018
% ------------------------------------------------------------------------
% This file is part of the weakly-supervised training method presented in:
%    Qizhu Li*, Anurag Arnab*, Philip H.S. Torr,
%    "Weakly- and Semi-Supervised Panoptic Segmentation,"
%    European Conference on Computer Vision (ECCV) 2018.
% Please consider citing the paper if you use this code.
% ------------------------------------------------------------------------
% This function visualises the original prediction, original prediction
% scores, processed semantic iterative GT, and generated instance
% iterative GT.
%  INPUT:
%  - opts : options generated by get_opts.m
%  - results : results generated by the preceding code
%
%  DEMO:
%  - See scripts/run_sub.m
% ------------------------------------------------------------------------

function visualise_results_cam_mandg(opts, results)

if ~opts.visualise_results
    return;
end

f = figure('units','normalized','outerposition',[0 0 1 1], 'Visible', 'off');

ax(1) = subplot(2, 2, 1, 'replace');
imshow(results.pred, opts.cmap);
c(1) = colorbar(ax(1));
c(1).Limits = [0, 19];
c(1).Ticks = 0.5:1:18.5;
c(1).TickLabels = opts.objectNames;
c(1).FontSize = 8;
title(ax(1), '1. CAM prediction');

ax(2) = subplot(2, 2, 2, 'replace');
imshow(results.mandg_pred, results.mandg_cmap);
title(ax(2), '2. MCG&Grabcut Masks');

ax(3) = subplot(2, 2, 3, 'replace');
imshow(results.final_pred, opts.cmap);
c(3) = colorbar(ax(3));
c(3).Limits = [0, 19];
c(3).Ticks = 0.5:1:18.5;
c(3).TickLabels = opts.objectNames;
c(3).FontSize = 8;
title(ax(3), '3. Processed semantic GT');

ax(4) = subplot(2, 2, 4, 'replace');
imshow(results.ins_pred, opts.cmap);
c(4) = colorbar(ax(4));
num_instances = numel(results.ins_info);
c(4).Limits = [1, num_instances + 1];
c(4).Ticks = 1.5:1:(num_instances+0.5);
c(4).TickLabels = strcat(strsplit(num2str(extractfield(results.ins_info, 'ins_id'))), ...
    {': '}, extractfield(results.ins_info, 'name'));
c(4).FontSize = 8;
title(ax(4), '2. Generated instance GT');

f.Visible = 'on';

end

author: qizhuli

repo: Weakly-Supervised-Panoptic-Segmentation


text:
function s = ymd_to_s_alexandrian ( y, m, d )

%*****************************************************************************80
%
%% YMD_TO_S_ALEXANDRIAN "prints" an Alexandrian YMD date into a string.
%
%  Format:
%
%    DayNumber MonthName YearNumber AX
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, M, D, the YMD date.
%
%    Output, string S, a representation of the date.
%
  sm = month_to_month_name_eg_civil ( m );

  s = sprintf ( '%d %s %d AX', d, sm, y );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [bb] = tribal2(pp,tt)
%TRIBAL2 compute the circumballs associated with a 2-simplex
%triangulation embedded in R^2 or R^3.
%   [BB] = TRIBAL2(PP,TT) returns the circumscribing balls
%   associated with the triangles in [PP,TT], such that BB =
%   [XC,YC,RC.^2].

%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 02/05/2018

    bb = pwrbal2(pp,zeros(size(pp,1),1),tt) ;

end

author: CHLNDDEV


text:
%{
    Copyright (C) Cambridge Electronic Design Limited 2014
    Author: James Thompson
    Web: www.ced.co.uk email: james@ced.co.uk, softhelp@ced.co.uk

    This file is part of CEDS64ML, a MATLAB interface to the SON64 library.

    CEDS64ML 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.

    CEDS64ML 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 CEDS64ML.  If not, see <http://www.gnu.org/licenses/>.
%}

function [ iOk, TimeDateOut ] = CEDS64TimeDate( fhand, TimeDateIn )
%CEDS64TIMEDATE Gets and sets the time and date saved in a file (usually the time that sampling began)
% Note this is not the same format as MATLAB's clock function
%   [ iOk, TimeDateOut ] = CEDS64TimeDate( fhand, TimeDateIn )
%   Inputs
%   fhand - Integer file handle
%   TimeDateIn - (Optional) The new time date as a 7-by-1 vector
%   Outputs
%   iOk = 0 if the operation completed without error, otherwise a negative
%   error code.
%   TimeDateOut - The old time date as a 7-by-1 vector
%   The structure of the time-date vector is:
%   element 1 - hundredths of seconds (0-99)
%   element 2 - seconds (0-59)
%   element 3 - minutes (0-59)
%   element 4 - hours (0-23)
%   element 5 - day of the month (1-31)
%   element 6 - month of the year (1-12)
%   element 7 - year (1980-2200)

%TimeDateOut = zeros(7, 1, 'int32');

if (nargin == 1 || nargin == 2) % get the old TimeDate
    dBufferIn = int64(0);
    dBufferOut = int64(0);
    [ iOk, A ] = calllib('ceds64int', 'S64TimeDate', fhand, dBufferOut, dBufferIn, -1);
    TimeDateOut(1) = bitand(bitshift(A, 0), 255);     % Hundredth of a second
    TimeDateOut(2) = bitand(bitshift(A, -8), 255);    % Second
    TimeDateOut(3) = bitand(bitshift(A, -16), 255);   % Minute
    TimeDateOut(4) = bitand(bitshift(A, -24), 255);   % Hour
    TimeDateOut(5) = bitand(bitshift(A, -32), 255);   % day
    TimeDateOut(6) = bitand(bitshift(A, -40), 255);   % month
    TimeDateOut(7) = bitand(bitshift(A, -48), 65535); % year
else
   iOk = -22;
   return;
end

% set the new TimeDate if we're given one
if (nargin == 2)
    if (length(TimeDateIn) ~= 7)
        iOk = -22;
        return;
    end
	dBufferIn = int64(0);
    dBufferOut = int64(0);
    dBufferIn = dBufferIn + TimeDateIn(1);
    dBufferIn = dBufferIn + bitshift(TimeDateIn(2), 8);
    dBufferIn = dBufferIn + bitshift(TimeDateIn(3), 16);
    dBufferIn = dBufferIn + bitshift(TimeDateIn(4), 24);
    dBufferIn = dBufferIn + bitshift(TimeDateIn(5), 32);
    dBufferIn = dBufferIn + bitshift(TimeDateIn(6), 40);
    dBufferIn = dBufferIn + bitshift(TimeDateIn(7), 48);
    [ iOk, A ] = calllib('ceds64int', 'S64TimeDate', fhand, dBufferOut, dBufferIn, 0);
end
end

author: brainstorm-tools


text:
function matches = dir_regex(path, expression, return_paths)
%DIR_REGEX Returns contents in the path matching a regular expression.
% Usage:
%   matches = dir_regex(path, expression)
%   matches = dir_regex(path, expression, true) % returns full paths
%
% See also: regexp, dir_ext, dir_paths

% Default: returns just filenames
if nargin < 3; return_paths = false; end

% Get directory listing
listing = dir_paths(path);

% Test regex
matches = listing(~areempty(regexp(listing, expression)));

% Singleton output
if iscell(matches) && numel(matches) == 1
    matches = matches{1};
end

% Return just the filenames
if ~return_paths
    matches = get_filename(matches);
end

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


text:
function [s, m_v, zero_idx] = get_central_ula_size(diffs)
%GET_CENTRAL_ULA_SIZE Obtains the size of the ULA centered at the origin
%from the difference coarray.
%Syntax:
%   s = GET_CENTRAL_ULA_SIZE(diffs);
%   [s, m_v, zero_idx] = GET_CENTRAL_ULA_SIZE(diffs);
%Input:
%   diffs - A vector of integer differences, or array design.
%Outputs:
%   s - Entire array size = 2*M_v - 1.
%   m_v - The size of the ULA after trimming the negative part.
%   zero_idx - Index of the zero element in diffs.

if isstruct(diffs)
    % convert array design to unique differences
    diffs = unique_differences(diffs.element_indices);
end

w = 0;
zero_idx = find(diffs == 0);
n_diff = length(diffs);
while zero_idx + w < n_diff && zero_idx - w > 1
    w = w + 1;
    if (diffs(zero_idx+w) - diffs(zero_idx+w-1) ~= 1) || (diffs(zero_idx-w+1) - diffs(zero_idx-w)) ~= 1
        w = w - 1;
        break;
    end
end
m_v = w + 1;
s = 2*m_v - 1;

end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function value = r8_power ( r, p )

%*****************************************************************************80
%
%% R8_POWER computes the P-th power of R.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, the base.
%
%    Input, integer P, the power, which may be negative.
%
%    Output, real VALUE, the value of R**P.
%
  mults = 0;
%
%  Force P to be an integer.
%
  p = floor ( p );

%
%  Special case.  R^0 = 1.
%
  if ( p == 0 )
    value = 1.0;
%
%  Special case.  Positive powers of 0 are 0.
%  For negative powers, we go ahead and compute it, hoping software will complain.
%
  elseif ( r == 0.0 )
    if ( 0 < p )
      value = 0.0;
    else
      value = r^p;
    end
  elseif ( 1 <= p )
    value = r^p;
  else
    value = r^p;
  end

  return
end


text:
function mr = mrCreateEmpty;
% Create an empty mr struct.
%
% mr = mrCreateEmpty;
%
% ras 07/05.
mr.data = [];
mr.name = '';
mr.path = '';
mr.hdr = [];
mr.info.scanner = '';
mr.info.subject = '';
mr.info.subjectSex = '';
mr.info.subjectAge = [];
mr.info.date = '';
mr.info.scanStart = [];
mr.info.examNum = [];
mr.info.protocol = '';
mr.info.coil = '';
mr.format = '';
mr.spaces = [];
mr.voxelSize = [];
mr.dims = [];
mr.extent = [];
mr.dimUnits = '';
mr.dataUnits = '';
mr.dataRange = [];
mr.phaseFlag = 0;  % by default, set to 1 otherwise
mr.params = [];
mr.comments = '';

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function V = LR_sort(W)
N = size(W, 2);
LLR = zeros(1, N);
for i = 1 : N
    if (W(1, i) ~= 0) && (W(2, i) ~= 0)
        LLR(i) = log(W(1, i)) - log(W(2, i));
    else
        if (W(1, i) == 0) && (W(2, i) ~= 0)
            LLR(i) = -inf;
        else
            if (W(1, i) ~= 0) && (W(2, i) == 0)
                LLR(i) = inf;
            end
        end
    end
end
[~, ordered]  = sort(LLR, 'descend');
V = W(:, ordered);
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
% Fig. 6.15   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
%close all;
clf

w=logspace(-2,2,100);
k=10;
num=k;
den=conv([1 0],[1 2 1]);
[m10,p10]=bode(num,den,w);
k=2;
num=k;
den=conv([1 0],[1 2 1]);
[m2,p2]=bode(num,den,w);
k=0.1;
num=k;
den=conv([1 0],[1 2 1]);
[mp1,pp1]=bode(num,den,w);
figure(1)
loglog(w,m10,w,m2,w,mp1,w,ones(size(w)));
xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.15 Frequency Response (a) magnitude');
text(1,20,'K=10');
text(0.2,3,'K=2');
text(0.9,0.01,'K=0.1');
bodegrid;
%pause;
figure(2)
semilogx(w,p10,w,p2,w,pp1,w,-180*ones(size(w)),'-');
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)')
title('Fig. 6.15 (b) phase');
bodegrid;

author: Sable

repo: mcbench-benchmarks


text:
function [beta_median beta_std beta_lbound beta_ubound sigma_median sigma_t_median sigma_t_lbound sigma_t_ubound]=stvol3estimates(beta_gibbs,sigma_gibbs,sigma_t_gibbs,n,T,cband)

% compute the median, variance, and credibility intervals for the posterior distribution of beta
beta_median=quantile(beta_gibbs,0.5,2);
beta_std=std(beta_gibbs,0,2);
beta_lbound=quantile(beta_gibbs,(1-cband)/2,2);
beta_ubound=quantile(beta_gibbs,1-(1-cband)/2,2);

% compute the results for sigma (long-run value)
sigma_median=reshape(quantile(sigma_gibbs,0.5,2),n,n);

% compute the rsults for sigma (sample values)
sigma_t_median=cell(n,n);
sigma_t_lbound=cell(n,n);
sigma_t_ubound=cell(n,n);
% loop over periods and entries
for ii=1:T
   for jj=1:n
      for kk=1:jj
      sigma_t_median{jj,kk}(ii,1)=quantile(sigma_t_gibbs{ii,1}(jj,kk,:),0.5,3);
      sigma_t_lbound{jj,kk}(ii,1)=quantile(sigma_t_gibbs{ii,1}(jj,kk,:),(1-cband)/2,3);
      sigma_t_ubound{jj,kk}(ii,1)=quantile(sigma_t_gibbs{ii,1}(jj,kk,:),1-(1-cband)/2,3);
      end
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function bnet = mk_limid(dag, node_sizes, varargin)
% MK_LIMID Make a limited information influence diagram
%
% BNET = MK_LIMID(DAG, NODE_SIZES, ...)
% DAG is the adjacency matrix for a directed acyclic graph.
% The nodes are assumed to be in topological order. Use TOPOLOGICAL_SORT if necessary.
% For decision nodes, the parents must explicitely include all nodes
% on which it can depends, in contrast to the implicit no-forgetting assumption of influence diagrams.
% (For details, see "Representing and solving decision problems with limited information",
%    Lauritzen and Nilsson, Management Science, 2001.)
%
% node_sizes(i) is the number of values node i can take on,
%   or the length of node i if i is a continuous-valued vector.
% node_sizes(i) = 1 if i is a utility node.
%
% The list below gives optional arguments [default value in brackets].
%
% chance   - the list of nodes which are random variables [1:N]
% decision - the list of nodes which are decision nodes [ [] ]
% utility  - the list of nodes which are utility nodes [ [] ]
% equiv_class - equiv_class(i)=j  means node i gets its params from CPD{j} [1:N]
%
% e.g., limid = mk_limid(dag, ns, 'chance', [1 3], 'utility', [2])

n = length(dag);

% default values for parameters
bnet.chance_nodes = 1:n;
bnet.equiv_class = 1:n;
bnet.utility_nodes = [];
bnet.decision_nodes = [];
bnet.dnodes = 1:n; % discrete

if nargin >= 3
  args = varargin;
  nargs = length(args);
  if ~isstr(args{1})
    if nargs >= 1, bnet.dnodes = args{1}; end
    if nargs >= 2, bnet.equiv_class = args{2}; end
  else
    for i=1:2:nargs
      switch args{i},
       case 'equiv_class', bnet.equiv_class = args{i+1};
       case 'chance',      bnet.chance_nodes = args{i+1};
       case 'utility',     bnet.utility_nodes = args{i+1};
       case 'decision',    bnet.decision_nodes = args{i+1};
       case 'discrete',    bnet.dnodes = args{i+1};
        otherwise,
	error(['invalid argument name ' args{i}]);
      end
    end
  end
end

bnet.limid = 1;

bnet.dag = dag;
bnet.node_sizes = node_sizes(:)';

bnet.cnodes = mysetdiff(1:n, bnet.dnodes);
% too many functions refer to cnodes to rename it to cts_nodes -
% We hope it won't be confused with chance nodes!

bnet.parents = cell(1,n);
for i=1:n
  bnet.parents{i} = parents(dag, i);
end

E = max(bnet.equiv_class);
mem = cell(1,E);
for i=1:n
  e = bnet.equiv_class(i);
  mem{e} = [mem{e} i];
end
bnet.members_of_equiv_class = mem;

bnet.CPD = cell(1, E);

% for e=1:E
%   i = bnet.members_of_equiv_class{e}(1); % pick arbitrary member
%   switch type{e}
%     case 'tabular',  bnet.CPD{e} = tabular_CPD(bnet, i);
%     case 'gaussian', bnet.CPD{e} = gaussian_CPD(bnet, i);
%     otherwise, error(['unrecognized CPD type ' type{e}]);
%   end
% end

directed = 1;
if ~acyclic(dag,directed)
  error('graph must be acyclic')
end

bnet.order = topological_sort(bnet.dag);


text:
%kvgpwl 'Create 2D Piecewise Linear Periodic Function (K1)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vgpwl.pane file
%
% Parameters:
% InputFile: trigger 'Trigger Input', optional: 'trigger input'
% Integer: r 'Number Rows ', default: 512: 'number of rows'
% Double: s 'Sampling Freq ', default: 1: 'sampling frequency'
% Integer: c 'Number Cols ', default: 512: 'number of columns'
% Double: l 'Minimum Value ', default: 0: 'signal minimum'
% Integer: b 'Number Bands ', default: 1: 'number of bands'
% Double: m 'Maximum Value ', default: 1: 'signal maximum'
% Double: xp 'X Period ', default: 1: 'signal period in X-direction'
% Double: yp 'Y Period ', default: 1: 'signal period in Y-direction'
% Double: xr 'X Rise Time ', default: 0: 'rise time in X-direction'
% Double: yr 'Y Rise Time ', default: 0: 'rise time in Y-direction '
% Double: xf 'X Fall Time ', default: 0: 'fall time in X-direction '
% Double: yf 'Y Fall Time ', default: 0: 'fall time in Y-direction '
% Double: xw 'X Pulse Width', default: 0: 'pulse width in X-direction '
% Double: yw 'Y Pulse Width ', default: 0: 'pulse width in Y-direction '
% MultiChoice: t 'Toggle: Float, Complex', default: 1: 'select data type'
% String: a 'Float Image Type', default: 'float': 'choose data type'
% String: a 'Complex Image Type', default: 'complex': 'lets user choose data type of image'
%    Choices are:
% OutputFile: o 'Output Image', required: 'resulting image'
%
% Example: o = kvgpwl(trigger, {'trigger','';'r',512;'s',1;'c',512;'l',0;'b',1;'m',1;'xp',1;'yp',1;'xr',0;'yr',0;'xf',0;'yf',0;'xw',0;'yw',0;'t',1;'a','float';'a','complex';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vgpwl - Create 2D Piecewise Linear Periodic Function  (K1)
%
%  DESCRIPTION
% creates a two dimensional piecewise linear image.  The image pixel values
% follow a square waveform. The user can specify the period, phase,
% minima and maxima, sampling frequency, rise and fall times, and image
% dimensions.
%
% NOTE: According to the Sampling Theorem, the sampling frequency should be
% at least twice that of the signal frequency to generate a signal
% properly.
% Note that the X-direction corresponds to across the image from left to
% right, and the Y-direction corresponds to down the image from top to
% bottom.
%
% The trigger input can be used to cause
% .I vgpwl
% to re-execute when used inside cantata. No image is read from this input.
%
%
%
%  EXAMPLES
% vgpwl -o plinear.xv -t float
% vgpwl -o test.xv -t complex -xr 10 -yr 10 -xw 5 -yw 7
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% can only generate VFF_TYP_FLOAT and VFF_TYP_COMPLEX type images.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kvgpwl(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvgpwl(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'trigger', '__input';'r', 512;'s', 1;'c', 512;'l', 0;'b', 1;'m', 1;'xp', 1;'yp', 1;'xr', 0;'yr', 0;'xf', 0;'yf', 0;'xw', 0;'yw', 0;'t', 1;'a', 'float';'a', 'complex';'o', '__output'};
maxval={1,2,2,2,0,2,0,2,2,1,1,1,1,1,1,0,0,0,0};
minval={1,2,2,2,0,2,0,2,2,1,1,1,1,1,1,0,0,0,0};
istoggle=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','Integer','Double','Integer','Double','Integer','Double','Double','Double','Double','Double','Double','Double','Double','Double','MultiChoice','String','String','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vgpwl"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%
% Calculation of the fractal dimension which is the slope on the
% log-log graph.
%
deriv = diff(log10(corint))./diff(log10(r));			% deriv= Vector of the appr. derivatives
r2 = r(1:(end-1));											% Forward difference approximation: deriv has one element less, and r must have the # of elements so r2

switch(dof)

    case 'newf'

        %
        % Calculation of the fractal dimension, by first calculating the
        % distances of depopulation "rd" and of saturation "rs".
        %
        rad = 0.8;				% 2rmax= linear size of the hypercube encompassing a given dataset
        ras = 7;
        v = find(r2 >= rad & r2 <= ras);					% v= Vector of the all the interevent distances that fall in the interval [rn,rs]
        lr = log10(r2(v));
        lc = log10(corint(v));

        [fd, Err] = polyfit(lr,lc,1);
        [ev, delta] = polyval(fd, log10(r), Err);

        %[fdlo,Err] = polyfit(log10(r2(v)), log10(corint(v)), 1);
        %[fd,Err] = polyfit(log10(r2(v)), log10(corint(v)), 1);
        %[deriv3, delta] = polyval(fd, log10(r2(v)), Err);

        Ha_Cax = gca;
        Hf_Cfig;
        hold on;
        Hl_gr2a = loglog(r2(v), corint(v), 'k.','Markersize',7);
        Hl_gr2b = plot(r,10.^ev,'k','Linewidth',1);
        xlabel('Interevent Distance R [km]');
        ylabel('Correlation Integral C(R)');
        set(Ha_Cax,'pos',[0.21 0.1 0.75 0.75]);

        if d == 2
            title('Correlation Integral versus the 2D Interevent Distance R');
        else
            title('Correlation Integral versus the 3D Interevent Distance R');
        end

        g = [];
        uicontrol('Units','normal','Position',[.0 .92 .15 .07],...
            'String','Scaling range', 'Callback','g = ginput(2);dof = ''newc''; dofdim')

        str1 = ['Range = ' num2str(rd,3) ' - ' num2str(rs,3) ' [km]'];
        str2 = ['D =  ' num2str(fd(1,1),3)];
        axes('pos',[0 0 1 1]); axis off; hold on;
        te1 = text(0.25, 0.8, str1 ,'Fontweight','bold');
        te2 = text(0.25, 0.75, str2, 'Fontweight', 'bold');

        %str3 = [' The scaling range is calculated as proposed by Nerenberg & Essex (1990): ' num2str(rd,3) ' - ' num2str(rs,4) ' [km] . If you wish to change the scaling range please click on the "scaling range" button'];
        %msg1 = msgbox( str3,'Fractal Dimension');

    case 'newc'

        rd = min(g(:,1)); rs = max(g(:,1));
        v = find(r2 >= rd & r2 <= rs);

        lr = log10(r2(v));
        lc = log10(corint(v));

        [fd, Err] = polyfit(lr,lc,1);
        [ev, delta] = polyval(fd, log10(r), Err);

        delete(Hl_gr2a);
        delete(Hl_gr2b);
        Hl_gr2a = loglog(r2(v), corint(v), 'k.','Markersize',7);
        Hl_gr2b = plot(r,10.^ev,'k','Linewidth',1);

        delete(te1);
        delete(te2);
        str1 = ['Range = ' num2str(rd,3) ' - ' num2str(rs,3) ' [km]'];
        str2 = ['D =  ' num2str(fd(1,1),3)];
        axes('pos',[0 0 1 1]); axis off; hold on;
        te1 = text(0.25, 0.8, str1 ,'Fontweight','bold');
        te2 = text(0.25, 0.75, str2, 'Fontweight', 'bold');

end

author: CelsoReyes


text:
function valid = is_valid_handle(hObj)
% valid = is_valid_handle(hObj) or is_valid_handle('get_new_init_key')
%   Check if a handle is valid (has the right data type and init_key matches)
%   Use is_valid_handle('get_new_init_key') to get new init_key from C++;

% a handle is a struct array with the following fields
%   (uint64) ptr      : the pointer to the C++ object
%   (double) init_key : caffe initialization key

persistent init_key;
if isempty(init_key)
  init_key = caffe_('get_init_key');
end

% is_valid_handle('get_new_init_key') to get new init_key from C++;
if ischar(hObj) && strcmp(hObj, 'get_new_init_key')
  init_key = caffe_('get_init_key');
  return
else
  % check whether data types are correct and init_key matches
  valid = isstruct(hObj) ...
    && isscalar(hObj.ptr) && isa(hObj.ptr, 'uint64') ...
    && isscalar(hObj.init_key) && isa(hObj.init_key, 'double') ...
    && hObj.init_key == init_key;
end

end

author: aimerykong

repo: deepImageAestheticsAnalysis


text:
%GIKINE Shoulder inverse kinematics of HAL-like right shoulder
%
% Computes the inverse kinematics of the right shoulder which is the
% same kinematically as a HAL object. This function is mainly useful
% for procedures which require many, many calls, so time can be saved
% by not referencing a HAL object.
%
% Copyright (C) Bryan Moutrie, 2013-2014
% Licensed under the GNU Lesser General Public License
% see full file for full statement
%
% Syntax:
%  (1) [q1, q2] = gikine(Tg, Tu)
%
% Outputs:
%  q1 : First family of solutions (mx3 matrix where m = size(Tu,3))
%  q2 : Second family of solutions (mx3 matrix where m = size(Tu,3))
%
% Inputs:
%  Tg : Transformation matrix of the shoulder frame. x, y, z point to
%        the right, above, and behind the person. Translation is
%        shoulder center of rotation
%  Tu : Transformation matrix of the upper arm frame. May be a 4x4xm
%        series of frames (or higher order, which is compressed to 3D)
%
% See also HAL HAL.gikine HAL.ikine wikine

% LICENSE STATEMENT:
%
% This file is part of pHRIWARE.
%
%% pHRIWARE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as
% published by the Free Software Foundation, either version 3 of
% the License, or (at your option) any later version.
%
% pHRIWARE 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 Lesser General Public
% License along with pHRIWARE.  If not, see <http://www.gnu.org/licenses/>.
%
% RTB LIBRARY:
%
% Copyright (C) 199q3-2014, by Peter I. Corke
% http://www.petercorke.com
% Released under the GNU Lesser General Public license,
% Modified 16/7/2014 (HAL is a subclass of SerialLink)

function [q1, q2] = gikine(Tg, Tu)

Rg = Tg(1:3,1:3);
gRu0 = [0 0 -1; 0 1 0; 1 0 0]';
u0R = (Rg*gRu0)';

u0Ru = reshape(u0R*reshape(Tu(1:3,1:3,:),3,[]),3,3,[]);

%-X,Z,Y Tait-Bryan angles

% Family 1:
q1(:,1) = squeeze(atan2(-u0Ru(3,2,:),u0Ru(2,2,:)));
q1(:,2) = squeeze(asin(-u0Ru(1,2,:)));
q1(:,3) = squeeze(atan2(-u0Ru(1,3,:),u0Ru(1,1,:)));

% Make -pi values pi so that they are deemed within joint limits
fix = q1 <= (-pi+1e-4);
q1(fix) = pi;

% Family 2:
q2(:,1) = squeeze(atan2(u0Ru(3,2,:),-u0Ru(2,2,:)));
q2(:,2) = pi-q1(:,2);
q2(:,3) = squeeze(atan2(u0Ru(1,3,:),-u0Ru(1,1,:)));
end

author: Allopart

repo: rbpf-gmapping


text:
function y = sturm(X,BC,F,G,R)
% STURM  Solve the Sturm-Liouville equation:
%        d( F*dY/dX )/dX - G*Y = R using linear finite elements.
%
%        Y = STURM(X,BC,'F','G','R')
%
%        INPUT:
%        X  - a one-dimensional grid-point array of length N.
%        BC - is a 2 by 3 matrix [A1, B1, C1 ; An, Bn, Cn]
%             specifying the boundary conditions, which are of
%             the form: A1*Y'(1) + B1*Y(1) = C1 and
%             An*Y'(n) + Bn*Y(n) = Cn. Dirichlet boundary
%             conditions can be applied by setting A1 and/or An
%             to zero.
%        F(X), G(X) and R(X) are user-supplied functions. These
%        must be provided as separate M-files (F.M, G.M and R.M),
%        and must return an array of same shape and length as X.
%
%        The result can be displayed by: plot(X,Y).
%
% Alex Pletzer: pletzer@pppl.gov (Aug. 97/July 99).
%

[n1,n2] = size(X);
x = X(:);
n=length(x);    % # of f e
n1 = n-1;
x = sort(x);
xmin = x(1); xmax = x(n);
dx = diff(x); x = x(1:n1);

% boundary conditions

[nbc1,nbc2] = size(BC);
if ([nbc1 nbc2] ~= [2 3]),
        disp('Error calling STURM: Wrong format for BC specification')
        disp('the second argument must be a 2 times 3 matrix')
        disp('Homgeneous Dirichlet boundary conditions assumed')
        BC = [0 1 0;0 1 0]
end

if (BC(1,:) == [0 0 0]) | (BC(2,:) == [0 0 0]),
        disp('Error calling STURM: Null BC specification')
        BC
        disp('Homgeneous Dirichlet boundary conditions assumed')
        BC = [0 1 0;0 1 0]
end

a1 = BC(1,1); b1 = BC(1,2); c1 = BC(1,3);
an = BC(2,1); bn = BC(2,2); cn = BC(2,3);

% 3pt gauss' grid

abcis = [-0.77459666924148,0.0,+0.77459666924148];
abcis = (abcis+1)/2;
weigh = [0.5555555555556,0.88888888888889,0.55555555555556];
weigh = weigh/2;

x1 = x + dx*abcis(1);
x2 = x + dx*abcis(2);
x3 = x + dx*abcis(3);

xg = reshape( [conj(x1');conj(x2');conj(x3')],1,3*n1);

% f, g  and s functions

fg = feval(F,xg);
gg = feval(G,xg);
sg = feval(R,xg);

% construct matrices

fg = conj(reshape(fg,3,n1)');
gg = conj(reshape(gg,3,n1)');
sg = conj(reshape(sg,3,n1)');
g1 = gg(:,1); g2 = gg(:,2); g3 = gg(:,3);

% "kinetic"

k1 = (fg*weigh')./dx;
k1 = [0;k1]; k2 = k1(2:n);
k2 = [k2;0];
kd = k2 + k1;   % diagonal
ko = - k1(2:n); % diagonal-1
k = sparse( diag(kd) + diag(ko,-1) + diag(ko,+1) );

% "potential"

v1 = ( gg*conj(weigh.*abcis.^2)'         ).*dx; % / * /

v2 = ( gg*conj(weigh.*abcis.*(1-abcis))' ).*dx; % / * \

v3 = ( gg*conj(weigh.*(1-abcis).^2)'     ).*dx; % \ * \

v1 = [0;v1]; v3 = [v3;0];

vd = v1 + v3;
vo = v2;
v = sparse( diag(vd) + diag(vo,-1) + diag(vo,+1) );

% total

w = k + v;

% source

s1 = - ( sg*conj(weigh.*abcis)'         ).*dx; % /
s2 = - ( sg*conj(weigh.*(1-abcis))'     ).*dx; % \
s1 = [0;s1]; s2 = [s2;0];
s = s2 + s1;

% boundary conditions

        % left endpoint
        if a1 == 0,
                disp('Dirichlet BC at left endpoint:')
                disp(['y(1) = ',num2str(c1/b1)])
                w(1,:) = zeros(1,n); w(1,1)=1; s(1)=c1/b1;
        else
                disp('General BC at left endpoint:')
                disp(['y_x(1) + ',num2str(b1/a1),' y(1) = ',num2str(c1/a1)])
                w(1,1) = w(1,1) - b1/a1; s(1) = s(1) - c1/a1;
        end

        % right endpoint
        if an == 0,
                disp('Dirichlet BC at right endpoint:')
                disp(['y(n) = ',num2str(cn/bn)])
                w(n,:) = zeros(1,n); w(n,n)=1; s(n)=cn/bn;
        else
                disp('General BC at right endpoint:')
                disp(['y_x(n) + ',num2str(bn/an),' y(n) = ',num2str(cn/an)])
                w(n,n) = w(n,n) + bn/an; s(n) = s(n) + cn/an;
        end

% solve linear system

y=w\s;

[m1,m2] = size(y);
if (m1 == n2) & (m2 == n1), y = conj(y'); end


text:
%APPCR1 Character recognition.

% Mark Beale, 12-15-93
% Copyright 1992-2005 The MathWorks, Inc.
% $Revision: 1.15.2.1 $  $Date: 2005/11/15 01:14:56 $

clf;
figure(gcf)

echo on

%    NEWFF   - Inititializes feed-forward networks.
%    TRAINGDX - Trains a feed-forward network with faster backpropagation.
%    SIM   - Simulates feed-forward networks.

%    CHARACTER RECOGNITION:

%    Using the above functions a feed-forward network is trained
%    to recognize character bit maps, in the presence of noise.

pause % Strike any key to continue...

%    DEFINING THE MODEL PROBLEM
%    ==========================

%    The script file PRPROB defines a matrix ALPHABET
%    which contains the bit maps of the 26 letters of the
%    alphabet.

%    This file also defines target vectors TARGETS for
%    each letter.  Each target vector has 26 elements with
%    all zeros, except for a single 1.  A has a 1 in the
%    first element, B in the second, etc.

[alphabet,targets] = prprob;
[R,Q] = size(alphabet);
[S2,Q] = size(targets);

pause % Strike any key to define the network...

%    DEFINING THE NETWORK
%    ====================

%    The character recognition network will have 25 TANSIG
%    neurons in its hidden layer.

S1 = 10;
net = newff(minmax(alphabet),[S1 S2],{'logsig' 'logsig'},'traingdx');
net.LW{2,1} = net.LW{2,1}*0.01;
net.b{2} = net.b{2}*0.01;

pause % Strike any key to train the network...

%    TRAINING THE NETWORK WITHOUT NOISE
%    ==================================

net.performFcn = 'sse';        % Sum-Squared Error performance function
net.trainParam.goal = 0.1;     % Sum-squared error goal.
net.trainParam.show = 20;      % Frequency of progress displays (in epochs).
net.trainParam.epochs = 5000;  % Maximum number of epochs to train.
net.trainParam.mc = 0.95;      % Momentum constant.

%    Training begins...please wait...

P = alphabet;
T = targets;

[net,tr] = train(net,P,T);

%    ...and finally finishes.

pause % Strike any key to train the network with noise...

%    TRAINING THE NETWORK WITH NOISE
%    ===============================

%    A copy of the network will now be made.  This copy will
%    be trained with noisy examples of letters of the alphabet.

netn = net;

netn.trainParam.goal = 0.6;    % Mean-squared error goal.
netn.trainParam.epochs = 300;  % Maximum number of epochs to train.

%    The network will be trained on 10 sets of noisy data.

pause % Strike any key to begin training...

%    Training begins...please wait...

T = [targets targets targets targets];
for pass = 1:10
  fprintf('Pass = %.0f\n',pass);
  P = [alphabet, alphabet, ...
      (alphabet + randn(R,Q)*0.1), ...
      (alphabet + randn(R,Q)*0.2)];

  [netn,tr] = train(netn,P,T);
  echo off
end
echo on

%    ...and finally finishes.

pause % Strike any key to finish training the network...

%    TRAINING THE SECOND NETWORK WITHOUT NOISE
%    =========================================

%    The second network is now retrained without noise to
%    insure that it correctly categorizes non-noizy letters.

netn.trainParam.goal = 0.1;    % Mean-squared error goal.
netn.trainParam.epochs = 500;  % Maximum number of epochs to train.
net.trainParam.show = 5;       % Frequency of progress displays (in epochs).

%    Training begins...please wait...

P = alphabet;
T = targets;

[netn,tr] = train(netn,P,T);

%    ...and finally finishes.

pause % Strike any key to test the networks...

%    TRAINING THE NETWORK
%    ====================

% SET TESTING PARAMETERS
noise_range = 0:.05:.5;
max_test = 100;
network1 = [];
network2 = [];
T = targets;

% PERFORM THE TEST
for noiselevel = noise_range
  fprintf('Testing networks with noise level of %.2f.\n',noiselevel);
  errors1 = 0;
  errors2 = 0;

  for i=1:max_test
    P = alphabet + randn(35,26)*noiselevel;

    % TEST NETWORK 1
    A = sim(net,P);
    AA = compet(A);
    errors1 = errors1 + sum(sum(abs(AA-T)))/2;

    % TEST NETWORK 2
    An = sim(netn,P);
    AAn = compet(An);
    errors2 = errors2 + sum(sum(abs(AAn-T)))/2;
    echo off
  end

  % AVERAGE ERRORS FOR 100 SETS OF 26 TARGET VECTORS.
  network1 = [network1 errors1/26/100];
  network2 = [network2 errors2/26/100];
end
echo on

pause % Strike any key to display the test results...

%    DISPLAY RESULTS
%    ===============

%    Here is a plot showing the percentage of errors for
%    the two networks for varying levels of noise.

clf
plot(noise_range,network1*100,'--',noise_range,network2*100);
title('Percentage of Recognition Errors');
xlabel('Noise Level');
ylabel('Network 1 - -   Network 2 ---');

%    Network 1, trained without noise, has more errors due
%    to noise than does Network 2, which was trained with noise.

echo off
disp('End of APPCR1')

author: Eurus-Holmes


text:
% Author: Jai Juneja (adapted from SLAM course by Joan Sola)
% Date: 12/02/2013
%
% For a point p_g in the global frame, determine the range-bearing
% measurement y in the robot's local frame. This first requires a
% transformation of the point p_g to the robot's local frame to give p_r.
% This is then converted to a range and bearing vector with getMeasurement.
%
% Inputs:
%   r = [x y alpha]'    :   Robot frame
%   p_g = [pg_x pg_y]'  :   Point in global frame
%
% Outputs:
%   y = [d a]           :   Range-bearing sensor measurement
%   Optional:
%   Y_r                 :   Jacobian of y wrt. r
%   Y_pg                :   Jacobian of y wrt. p_g

function [y, Y_r, Y_pg] = scanPoint(r, p_g)

    if nargout == 1
        p_r = transToLocal(r, p_g);
        y = getMeasurement(p_r);    % Obtain range-bearing measurement y of p_r
    else
        % Compute Jacobians:
        [p_r, PR_r, PR_pg] = transToLocal(r, p_g);
        [y, Y_pr] = getMeasurement(p_r);

        % From the chain rule, we deduce:
        Y_r = Y_pr * PR_r;
        Y_pg = Y_pr * PR_pg;
    end
end

author: jaijuneja

repo: ekf-slam-matlab


text:
% Master m-file: seis_var.m
%
% Masterscript for the study of seismicity rate changes and variations
%
% Author: J. Woessner, woessner@seismo.ifg.ethz.ch
% last update: 26.06.02
%
basfig=figure_w_normalized_uicontrolunits('tag','welcomefig','Units','normalized','Position',[0.05 0.5 0.4 0.4],'Name','Detecting seismicity variations',...
    'NumberTitle','off');

% Buttons
seis_shift=uicontrol('Style', 'pushbutton', 'String', 'Magnitude characteristics','Units','normalized',...
    'Position',[0.1 0.10 0.3 0.06],'Callback','seisshift');
timestretch=uicontrol('Style', 'pushbutton', 'String', 'Time stretch a catalog','Units','normalized',...
    'Position',[0.1 0.80 0.3 0.06],'Callback','ftimestretch');
exit=uicontrol('Style', 'pushbutton', 'String', 'Exit','Units','normalized',...
    'Position',[0.7 0.1 0.15 0.06],'Callback','delete(findobj(''tag'',''welcomefig'')); close all');

% Figures
%topofig=figure_w_normalized_uicontrolunits('tag','Topography','Name','Topography plot','Nextplot','add','Units', 'normalized',...
%    'Position',[0.5 0.5 0.4 0.4],'Numbertitle','off');
%axtopofig=axes('tag','axtopo','NextPlot','add','box','on');