text:
%% Sobel Derivatives
%
% In this demo, we show how to:
%
% * Use the OpenCV function |cv.Sobel| to calculate the derivatives of an
%   image
% * Use the OpenCV function |cv.Scharr| to calculate a more accurate
%   derivative for a kernel of size 3x3
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/d2/d2c/tutorial_sobel_derivatives.html>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/ImgTrans/Sobel_Demo.cpp>
%

%% Theory
%
% NOTE: The explanation below belongs to the book *Learning OpenCV* by
% Bradski and Kaehler.
%
% We have previously seen applicative examples of convolutions. One of the
% most important convolutions is the computation of derivatives in an image
% (or an approximation to them). Why may be important the calculus of the
% derivatives in an image? Let's imagine we want to detect the _edges_ present
% in the image. For instance:
%
% <<https://docs.opencv.org/3.2.0/Sobel_Derivatives_Tutorial_Theory_0.jpg>>
%
% You can easily notice that in an _edge_, the pixel intensity _changes_ in a
% notorious way. A good way to express _changes_ is by using _derivatives_. A
% high change in gradient indicates a major change in the image.
%
% To be more graphical, let's assume we have a 1D-image. An edge is shown by
% the "jump" in intensity in the plot below:
%
% <<https://docs.opencv.org/3.2.0/Sobel_Derivatives_Tutorial_Theory_Intensity_Function.jpg>>
%
% The edge "jump" can be seen more easily if we take the first derivative
% (actually, here appears as a maximum)
%
% <<https://docs.opencv.org/3.2.0/Sobel_Derivatives_Tutorial_Theory_dIntensity_Function.jpg>>
%
% So, from the explanation above, we can deduce that a method to detect edges
% in an image can be performed by locating pixel locations where the gradient
% is higher than its neighbors (or to generalize, higher than a threshold).
%
%% Sobel Operator
%
% The Sobel Operator is a discrete differentiation operator. It computes an
% approximation of the gradient of an image intensity function. The Sobel
% Operator combines Gaussian smoothing and differentiation.
%
% Assuming that the image to be operated is $I$, we calculate two derivatives:
%
% * *Horizontal changes*: This is computed by convolving $I$ with a kernel
%   $G_{x}$ with odd size. For example for a kernel size of 3, $G_{x}$ would
%   be computed as:
%
% $$G_{x} = \left[{\matrix{
%               -1 & 0 & +1 \cr
%               -2 & 0 & +2 \cr
%               -1 & 0 & +1
%           }}\right] * I$$
%
% * *Vertical changes*: This is computed by convolving $I$ with a kernel
%   $G_{y}$ with odd size. For example for a kernel size of 3, $G_{y}$ would
%   be computed as:
%
% $$G_{y} = \left[{\matrix{
%               -1 & -2 & -1 \cr
%                0 &  0 &  0 \cr
%               +1 & +2 & +1
%           }}\right] * I$$
%
% At each point of the image we calculate an approximation of the _gradient_
% in that point by combining both results above:
%
% $$G = \sqrt{ G_{x}^{2} + G_{y}^{2} }$$
%
% Although sometimes the following simpler equation is used:
%
% $$G = |G_{x}| + |G_{y}|$$
%
% Note: When the size of the kernel is |3|, the Sobel kernel shown above may
% produce noticeable inaccuracies (after all, Sobel is only an approximation
% of the derivative). OpenCV addresses this inaccuracy for kernels of size 3
% by using the |cv.Scharr| function. This is as fast but more accurate than
% the standard Sobel function. It implements the following kernels:
%
% $$G_{x} = \left[{\matrix{
%                -3 & 0 &  +3 \cr
%               -10 & 0 & +10 \cr
%                -3 & 0 &  +3
%           }}\right]$$
%
% $$G_{y} = \left[{\matrix{
%               -3 & -10 & -3 \cr
%                0 &   0 &  0 \cr
%               +3 & +10 & +3
%           }}\right]$$
%
% Note: You can check out more information of this function in the OpenCV
% reference (|cv.Scharr|). Also, in the sample code below, you will notice
% that above the code for |cv.Sobel| function there is also code for the
% |cv.Scharr| function commented. Enabling it should give you an idea of how
% this function works.
%

%% Code
%
% The program below applies the _Sobel Operator_ and generates as output an
% image with the detected _edges_ bright on a darker background.
%

%%
% load source image
src = cv.imread(fullfile(mexopencv.root(),'test','lena.jpg'), 'Color',true);

%%
% apply a Gaussian blur to reduce the noise
src = cv.GaussianBlur(src, 'KSize',[3 3]);

%%
% convert filtered image to grayscale
gray = cv.cvtColor(src, 'RGB2GRAY');

%%
% calculate the derivatives in x and y directions
% (input is 8-bit, we set the output image depth to 16-bit to avoid overflow)
if true
    gradx = cv.Sobel(gray, 'XOrder',1, 'YOrder',0, 'DDepth','int16');
    grady = cv.Sobel(gray, 'XOrder',0, 'YOrder',1, 'DDepth','int16');
else
    gradx = cv.Scharr(gray, 'XOrder',1, 'YOrder',0, 'DDepth','int16');
    grady = cv.Scharr(gray, 'XOrder',0, 'YOrder',1, 'DDepth','int16');
end

%%
% take absolute value and convert our partial results back to 8-bit
gradxabs = cv.convertScaleAbs(gradx);
gradyabs = cv.convertScaleAbs(grady);

%%
% approximate the gradient by adding both directional gradients
% (this is not an exact calculation, but it is good for our purposes)
grad = cv.addWeighted(gradxabs,0.5, gradyabs,0.5, 0.0);

%%
% show result
imshow(grad)
title('Sobel Demo - Simple Edge Detector')

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
classdef LineSegment < handle
    properties
        p1_;
        p2_;
        R_;
        obs_;
        ellipsoid_;
        polyhedron_;
        local_bbox_; % Bounding Box of map
        epsilon_;
    end
    methods
        function obj = LineSegment(p1, p2)
            obj.p1_ = p1;
            obj.p2_ = p2;
            obj.epsilon_ = 1e-10;
        end

        function dilate(obj, radius)
            obj.find_ellipsoid(radius);
            obj.find_polyhedron();
            obj.add_local_bbox(obj.polyhedron_);
        end

        function set_local_bbox(obj, local_bbox)
            obj.local_bbox_ = local_bbox;
        end

        function set_obs(obj, obs)
            Vs = Polyhedron();
            obj.add_local_bbox(Vs);
            obj.obs_ = Vs.points_inside(obs);
        end

        function add_local_bbox(obj, Vs)
            r = norm(obj.p2_ - obj.p1_)/2;
            dir = (obj.p2_ - obj.p1_)/norm(obj.p2_ - obj.p1_);
            dir_h = [dir(2) -dir(1) 0];
            if(norm(dir_h) == 0)
                    dir_h = [-1 0 0];
            end
            dir_h = dir_h/norm(dir_h);

            % along x
            pp1 = obj.p1_ + dir_h*(r);
            pp2 = obj.p1_ - dir_h*(r);
            Vs.add(Hyperplane(pp1, dir_h));
            Vs.add(Hyperplane(pp2, -dir_h));

            % along y
            pp3 = obj.p2_ + dir*(r);
            pp4 = obj.p1_ - dir*(r);
            Vs.add(Hyperplane(pp3, dir));
            Vs.add(Hyperplane(pp4, -dir));

            % along z,
            dir_v = [0 0 0];
            dir_v(1) =  dir(2) * dir_h(3) - dir(3) * dir_h(2);
            dir_v(2) =  dir(3) * dir_h(1) - dir(1) * dir_h(3);
            dir_v(3) =  dir(1) * dir_h(2) - dir(2) * dir_h(1);
            pp5 = obj.p1_ + dir_v*(r);
            pp6 = obj.p1_ - dir_v*(r);
            Vs.add(Hyperplane(pp5, dir_v));
            Vs.add(Hyperplane(pp6, -dir_v));
        end

        % 3D
        function find_ellipsoid(obj, offset_x)
            f = norm(obj.p1_ - obj.p2_)/2;
            C = [f 0 0; 0 f 0; 0 0 f];
            % f is the radius of the ellipse, the initialized ellipse is a circle
            % C =  [f  0  0]
            %      [0  f  0]
            %      [0  0  f]

            % h = (R^T)*x, (1 0 0)^T = R^T*(p2_ - p1_)  => (p2_ - p1_) = R*(1 0 0)^T
            Ri = rotationMatrix([1 0 0]', (obj.p2_ - obj.p1_)');
            C = Ri * C * (Ri');
            axes = [f f f];
            obj.R_ = Ri;

            E = Ellipsoid(C, (obj.p1_ + obj.p2_)/2 );

            obs = E.points_inside(obj.obs_);
            obs_inside = obs;

            % decide short axes-1
            while (~isempty(obs_inside))
                cp = E.closest_point(obs_inside);
                rcp = (Ri')*(cp - E.d_)';

                % Generate a new ellipse
                if(abs(rcp(1)) < axes(1))
                    newy = sqrt((rcp(2))^2 + (rcp(3))^2);
                    axes(2) = newy / sqrt(1 - (rcp(1)/axes(1))^2);
                end

                new_C = [axes(1) 0 0; 0 axes(2) 0; 0 0 axes(2)];
                E.setC(Ri * new_C * (Ri'));

                % Delete all the points outside the new ellipse
                obs_new = [];
                [len, ~] = size(obs_inside);
                for i = 1 : len
                % constexpr decimal_t epsilon_ = 1e-10;
                    if(1 - E.dist(obs_inside(i,:)) > obj.epsilon_ )
                        obs_new = [obs_new; obs_inside(i,:)];
                    end
                end
                obs_inside = obs_new;
            end

            % decide short axes-2
            C = [axes(1) 0 0; 0 axes(2) 0; 0 0 axes(3)];
            E.setC(Ri * C * (Ri'));

            obs_inside = E.points_inside(obs);
            while (~isempty(obs_inside))
                cp = E.closest_point(obs_inside);
                rcp = (Ri')*(cp - E.d_)';
                dd = sqrt( 1 - (rcp(1)/axes(1))^2 - (rcp(2)/axes(2))^2 );
                if(dd > obj.epsilon_)
                    axes(3) = abs(rcp(3)) / dd;
                end

                new_C = [axes(1) 0 0; 0 axes(2) 0; 0 0 axes(3)];
                E.setC(Ri * new_C * (Ri'));

                % Delete all the points outside the new ellipse
                obs_new = [];
                [len, ~] = size(obs_inside);
                for i = 1 : len
                % constexpr decimal_t epsilon_ = 1e-10;
                    if(1 - E.dist(obs_inside(i,:)) > obj.epsilon_ )
                        obs_new = [obs_new; obs_inside(i,:)];
                    end
                end
                obs_inside = obs_new;
            end

            E.axes_ = axes;
            obj.ellipsoid_ = E;

        end

        function find_polyhedron(obj)
            Vs = Polyhedron();
            obs_remain = obj.obs_;
            while(length(obs_remain))
                plane = obj.ellipsoid_.closest_hyperplane(obs_remain);
                Vs.add(plane);
                obs_tmp = [];
                [len, ~] = size(obs_remain);
                for i = 1 : len
                    p = obs_remain(i,:);
                    if(plane.signed_dist(p) < 0)
                        obs_tmp = [obs_tmp; p];
                    end
                end
                obs_remain = obs_tmp;
            end
            obj.polyhedron_ = Vs;
        end
    end
end

author: LenaShengzhen

repo: AerialRobotics


text:
function indx = r8row_sort_heap_index_a ( m, n, a )

%*****************************************************************************80
%
%% R8ROW_SORT_HEAP_INDEX_A does an indexed heap ascending sort of an R8ROW.
%
%  Discussion:
%
%    An R8ROW is an M by N array of R8's, regarded as an array of M rows,
%    each of length N.
%
%    The sorting is not actually carried out.  Rather an index array is
%    created which defines the sorting.  This array may be used to sort
%    or index the array, or to sort or index related arrays keyed on the
%    original array.
%
%    A(I1,*) < A(I1,*) if the first nonzero entry of A(I1,*)-A(I2,*)
%    is negative.
%
%    Once the index array is computed, the sorting can be carried out
%    "implicitly:
%
%      A(INDX(1:M),1:N) is sorted.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in each column of A.
%
%    Input, integer N, the number of columns in A.
%
%    Input, real A(M,N), the array.
%
%    Output, integer INDX(M), the sort index.  The I-th element of the sorted
%    array is row INDX(I).
%
  indx(1:m) = 1 : m;

  if ( m <= 1 )
    return
  end

  l = floor ( m / 2 ) + 1;
  ir = m;

  while ( 1 )

    if ( 1 < l )

      l = l - 1;
      indxt = indx(l);
      row(1:n) = a(indxt,1:n);

    else

      indxt = indx(ir);
      row(1:n) = a(indxt,1:n);
      indx(ir) = indx(1);
      ir = ir - 1;

      if ( ir == 1 )
        indx(1) = indxt;
        break
      end

    end

    i = l;
    j = l + l;

    while ( j <= ir )

      if ( j < ir )

        if ( r8row_compare ( m, n, a, indx(j), indx(j+1) ) < 0 )
          j = j + 1;
        end

      end

      if ( r8vec_compare ( n, row, a(indx(j),1:n) ) < 0 )
        indx(i) = indx(j);
        i = j;
        j = j + j;
      else
        j = ir + 1;
      end

    end

    indx(i) = indxt;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function rez = set_cutoff(rez)
% after everything else is done, this function takes spike trains and cuts off
% any noise they might have picked up at low amplitude values
% We look for bimodality in the amplitude plot, thus setting an individual threshold
% for each neuron.
% Also, this function calls "good" and "bad" clusters based on the auto-correlogram

ops = rez.ops;
dt = 1/1000; % step size for CCG binning

Nk = numel(rez.mu); % number of templates

% sort by firing rate first
rez.good = zeros(Nk, 1);
rez.est_contam_rate = ones(Nk, 1);

for j = 1:Nk
    ix = find(rez.st3(:,2)==j); % find all spikes from this neuron
    ss = rez.st3(ix,1)/ops.fs; % convert to seconds
    if numel(ss)==0
        continue; % break if there are no spikes
    end

    vexp = rez.st3(ix,4); % vexp is the relative residual variance of the spikes

    Th = ops.Th(1); % start with a high threshold

    fcontamination = 0.1; % acceptable contamination rate

    rez.est_contam_rate(j) = 1;
    while Th>=ops.Th(2)
      % continually lower the threshold, while the estimated unit contamination is low
        st = ss(vexp>Th); % take spikes above the current threshold
        if isempty(st)
            Th = Th - .5; % if there are no spikes, we need to keep lowering the threshold
            continue;
        end
        [K, Qi, Q00, Q01, rir] = ccg(st, st, 500, dt); % % compute the auto-correlogram with 500 bins at 1ms bins
        Q = min(Qi/(max(Q00, Q01))); % this is a measure of refractoriness
        R = min(rir); % this is a second measure of refractoriness (kicks in for very low firing rates)
        if Q>fcontamination || R>.05 % if the unit is already contaminated, we break, and use the next higher threshold
           break;
        else
            if Th==ops.Th(1) && Q<.05
              % only on the first iteration, we consider if the unit starts well isolated
              % if it does, then we put much stricter criteria for isolation
              % to make sure we don't settle for a relatively high contamination unit
                fcontamination = min(.05, max(.01, Q*2));
%                 fcontamination = .05;

                % if the unit starts out contaminated, we will settle with the higher contamination rate
            end
            rez.good(j) = 1; % this unit is good, because we will stop lowering the threshold when it becomes bad
            Th = Th - .5; % try the next lower threshold
        end
    end
    Th = Th + .5; % we exited the loop because the contamination was too high. We revert to the higher threshold
    st = ss(vexp>Th); % take spikes above the current threshold
    if isempty(st)
        rez.est_contam_rate(j) = 1;
    else
        [K, Qi, Q00, Q01, rir] = ccg(st, st, 500, dt); % % compute the auto-correlogram with 500 bins at 1ms bins
        Q = min(Qi/(max(Q00, Q01))); % this is a measure of refractoriness
        rez.est_contam_rate(j) = Q; % this score will be displayed in Phy
    end

    rez.Ths(j) = Th; % store the threshold for potential debugging
    rez.st3(ix(vexp<=Th), 2) = 0; % any spikes below the threshold get discarded into a 0-th cluster.
end

% we sometimes get NaNs, why? replace with full contamination
rez.est_contam_rate(isnan(rez.est_contam_rate)) = 1;

% remove spikes from the 0th cluster
rez = remove_spikes(rez,rez.st3(:,2)==0,'below_cutoff');

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function [geom_data, valid_ids] = Read_geom_files(users, model_param_data_dir)

    geom_data = [];
    valid_ids = [];

    load('../../pca_generation/pdm_68_aligned_wild.mat');

    for i=1:numel(users)

        geom_file = [model_param_data_dir, '/' users{i} '.csv'];

        if(i == 1)
            tab = readtable(geom_file);
            column_names = tab.Properties.VariableNames;
            valid_ind = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'confidence'));
            shape_inds = cellfun(@(x) ~isempty(x) && x==1, strfind(column_names, 'p_'));
        end

        res = dlmread(geom_file, ',', 1, 0);

        % Check the confidence of detection
        valid = res(:, valid_ind) > 0.7;
        res = res(:, shape_inds);

        % Do not consider global parameters
        res = res(:, 7:end);

        actual_locs = res * V';
        res = cat(2, actual_locs, res);

        valid_ids = cat(1, valid_ids, valid);

        geom_data = cat(1, geom_data, res);

    end
end

author: TadasBaltrusaitis

repo: OpenFace


text:
%% Calculate errors between the corrected image and the ground truth image.
%
% Copyright (c) 2018-present, Mahmoud Afifi
% York University, Canada
% mafifi@eecs.yorku.ca | m.3afifi@gmail.com
%
% This source code is licensed under the license found in the
% LICENSE file in the root directory of this source tree.
% All rights reserved.
%
% Please cite the following work if this program is used:
% Mahmoud Afifi, Brian Price, Scott Cohen, and Michael S. Brown,
% "When color constancy goes wrong: Correcting improperly white-balanced
% images", CVPR 2019.
%
% Input:
%   -corrected: corrected image.
%   -gt: ground-truth image.
%   -color_chart_area: If there is a color chart in the image, that is
%   masked out from both images, this variable represents the number of
%   pixels of the color chart.
%   -opt: determines the required error metric(s) to be reported.
%         Options:
%           opt = 1 delta E 2000 (default).
%           opt = 2 delta E 2000 and mean squared error (MSE)
%           opt = 3 delta E 2000, MSE, and mean angular eror (MAE)
%           opt = 4 delta E 2000, MSE, MAE, and delta E 76
%
% Output:
%   -varargout: a cell contains error between corrected and gt images.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%

function varargout = evaluate_cc(corrected, gt, color_chart_area, opt)

if nargin == 3
    opt = 1;
end
switch opt
    case 1
        varargout{1} = calc_deltaE2000(corrected,gt,color_chart_area);
    case 2
        varargout{1} = calc_deltaE2000(corrected,gt,color_chart_area);
        varargout{2} =calc_mse(corrected,gt,color_chart_area);
    case 3
        varargout{1} = calc_deltaE2000(corrected,gt,color_chart_area);
        varargout{2} =calc_mse(corrected,gt,color_chart_area);
        varargout{3} =calc_mae(reshape(corrected,[],3),reshape(gt,[],3),...
            color_chart_area);
    case 4
        varargout{1} = calc_deltaE2000(corrected,gt,color_chart_area);
        varargout{2} =calc_mse(corrected,gt,color_chart_area);
        varargout{3} =calc_mae(reshape(corrected,[],3),reshape(gt,[],3),...
            color_chart_area);
        varargout{4} = calc_deltaE(corrected,gt,color_chart_area);
    otherwise
        error('Error in evaluate_cc function');
end
end

author: mahmoudnafifi


text:
  function [xo, yo] = plot_ellipse(cx, cy, rx, ry, theta, varargin)
%|function [xo, yo] = plot_ellipse(cx, cy, rx, ry, theta, [options])
%| plot an ellipse
%| option
%|	'n'		# of points for half of ellipse.  default: 301
%|	'c'		line type.  default: 'b-'
%|	'hold'	1|0	if 1, then add to current plot.  default: 0

if nargin == 1 && streq(cx, 'test'), plot_ellipse_test, return, end
if nargin < 5, ir_usage, end

arg.n = 301;
arg.c = 'b-';
arg.hold = false;
arg = vararg_pair(arg, varargin);

xe = linspace(-rx, rx, arg.n)';
yp = ry * sqrt(1 - (xe/rx).^2);
ym = -yp;
xp = cx + (cos(theta) * xe - sin(theta) * yp);
yp = cy + (sin(theta) * xe + cos(theta) * yp);
xm = cx + (cos(theta) * xe - sin(theta) * ym);
ym = cy + (sin(theta) * xe + cos(theta) * ym);

xo = [xp; flipud(xm)];
yo = [yp; flipud(ym)];

if ~nargout
	if im
		if arg.hold, hold on, end
		plot(xo, yo, arg.c)
		if arg.hold, hold off, end
	end
	clear xo yo
end

function plot_ellipse_test
plot_ellipse(5, 10, 10, 20, 1)
[xo yo] = plot_ellipse(5, 10, 10, 20, 1);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function vw = viewSetUI(vw,param,val,varargin)
%Organize methods for setting view parameters.
%
% This function is wrapped by viewSet. It should not be called by anything
% else other than viewSet.
%
% This function retrieves information from the view that relates to a
% specific component of the application.
%
% We assume that input comes to us already fixed and does not need to be
% formatted again.

if notDefined('vw'),  error('No view defined.'); end
if notDefined('param'), error('No parameter defined'); end
if notDefined('val'),   val = []; end

mrGlobals;

switch param

    case 'initdisplaymodes'
        vw = resetDisplayModes(vw);
    case 'ui'
        vw.ui = val;
    case 'anatomymode'
        vw.ui.anatMode = val;
    case 'coherencemode'
        vw.ui.coMode = val;
    case 'correlationmode'
        vw.ui.corMode = val;
    case 'phasemode'
        vw.ui.phMode = val;
    case 'fignum'
        vw.ui.figNum = val;
    case 'windowhandle'
        vw.ui.windowHandle = val;
    case 'mainaxishandle'
        vw.ui.mainAxisHandle = val;
    case 'colorbarhandle'
        vw.ui.colorbarHandle = val;
    case 'cbarrange'
        vw.ui.cbarRange = val;
    case 'amplitudemode'
        vw.ui.ampMode = val;
    case 'uiimage'
        vw.ui.image = val;
    case 'projectedamplitudemode'
        vw.ui.projampMode = val;
    case 'mapmode'
        vw.ui.mapMode = val;
    case 'displaymode'
        vw.ui.displayMode = val;
    case 'phasecma'
        nGrays = viewGet(vw, 'curnumgrays');

        % allow transposed version (3 x n) instead of usual matlab cmap order (n x 3)
        if size(val, 2) > 3 && size(val, 1)==3, val = val'; end
        if max(val(:)) > 1, val = val ./ 255;               end
        vw.ui.phMode.cmap((nGrays+1):end,:) = val ;

    case 'locs'
        % cursor location as [axi cor sag]
        vw.loc = val;
        if checkfields(vw, 'ui', 'sliceNumFields')  % set UI fields
            for n = 1:3
                set(vw.ui.sliceNumFields(n), 'String', num2str(val(n)));
            end
        end
    case 'crosshairs'
        vw.ui.crosshairs = val;
    case 'flipud'
        % Boolean indicating whether to invert the image u/d in
        % the graphical user interface.  It is sometimes convenient to do
        % this in the Inplane view if the top of the slice corresponds to
        % the bottom of the brain.
        % Example:
        %   vw = viewSet(vw, 'flip updown', true);
        if checkfields(vw, 'ui'),   vw.ui.flipUD = val; end
    case 'zoom'
        vw.ui.zoom = val;
    otherwise
        error('Unknown view parameter %s.', param);

end %switch

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Calculate border coverage for detected fields
%
% This function calculates firing map border coverage that is further
% used in calculation of a border score.
%
%  USAGE
%   coverage = analyses.borderCoverage(fields, <options>)
%   fields          Array of structures with information about detected fields.
%   <options>       Optional list of property-value pairs (see table below)
%
%   ==============================================================================================
%    Properties    Values
%   ----------------------------------------------------------------------------------------------
%    'searchWidth'  If map is not perfect, but contains NaN values along borders, then
%                   search for border pixels can have NaNs. To mitigate this, we check
%                   searchWidth rows/columns near border and if the closest to the border pixel
%                   equals to NaN, we search for first non-NaN value in searchWidth rows-columns.
%                   This argument is optional and default value is 8.
%    'walls'        Definition of walls along which the border score is calculated. Provided by
%                   a string which contains characters that stand for walls:
%                   T - top wall (we assume the bird-eye view on the arena)
%                   R - right wall
%                   B - bottom wall
%                   L - left wall
%                   Characters are case insensitive. Default value is 'TRBL' meaning that border
%                   score is calculated along all walls. Any combination is possible, e.g.
%                   'R' to calculate along right wall, 'BL' to calculate along two walls, e.t.c.
%   ==============================================================================================
%   coverage        Border coverage, ranges from 0 to 1.
%
function coverage = borderCoverage(fields, varargin)
    coverage = 0;

    inp = inputParser;
    defaultSearchWidth = 8;
    defaultWalls = 'TRBL'; % top, right, bottom, left

    checkSearchWidth = @(x) isnumeric(x) && isscalar(x) && (x > 0);
    checkWalls = @internCheckWalls;

    addRequired(inp, 'fields');
    addParameter(inp, 'searchWidth', defaultSearchWidth, checkSearchWidth);
    addParameter(inp, 'walls', defaultWalls, checkWalls);

    inp.KeepUnmatched = true;
    parse(inp, fields, varargin{:});

    % get parsed results
    walls = inp.Results.walls;
    searchWidth = inp.Results.searchWidth;

    % find out what walls are present
    dict = {'B', 'L', 'R', 'T'}; % already sorted
    u = unique(walls);
    t = ['^' sprintf('%c{0,%d}', [u; histc(walls, u)]) '$'];
    wallsIdx = find(~cellfun('isempty', regexpi(dict, t))); % indices in dict of present walls

    for i=1:length(fields)
        curField = fields(i).map;

        % coverage for left wall
        if any(ismember(wallsIdx, 2))
            aux_map = curField(:, 1:searchWidth);
            [covered, norm] = wall_field(aux_map);
            if (covered/norm > coverage)
                coverage = covered/norm;
            end
        end

        % coverage for right wall
        if any(ismember(wallsIdx, 3))
            aux_map = curField(:, end:-1:end+1-searchWidth);
            [covered, norm] = wall_field(aux_map);
            if (covered/norm > coverage)
                coverage=covered/norm;
            end
        end

        % coverage for bottom wall, since we are dealing with data that came from a camera
        % 'bottom' is actually at the top of the matrix curField.
        if (any(ismember(wallsIdx, 1)))
            aux_map = curField(1:searchWidth, :)';
            [covered, norm] = wall_field(aux_map);
            if (covered/norm > coverage)
                coverage=covered/norm;
            end
        end

        % coverage for top wall
        if (any(ismember(wallsIdx, 4)))
            aux_map = curField(end:-1:end+1-searchWidth, :)';
            [covered, norm] = wall_field(aux_map);
            if (covered/norm > coverage)
                coverage = covered/norm;
            end
        end
    end
end

% 'covered' pixels will have distance to border 0.
% Essentially we need to calculate number of elements,
% that equal to zero and take NaNs into account.
function [covered, norm] = wall_field(map)
    ly = size(map, 1);

    D = bwdist(map);
    nanIndices = find(isnan(map(:, 1)));
    numNans = length(nanIndices);
    for i = 1:numNans
        testRow = nanIndices(i);
        nonNan = find(~isnan(map(testRow, :)), 1, 'first');
        if ~isempty(nonNan)
            numNans = numNans - 1;
            D(testRow, 1) = D(testRow, nonNan);
        end
    end

    norm = ly - numNans;
    covered = nansum(D(:, 1) == 0) - numNans;
end

% check input argument that defines walls
% We need to figure out if argument contains one of the characters
% from the dictionary.
% Rise descriptive exception in case of an error.
function res = internCheckWalls(walls)
    res = true;
    if ~ischar(walls)
        error('BNT:args:notChar', 'Argument is not a string.')
    end
    if length(walls) > 4
        error('BNT:args:length', 'Argument can not exceed 4 characters (default ''TLBR'')');
    end

    dict = {'T', 'R', 'B', 'L'};

    % http://stackoverflow.com/questions/19343339/matlab-find-the-indices-of-a-cell-array-of-strings-with-characters-all-containe
    u = unique(walls);
    t = ['^' sprintf('%c{0,%d}', [u; histc(walls, u)]) '$'];
    s = cellfun(@sort, dict, 'uni', 0);
    idx = find(~cellfun('isempty', regexp(s, t)), 1);
    %res = ~isempty(idx); % if it is not empty, then we have characters from dict in walls
    if isempty(idx)
        error('BNT:args:noValidChars', 'There is no information about the walls in the argument.');
    end
end

author: kavli-ntnu


text:
%% simplexImIntersect
% Below is a demonstration of the features of the |simplexImIntersect| function

%% Syntax
% |[regionlabel]=simplexImIntersect(F1,V1,V2,voxelSize);|

%% Description
% This function uses the input simplex (or simplices) defined by F1
% (faces), V1 (vertices), and C1 (boundary labels) to label the vertices V2
% according to what simplex regions they are contained in. The labels are
% stored in the output regionLabel. A vertex inside a simplex region is
% labelled with the boundary label of the simplex region. Vertices outside
% of all simplices are labelled with NaN. Points are labelled based on the
% patch2Im function, i.e. the following steps are used:
% 1) The simplex is converted to an image where voxels intensities denote
% simplex region labels. (see patch2Im).
% 2) The vertices in V2 are converted to image coordinates in this image
% 3) The image coordinates are converted to image voxel indices
% 4) Voxel indices are used to retrieve the labels from the simplex image.
%
% The image constructed uses the optional input voxelSize (the default if
% not provided is half of the mean edge size of the input simplex). The
% completeness/accuracy of the labelling depends on the voxel size. Some
% points are labelled as 0 which means they are found inside boundary
% voxels and cannot, based on the current voxel size, be assigned as
% outside or inside a particular region. All vertices labelled as NaN or a
% value >0 are labelled correctly. However some of the vertices labelled as
% 0 may actually be inside a simplex region or outside all simplex regions,
% i.e. 0 denotes that their status is unknown given the voxel size used.
%
% See also |patch2im|

%% Examples
clear; close all; clc;

%%
% Plot settings
faceAlpha1=1;
faceAlpha2=0.3;
fontSize=15;
markerSize=35;

%%
% Create example geometries
[F1,V1,~]=geoSphere(3,1); %First simplex
[F2,V2,~]=geoSphere(3,0.75); %Second simplex
V2(:,2)=V2(:,2)+1; %Shift second sphere

%%
% Determine optional voxel size input from mean edge size
[D1]=patchEdgeLengths(F1,V1);
[D2]=patchEdgeLengths(F2,V2);
d=mean([D1(:);D2(:)]);
voxelSize=d/2;

%% Find points outside of a simplex

[regionLabel]=simplexImIntersect(F1,V1,[],V2,voxelSize);
logicOut=isnan(regionLabel);
logicFacesOut=all(logicOut(F2),2);

%%
% Visualize results

cFigure;
subplot(1,3,1); hold on;
title('Intersecting sets')
gpatch(F1,V1,'g','none',faceAlpha2);
gpatch(F2,V2,'kw','k',faceAlpha1);
axisGeom;
view(-90,0);
camlight headlight;

subplot(1,3,2); hold on;
title('Point labels')
gpatch(F1,V1,'g','none',faceAlpha2);
gpatch(F2,V2,'kw','none',faceAlpha1);
scatterV(V2,markerSize,regionLabel,'filled');
colormap gjet; icolorbar;
axisGeom;
view(-90,0);
camlight headlight;

subplot(1,3,3); hold on;
title('Cropped geometry')
gpatch(F1,V1,'g','none',faceAlpha2);
gpatch(F2(logicFacesOut,:),V2,'kw','k',faceAlpha1);
axisGeom;
view(-90,0);
camlight headlight;

drawnow;

%%
% Defining a multi boundary set

%Example surface set 1
r=2; %Sphere radius
rc=3; %Central radius
nr=15;
nc=25;
ptype='quad';
[Fs1,Vs1]=patchTorus(r,nr,rc,nc,ptype);
[Fs2,Vs2]=quadSphere(2,r,2);
Vs2(:,2)=Vs2(:,2)-5;
[Fs3,Vs3]=quadSphere(2,r/2,2);
Vs3(:,2)=Vs3(:,2)-5;
[Fs4,Vs4]=quadSphere(3,r/2,2);
Vs4(:,1)=Vs4(:,1)+2;
Vs4(:,2)=Vs4(:,2)+2;

[F1,V1,C1]=joinElementSets({Fs1,Fs2,Fs3,Fs4},{Vs1,Vs2,Vs3,Vs4});

%Example surface 2
[F2,V2]=stanford_bunny; %Bunny surface
[F2,V2]=subtri(F2,V2,1); %Refine surface
r=sqrt(sum(V2.^2,2)); %Radii from centre of bunny to "normalize" shape
V2=V2./max(r); %"Normalize"
V2=V2*11; %Scale up

%%

[regionLabel]=simplexImIntersect(F1,V1,C1,V2,voxelSize);
logicOut=isnan(regionLabel);
logicFacesOut=all(logicOut(F2),2);

%%
% Visualize results

cFigure;
subplot(1,3,1); hold on;
title('Intersecting sets')
gpatch(F1,V1,C1,'none',faceAlpha2);
gpatch(F2,V2,'kw','none',faceAlpha1);
axisGeom;
camlight headlight;

subplot(1,3,2); hold on;
title('Point labels')
gpatch(F1,V1,C1,'none',faceAlpha2);
gpatch(F2,V2,'kw','none',faceAlpha1);
scatterV(V2,markerSize,regionLabel,'filled');
colormap gjet; icolorbar;
axisGeom;
camlight headlight;

subplot(1,3,3); hold on;
title('Cropped geometry')
gpatch(F1,V1,C1,'none',faceAlpha2);
gpatch(F2(logicFacesOut,:),V2,'kw','none',faceAlpha1);
axisGeom;
camlight headlight;

drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function X = AnaMMV(Y, HforVec, AforVec, err, Hnormfac, Anormfac)

% Solves problem of the form:
% min ||AX||_2,1 s.t. ||Y-HX|| < err

% Inputs
% Y - n x r observation vectors
% Hforvec - n x N measurement operator
% Aforvec - M X N sparsifyng operator
% err - error tolerance
% Hnormfac - maximum eigenvalue of measurement operator
% Anormfac - maximum eigenvalue of sparsifying transform

% Outputs
% X - N x r input vectors to be recovered

% Converting input measurement matrix to operator
explicitH = ~(ischar(HforVec) || isa(HforVec, 'function_handle'));
if (explicitH)
    HOp = opMatrix(HforVec); clear HforVec
else
    HOp = HforVec;
end
% Converting input sparsifying matrix to operator
explicitA = ~(ischar(AforVec) || isa(AforVec, 'function_handle'));
if (explicitA)
    AOp = opMatrix(AforVec); clear HforVec
else
    AOp = AforVec;
end

r = size(Y,2); % number of observations/inputs

H = opMatWrapper(HOp, r); % wrapper for handling matrix inputs
A = opMatWrapper(AOp, r); % wrapper for handling matrix inputs

alpha = 1.05*(Hnormfac^2);
c = 1.05*(Anormfac^2);

X = H(Y,2); % Initialize X
Z = A(X,1).*0; % Initialize Z
N = size(X,1); % length of each input vector
M = size(Z,1); % length of transform domain sparse vector

maxIter = 100; % Define the maximum number of iterations
tol = 1e-4; % tolerance level
decfac = 0.5; % decrease factor for lambda

lambdaInit = decfac*max(max((abs(X)))); lambda = lambdaInit;

while lambda > lambdaInit*tol
    iter = 0;
    while iter < maxIter
        iter = iter + 1;

        W = A(X,1);
        for i = 1:M
            D(i,:) = (1/norm(W(i,:)));
        end

        B = X + (1/alpha)*H(Y-H(X,1),2);
        Z = diag(1./((alpha/lambda)*(1./D) + c))*(c*Z + A(B - A(Z,2),1));
        X = B - A(Z,2);

        if norm(Y-H(X,1),'fro') < err
            break;
        end
    end
    lambda = decfac*lambda;
end

end

author: Sable

repo: mcbench-benchmarks


text:
function result = set_joinepos(varargin)
% Join epoched signals across epochs.
% Result = set_joinepos(Set1, Set2, ...)
%
% In:
%   SetK   : The k'th data set to join.
%
% Out:
%   Result : A new data set with trials from all sets joined. The following changes are made:
%            * .data and all other time-series fields are joined across trials (3rd dimension)
%            * .epoch is joined and its .event field is updated appropriately
%            * .event is joined and its .epoch and .latency fields are updated appropriately
%            * .xmax/.trials are updated
%
% Notes:
%   This function returns a new data set with meta-data set to that of the first input set, and the
%   trials joined across all sets. No checks for meta-data consistency are done. There is a
%   heavy-duty function for merging inconsistent sets called set_merge, which can merge cats and
%   dogs. This function does not attempt to keep miscellaneous EEGLAB meta-data consistent, including:
%   setname,filename,filepath,subject,group,condition,session,comments,urevent,reject,stats,history,etc
%
%   To concatenate continuous sets across time, use set_concat.
%
% Examples:
%   % merge data sets eegA, eegB and eegC across epochs
%   eeg = set_joinepos(eegA,eegB,eegC)
%
% See also:
%   set_concat, set_merge
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-03-31
dp;

% set_joinepos_version<1.0> -- for the cache

if ~exp_beginfun('editing') return; end

declare_properties('name','JoinTrials','independent_channels',true,'independent_trials',true);

% input validation
for k=1:length(varargin)
    utl_check_fields(varargin{k},{'data','event','epoch','pnts','xmin','xmax','srate'},'input','signal');
    if ~isempty(varargin{k}.event)
        if ~isfield(varargin{k}.event,'latency')
            error('A dataset passed to set_joineos is lacking the .event.latency field.'); end
        if ~isfield(varargin{k}.event,'epoch')
            error('A dataset passed to set_joineos is lacking the .event.epoch field.'); end
        latency_numels = cellfun('prodofsize',{varargin{k}.event.latency});
        if any(latency_numels == 0)
            error('One or more of the events in the given data set have an empty .latency field, which is not permitted.'); end
        if any(latency_numels ~= 1)
            error('One or more of the events in the given data set have a .latency value that is not a scalar, which is not permitted.'); end
    end
    if ~isempty(varargin{k}.epoch)
        if ~isfield(varargin{k}.epoch,'event')
            error('A dataset passed to set_joineos is lacking the .epoch.event field.'); end
        epoch_numels = cellfun('prodofsize',{varargin{k}.event.epoch});
        if any(epoch_numels == 0)
            error('One or more of the events in the given data set have an empty .epoch field, which is not permitted.'); end
        if any(epoch_numels ~= 1)
            error('One or more of the events in the given data set have an .epoch value that is not a scalar, which is not permitted.'); end
    end
end

if ~isempty(varargin)
    result = varargin{1};
    if length(varargin) > 1
        % concatenate time series fields
        for field = utl_timeseries_fields(result)
            data = cellfun(@(x)x.(field{1}),varargin,'UniformOutput',false);
            try
                result.(field{1}) = cat(3,data{:});
            catch e
                % concatenation failed: produce a reasonable error message
                if ~isempty(data)
                    sizes1 = cellfun('size',data,1); sizes1 = sizes(sizes1 ~= 0);
                    sizes2 = cellfun('size',data,2); sizes2 = sizes(sizes2 ~= 0);
                    if length(sizes1) > 1 && ~all(sizes1==sizes1(1))
                        error('The time-series field .%s must have the same number of channels in each data set.',field{1}); end
                    if length(sizes2) > 1 && ~all(sizes2==sizes2(1))
                        error('The time-series field .%s must have the same number of tine points in each data set.',field{1}); end
                end
                size_info = hlp_tostring(cellfun(@size,data,'UniformOutput',false));
                error('Concatenation of time-series fields failed with error: %s (the data sizes were %s -- make sure that they are mutually compatible).',e.message,size_info);
            end
            if isempty(result.(field{1}))
                result.(field{1}) = []; end
        end
        % count events, epochs and samples in each set
        event_count = cellfun(@(x)length(x.event),varargin);
        epoch_count = cellfun(@(x)length(x.epoch),varargin);
        sample_count = cellfun(@(x)x.pnts,varargin).*epoch_count;
        % concatenate .event and .epoch fields
        event = cellfun(@(x)x.event,varargin,'UniformOutput',false); result.event = [event{:}];
        epoch = cellfun(@(x)x.epoch,varargin,'UniformOutput',false); result.epoch = [epoch{:}];
        % shift event latencies based on cumulative sample counts
        if ~isempty(result.event)
            [result.event.latency] = arraydeal([result.event.latency]+replicate(cumsum(sample_count)-sample_count,event_count));
            % shift event/epoch cross-references based on cumulative counts
            [result.event.epoch] = arraydeal([result.event.epoch]+replicate(cumsum(epoch_count)-epoch_count,event_count));
            [result.epoch.event] = chopdeal([result.epoch.event]+replicate(cumsum(event_count)-event_count,event_count),cellfun('length',{result.epoch.event}));
        end
        % update misc fields
        [result.nbchan,result.pnts,result.trials,extra_dims] = size(result.data); %#ok<NASGU>
        result.xmax = result.xmin + (result.pnts-1)/result.srate;
    end
else
    result = struct('setname','','filename','','filepath','','subject','','group','','condition','','session',[],'comments','','nbchan',0,...
        'trials',0,'pnts',0,'srate',1,'xmin',0,'xmax',0,'times',[],'data',[],'icaact',[],'icawinv',[],'icasphere',[],'icaweights',[], ...
        'icachansind',[],'chanlocs',[],'urchanlocs',[],'chaninfo',[],'ref',[],'event',[],'urevent',[],'eventdescription',{{}}, ...
        'epoch',[],'epochdescription',{{}},'reject',[],'stats',[],'specdata',[],'specicaact',[],'splinefile','','icasplinefile','', ...
        'dipfit',[],'history','','saved','no','etc',[]);
end

exp_endfun;

function result = replicate(values,counts)
% Replicate each element Values(k) by Count(k) times.
result = zeros(1,sum(counts));
k = 0;
for p=find(counts)
    result(k+(1:counts(p))) = values(p);
    k = k+counts(p);
end

repo: REST


text:
close all;
clear;
clc;

% [lat, lon, alt, vel, att, pos_std, vel_std, att_std, ins_status, pos_type, time] = inspvaxa2pvax('./0115/BY.txt');
[lat, lon, alt, vel, att, pos_std, vel_std, att_std, ins_status, pos_type, time] = inspvaxa2pvax('./0122/BY_INSPVAX.txt');
[neu] = pos2xyz(lat, lon, alt);

%%
sg_xyz = neu(find(pos_type=='INS_PSRSP'), :);
dg_xyz = neu(find(pos_type=='INS_PSRDIFF'), :);
float_xyz = neu(find(pos_type=='INS_RTKFLOAT'), :);
fixed_xyz = neu(find(pos_type=='INS_RTKFIXED'), :);

figure;
hold on; grid on;
legend_str = string([]);
if ~isempty(sg_xyz)
    plot(sg_xyz(:, 1), sg_xyz(:, 2), 'r.', 'LineWidth', 3);
    legend_str = [legend_str; 'Standalone'];
end
if ~isempty(dg_xyz)
    plot(dg_xyz(:, 1), dg_xyz(:, 2), 'm.', 'LineWidth', 3);
    legend_str = [legend_str; 'DGNSS'];
end
if ~isempty(float_xyz)
    plot(float_xyz(:, 1), float_xyz(:, 2), 'b.', 'LineWidth', 3);
    legend_str = [legend_str; 'RTK Float'];
end
if ~isempty(fixed_xyz)
    plot(fixed_xyz(:, 1), fixed_xyz(:, 2), 'g.', 'LineWidth', 3);
    legend_str = [legend_str; 'RTK Fixed'];
end
plot(neu(:, 1), neu(:, 2), 'k');
axis equal;
xlabel('东向距离(m)');
ylabel('北向距离(m)');
legend(legend_str);

set(gcf, 'Units', 'normalized', 'Position', [0.025, 0.05, 0.95, 0.85]);

%%
figure('name', '高度与位置方差曲线');
subplot(2,1,1);
plot(neu(:, 3), 'LineWidth', 1.5); grid on;
xlim([0 length(ins_status)]);
xlabel('时间(s)');
ylabel('高度(m)');
title('高度曲线');
subplot(2,1,2);
plot(pos_std(:, 1), 'LineWidth', 1.5); hold on; grid on;
plot(pos_std(:, 2), 'LineWidth', 1.5);
plot(pos_std(:, 3), 'LineWidth', 1.5);
xlim([0 length(ins_status)]);
xlabel('时间(s)');
ylabel('位置方差(m)');
legend('东向', '北向', '天向', 'Orientation', 'horizontal');
title('位置方差曲线');

set(gcf, 'Units', 'normalized', 'Position', [0.025, 0.05, 0.95, 0.85]);

%%
figure('name', '水平速度曲线');
subplot(2,1,1);
plot(sqrt(vel(:, 1).^2 + vel(:, 2).^2)*3.6, 'LineWidth', 1.5); hold on; grid on;
xlim([0 length(ins_status)]);
ylabel('水平速度(km/h)');
title('水平速度');
subplot(2,1,2);
plot(sqrt(vel_std(:, 1).^2 + vel_std(:, 2).^2), 'LineWidth', 1.5); hold on; grid on;
xlim([0 length(ins_status)]);
xlabel('时间(s)');
ylabel('水平速度方差(m/s)');

set(gcf, 'Units', 'normalized', 'Position', [0.025, 0.05, 0.95, 0.85]);

%%
figure('name', '东北天速度曲线');
subplot(2,1,1);
plot(vel(:, 1), 'LineWidth', 1.5); hold on; grid on;
plot(vel(:, 2), 'LineWidth', 1.5);
plot(vel(:, 3), 'LineWidth', 1.5);
xlim([0 length(ins_status)]);
ylabel('东北天速度(m/s)');
title('东北天速度曲线');
legend('东向', '北向', '天向', 'Orientation', 'horizontal');
subplot(2,1,2);
plot(vel_std(:, 1), 'LineWidth', 1.5); hold on; grid on;
plot(vel_std(:, 2), 'LineWidth', 1.5);
plot(vel_std(:, 3), 'LineWidth', 1.5);
xlim([0 length(ins_status)]);
xlabel('时间(s)');
ylabel('东北天速度方差(m/s)');
legend('东向', '北向', '天向', 'Orientation', 'horizontal');

set(gcf, 'Units', 'normalized', 'Position', [0.025, 0.05, 0.95, 0.85]);

%%
figure('name', '姿态曲线');
subplot(3,1,1);
plot(att(:, 1), 'r', 'LineWidth', 1.5); hold on; grid on;
plot(att(:, 2), 'b', 'LineWidth', 1.5);
xlim([0 length(ins_status)]);
ylabel('水平姿态(°)');
legend('俯仰', '横滚', 'Orientation', 'horizontal');
subplot(3,1,2);
plot(att(:, 3), 'm', 'LineWidth', 1.5); grid on;
xlim([0 length(ins_status)]);
ylabel('航向(°)');
legend('航向', 'Orientation', 'horizontal');
subplot(3,1,3);
plot(att_std(:, 1), 'r', 'LineWidth', 1.5); hold on; grid on;
plot(att_std(:, 2), 'b', 'LineWidth', 1.5);
plot(att_std(:, 3), 'm', 'LineWidth', 1.5);
xlim([0 length(ins_status)]);
xlabel('时间(s)');
ylabel('方差(°)');
legend('俯仰', '横滚', '航向', 'Orientation', 'horizontal');

set(gcf, 'Units', 'normalized', 'Position', [0.025, 0.05, 0.95, 0.85]);

%%
wm = webmap('World Imagery');
wmline(lat, lon, 'Color', 'blue', 'Width', 3, 'OverlayName', 'GNSS');
% 线性可选颜色：
% 'red', 'green', 'blue', 'white', 'cyan', 'magenta', 'yellow', 'black'

author: yandld


text:
function qr=rotqc2qr(qc)
%ROTQC2QR converts a matrix of complex quaternion row vectors into real form
%
% Inputs:
%
%     QC(2m,n)   mxn matrix of complex-valued quaternions
%
% Outputs:
%
%     QR(4m,n)   mxn matrix of real-valued quaternions
%
% The complex-valued quaternion [r+j*b  a+j*c] becomes [r a b c]

%
%      Copyright (C) Mike Brookes 2000-2006
%      Version: $Id: rotqc2qr.m,v 1.2 2007/11/18 19:38:40 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[m,n]=size(qc);
i=(1:2:2*m)-mod(0:m-1,2);
qr=zeros(2*m,n);
qr(i,:)=real(qc);
qr(i+2,:)=imag(qc);

author: decouples


text:
function circlemulti(centers,radius,nop,style)

% H=CIRCLEMULTI(CENTER,RADIUS,NOP,STYLE)
% Draws multi circles specified by centres as rows of the matrix centres.
% example: circlemulti(10*rand(10,2), 3,20,'--r')
hold on
sc=size(centers,1);
if length(radius)==1
    radius = repmat(radius,sc,1);
end
for ii=1:sc
    circle(centers(ii,:),radius(ii),nop,style);
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function val = m01Q1001(F10, F01)
%------------------------------------------------------------------------------
%
% Integrates interpolated gridfunction {F10 U F01} times y,
% interpolation is assumed piecewise constant.
% The result corresponds to a first order moment.
%
% See also: m01
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: December 12, 2000.
%  2000 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
[hx, hy] = Q1001gridfdims(F10, F01);
Qhxy = 2 * hx * hy;
%
[F10ycp, F01ycp] = Q1001ycpowp(F10, F01, 1, 0.0);
val = sum(sum( F10ycp.*F10 ));
clear F10ycp;
%
val = (val + sum(sum( F01ycp.*F01 ))) * Qhxy;
clear F01ycp;
%------------------------------------------------------------------------------


text:
function pstruct = tapas_gaussian_obs_namep(pvec)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pstruct = struct;

pstruct.ze = pvec(1);

return;

author: translationalneuromodeling

repo: tapas


text:
function list = qMRlistModel
if nargout, list = list_models; end
for im = list_models'
    header=qMRinfo(im{:}); % not compatible with octave
    if moxunit_util_platform_is_octave, return; end
    % get the first line
    ind = strfind(header,sprintf('\n'));
    header = header(1:ind-1);
    % remove space
    header = header(find(~isspace(header),1,'first'):end);
    % add tab after colons
    indcolons = strfind(header,':');
    header = ['<strong>' header(1:indcolons) '</strong>' repmat(' ',[1 max(1,30 - indcolons)]) header((indcolons+find(~isspace(header(indcolons+1:end)),1,'first')):end)];
    disp(header)
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [blockVectorX,lambda,varargout] = ...
    lobpcg(blockVectorX,operatorA,varargin)
%LOBPCG solves Hermitian partial eigenproblems using preconditioning
%
% [blockVectorX,lambda]=lobpcg(blockVectorX,operatorA)
%
% outputs the array of algebraic smallest eigenvalues lambda and
% corresponding matrix of orthonormalized eigenvectors blockVectorX of the
% Hermitian (full or sparse) operator operatorA using input matrix
% blockVectorX as an initial guess, without preconditioning, somewhat
% similar to
%
% opts.issym=1;opts.isreal=1;K=size(blockVectorX,2);
% [blockVectorX,lambda]=eigs(operatorA,K,'SR',opts);
%
% for real symmetric operator operatorA, or
%
% K=size(blockVectorX,2);[blockVectorX,lambda]=eigs(operatorA,K,'SR');
% for Hermitian operator operatorA.
%
% [blockVectorX,lambda,failureFlag]=lobpcg(blockVectorX,operatorA)
% also returns a convergence flag.
% If failureFlag is 0 then all the eigenvalues converged; otherwise not all
% converged.
%
% [blockVectorX,lambda,failureFlag,lambdaHistory,residualNormsHistory]=...
% lobpcg(blockVectorX,'operatorA','operatorB','operatorT',blockVectorY,...
% residualTolerance,maxIterations,verbosityLevel);
%
% computes smallest eigenvalues lambda and corresponding eigenvectors
% blockVectorX of the generalized eigenproblem Ax=lambda Bx, where
% Hermitian operators operatorA and operatorB are given as functions, as
% well as a preconditioner, operatorT. The operators operatorB and
% operatorT must be in addition POSITIVE DEFINITE. To compute the largest
% eigenpairs of operatorA, simply apply the code to operatorA multiplied by
% -1. The code does not involve ANY matrix factorizations of operratorA and
% operatorB, thus, e.g., it preserves the sparsity and the structure of
% operatorA and operatorB.
%
% residualTolerance and maxIterations control tolerance and max number of
% steps, and verbosityLevel = 0, 1, or 2 controls the amount of printed
% info. lambdaHistory is a matrix with all iterative lambdas, and
% residualNormsHistory are matrices of the history of 2-norms of residuals
%
% Required input:
%  * blockVectorX (class numeric) - initial approximation to eigenvectors,
%    full or sparse matrix n-by-blockSize. blockVectorX must be full rank.
%  * operatorA (class numeric, char, or function_handle) - the main operator
%    of the eigenproblem, can be a matrix, a function name, or handle
%
% Optional function input:
%   * operatorB (class numeric, char, or function_handle) - the second
%     operator, if solving a generalized eigenproblem, can be a matrix,
%      a function name, or handle; by default if empty, operatorB=I.
%   * operatorT  (class char or function_handle) - the preconditioner,
%     by default operatorT(blockVectorX)=blockVectorX.
%
% Optional constraints input:
%   blockVectorY (class numeric) - a full or sparse n-by-sizeY matrix of
%   constraints, where sizeY < n. blockVectorY must be full rank.
%   The iterations will be performed in the (operatorB-)
%   orthogonal complement of the column-space of blockVectorY.
%
% Optional scalar input parameters:
%   residualTolerance (class numeric) - tolerance, by default,
%   residualTolerance=n*sqrt(eps) maxIterations - max number of iterations,
%   by default, maxIterations = min(n,20) verbosityLevel - either 0 (no
%   info), 1, or 2 (with pictures); by default, verbosityLevel = 0.
%
% Required output: blockVectorX and lambda (both class numeric) are
% computed blockSize eigenpairs, where blockSize=size(blockVectorX,2)
% for the initial guess blockVectorX if it is full rank.
%
% Optional output: failureFlag (class integer), lambdaHistory (class numeric)
% and residualNormsHistory (class numeric) are described above.
%
% Functions operatorA(blockVectorX), operatorB(blockVectorX) and
% operatorT(blockVectorX) must support blockVectorX being a matrix, not
% just a column vector.
%
% Every iteration involves one application of operatorA and operatorB, and
% one of operatorT.
%
% Main memory requirements: 6 (9 if isempty(operatorB)=0) matrices of the
% same size as blockVectorX, 2 matrices of the same size as blockVectorY
% (if present), and two square matrices of the size 3*blockSize.
%
% In all examples below, we use the Laplacian operator in a 20x20 square
% with the mesh size 1 which can be generated in MATLAB by running
% A = delsq(numgrid('S',21)); n=size(A,1);
% or in MATLAB and Octave by
% [~,~,A] = laplacian([19,19]); n=size(A,1);
% see http://www.mathworks.com/matlabcentral/fileexchange/27279
%
% The following Example:
%
% [blockVectorX,lambda,failureFlag]=lobpcg(randn(n,8),A,1e-5,50,2);
%
% attempts to compute 8 first eigenpairs without preconditioning,
% but not all eigenpairs converge after 50 steps, so failureFlag=1.
%
% The next Example:
%
% blockVectorY=[];lambda_all=[];
% for j=1:4
%   [blockVectorX,lambda]=...
%                    lobpcg(randn(n,2),A,blockVectorY,1e-5,200,2);
%   blockVectorY=[blockVectorY,blockVectorX];
%   lambda_all=[lambda_all' lambda']'; pause;
% end
%
% attemps to compute the same 8 eigenpairs by calling the code 4 times
% with blockSize=2 using orthogonalization to the previously founded
% eigenvectors.
%
% The following Example:
%
% R=ichol(A,struct('michol','on')); precfun = @(x)R\(R'\x);
% [blockVectorX,lambda,failureFlag]=lobpcg(randn(n,8),A,[],@(x)precfun(x),1e-5,60,2);
%
% computes the same eigenpairs in less then 25 steps, so that failureFlag=0
% using the preconditioner function "precfun", defined inline. If "precfun"
% is defined as a MATLAB function in a file, the function handle
% @(x)precfun(x) can be equivalently replaced by the function name 'precfun'
% Running
%
% [blockVectorX,lambda,failureFlag]=...
%          lobpcg(randn(n,8),A,speye(n),@(x)precfun(x),1e-5,50,2);
%
% produces similar answers, but is somewhat slower and needs more memory as
% technically a generalized eigenproblem with B=I is solved here.
%
% The following Example for a mostly diagonally dominant sparse matrix A
% demonstrates different types of preconditioning, compared to the standard
% use of the main diagonal of A:
%
% clear all; close all;
% n = 1000; M = spdiags([1:n]',0,n,n); precfun=@(x)M\x;
% A=M+sprandsym(n,.1); Xini=randn(n,5); maxiter=15; tol=1e-5;
% [~,~,~,~,rnp]=lobpcg(Xini,A,tol,maxiter,1);
% [~,~,~,~,r]=lobpcg(Xini,A,[],@(x)precfun(x),tol,maxiter,1);
% subplot(2,2,1), semilogy(r'); hold on; semilogy(rnp',':>');
% title('No preconditioning (top)'); axis tight;
% M(1,2) = 2; precfun=@(x)M\x; % M is no longer symmetric
% [~,~,~,~,rns]=lobpcg(Xini,A,[],@(x)precfun(x),tol,maxiter,1);
% subplot(2,2,2),  semilogy(r'); hold on; semilogy(rns','--s');
% title('Nonsymmetric preconditioning (square)'); axis tight;
% M(1,2) = 0; precfun=@(x)M\(x+10*sin(x)); % nonlinear preconditioning
% [~,~,~,~,rnl]=lobpcg(Xini,A,[],@(x)precfun(x),tol,maxiter,1);
% subplot(2,2,3),  semilogy(r'); hold on; semilogy(rnl','-.*');
% title('Nonlinear preconditioning (star)'); axis tight;
% M=abs(M-3.5*speye(n,n)); precfun=@(x)M\x;
% [~,~,~,~,rs]=lobpcg(Xini,A,[],@(x)precfun(x),tol,maxiter,1);
% subplot(2,2,4),  semilogy(r'); hold on; semilogy(rs','-d');
% title('Selective preconditioning (diamond)'); axis tight;
%
% Revision 4.16 adds support for distributed or codistributed arrays
% available in MATLAB BigData toolbox, e.g.,
%
% A = codistributed(diag(1:100)); B = codistributed(diag(101:200));
% [blockVectorX,lambda]=lobpcg(randn(100,2),A,1e-5,5,2)
%
% Revision 4.17 adds support for single precision, e.g.,
% A = diag(1:100); B = single(diag(101:200));
% [blockVectorX,lambda]=lobpcg(randn(100,2),A,1e-5,15,2);
% A = diag(1:100); B = diag(101:200);
% [blockVectorX,lambda]=lobpcg(randn(100,2,'single'),A,1e-5,15,2);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This main function LOBPCG is a version of
% the preconditioned conjugate gradient method (Algorithm 5.1) described in
% A. V. Knyazev, Toward the Optimal Preconditioned Eigensolver:
% Locally Optimal Block Preconditioned Conjugate Gradient Method,
% SIAM Journal on Scientific Computing 23 (2001), no. 2, pp. 517-541.
% http://dx.doi.org/10.1137/S1064827500366124
%
% Known bugs/features:
%
% - an excessively small requested tolerance may result in often restarts
% and instability. The code is not written to produce an eps-level
% accuracy! Use common sense.
%
% - the code may be very sensitive to the number of eigenpairs computed,
% if there is a cluster of eigenvalues not completely included, cf.
%
% operatorA=diag([1 1.99 2:99]);
% [blockVectorX,lambda]=lobpcg(randn(100,1),operatorA,1e-10,80,2);
% [blockVectorX,lambda]=lobpcg(randn(100,2),operatorA,1e-10,80,2);
% [blockVectorX,lambda]=lobpcg(randn(100,3),operatorA,1e-10,80,2);
%
% - using a nonsymmetric preconditioner is possible, but may be unstable
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The main distribution site:
% https://github.com/lobpcg/blopex
%
% A C-version of this code is a part of the
% https://github.com/lobpcg/blopex
% package and is directly available, e.g., in SLEPc and HYPRE.
%
% A python version of this code is in
% https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.lobpcg.html
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   License:  MIT / Apache-2.0
%   Copyright (c) 2000-2019 A.V. Knyazev, Andrew.Knyazev@ucdenver.edu
%   $Revision: 1.2 $  $Date: 13-June-2019
%   This revision is tested in 9.6.0.1114505 (R2019a) Update 2, but is
%   expected to work on any >R2007b MATLAB.
%   Revision 4.13 tested in MATLAB 6.5-7.13.
%   Revision 4.13 tested and available in Octave 3.2.3-3.4.2, see
%   https://octave.sourceforge.io/linear-algebra/function/lobpcg.html
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Begin
% Function gather defined to be identity if nonexistent, before 2016a
if exist("gather", "file") == 2
    mygather=@(x)gather(x);
else
    mygather=@(x)x;
end
% constants
CONVENTIONAL_CONSTRAINTS = 1;
SYMMETRIC_CONSTRAINTS = 2;
%Initial settings
failureFlag = 1;
if nargin < 2
    error('BLOPEX:lobpcg:NotEnoughInputs',...
        strcat('There must be at least 2 input agruments: ',...
        'blockVectorX and operatorA'));
end
if nargin > 8
    warning('BLOPEX:lobpcg:TooManyInputs',...
        strcat('There must be at most 8 input agruments ',...
        'unless arguments are passed to a function'));
end
if ~isnumeric(blockVectorX)
    error('BLOPEX:lobpcg:FirstInputNotNumeric',...
        'The first input argument blockVectorX must be numeric');
end
[n,blockSize]=size(blockVectorX);
if blockSize > n
    error('BLOPEX:lobpcg:FirstInputFat',...
        'The first input argument blockVectorX must be tall, not fat');
end
if n < 6
    error('BLOPEX:lobpcg:MatrixTooSmall',...
        'The code does not work for matrices of small sizes');
end
if isa(operatorA,'numeric')
    nA = size(operatorA,1);
    if any(size(operatorA) ~= nA)
        error('BLOPEX:lobpcg:MatrixNotSquare',...
            'operatorA must be a square matrix or a string');
    end
    if size(operatorA) ~= n
        error('BLOPEX:lobpcg:MatrixWrongSize',...
            ['The size ' int2str(size(operatorA))...
            ' of operatorA is not the same as ' int2str(n)...
            ' - the number of rows of blockVectorX']);
    end
end
count_string = 0;
operatorT = [];
operatorB = [];
residualTolerance = [];
maxIterations = [];
verbosityLevel = [];
blockVectorY = []; sizeY = 0;
for j = 1:nargin-2
    if isequal(size(varargin{j}),[n,n])
        if isempty(operatorB)
            operatorB = varargin{j};
        else
            error('BLOPEX:lobpcg:TooManyMatrixInputs',...
                strcat('Too many matrix input arguments. ',...
                'Preconditioner operatorT must be an M-function'));
        end
    elseif isequal(size(varargin{j},1),n) && size(varargin{j},2) < n
        if isempty(blockVectorY)
            blockVectorY = varargin{j};
            sizeY=size(blockVectorY,2);
        else
            error('BLOPEX:lobpcg:WrongConstraintsFormat',...
                'Something wrong with blockVectorY input argument');
        end
    elseif ischar(varargin{j}) || isa(varargin{j},'function_handle')
        if count_string == 0
            if isempty(operatorB)
                operatorB = varargin{j};
                count_string = count_string + 1;
            else
                operatorT = varargin{j};
            end
        elseif count_string == 1
            operatorT = varargin{j};
        else
            warning('BLOPEX:lobpcg:TooManyStringFunctionHandleInputs',...
                'Too many string or FunctionHandle input arguments');
        end
    elseif isequal(size(varargin{j}),[n,n])
        error('BLOPEX:lobpcg:WrongPreconditionerFormat',...
            'Preconditioner operatorT must be an M-function');
    elseif max(size(varargin{j})) == 1
        if isempty(residualTolerance)
            residualTolerance = varargin{j};
        elseif isempty(maxIterations)
            maxIterations = varargin{j};
        elseif isempty(verbosityLevel)
            verbosityLevel = varargin{j};
        else
            warning('BLOPEX:lobpcg:TooManyScalarInputs',...
                'Too many scalar parameters, need only three');
        end
    elseif isempty(varargin{j})
        if isempty(operatorB)
            count_string = count_string + 1;
        elseif ~isempty(operatorT)
            count_string = count_string + 1;
        elseif ~isempty(blockVectorY)
            error('BLOPEX:lobpcg:UnrecognizedEmptyInput',...
                ['Unrecognized empty input argument number ' int2str(j+2)]);
        end
    else
        error('BLOPEX:lobpcg:UnrecognizedInput',...
            ['Input argument number ' int2str(j+2) ' not recognized.']);
    end
end
if verbosityLevel
    if issparse(blockVectorX)
        fprintf(['The sparse initial guess with %i colunms '...
            'and %i raws is detected  \n'],n,blockSize);
    else
        fprintf(['The full initial guess with %i colunms '...
            'and %i raws is detected  \n'],n,blockSize);
    end
    if ischar(operatorA)
        fprintf('The main operator is detected as an M-function %s \n',...
            operatorA);
    elseif isa(operatorA,'function_handle')
        fprintf('The main operator is detected as an M-function %s \n',...
            func2str(operatorA));
    elseif issparse(operatorA)
        fprintf('The main operator is detected as a sparse matrix \n');
    else
        fprintf('The main operator is detected as a full matrix \n');
    end
    if isempty(operatorB)
        fprintf('Solving standard eigenvalue problem, not generalized \n');
    elseif ischar(operatorB)
        fprintf(['The second operator of the generalized eigenproblem \n'...
            'is detected as an M-function %s \n'],operatorB);
    elseif isa(operatorB,'function_handle')
        fprintf(['The second operator of the generalized eigenproblem \n'...
            'is detected as an M-function %s \n'],func2str(operatorB));
    elseif issparse(operatorB)
        fprintf(strcat('The second operator of the generalized',...
            'eigenproblem \n is detected as a sparse matrix \n'));
    else
        fprintf(strcat('The second operator of the generalized',...
            'eigenproblem \n is detected as a full matrix \n'));
    end
    if isempty(operatorT)
        fprintf('No preconditioner is detected \n');
    elseif ischar(operatorT)
        fprintf('The preconditioner is detected as an M-function %s \n',...
            operatorT);
    elseif isa(operatorT,'function_handle')
        fprintf('The preconditioner is detected as an M-function %s \n',...
            func2str(operatorT));
    end
    if isempty(blockVectorY)
        fprintf('No matrix of constraints is detected \n')
    elseif issparse(blockVectorY)
        fprintf('The sparse matrix of %i constraints is detected \n',sizeY);
    else
        fprintf('The full matrix of %i constraints is detected \n',sizeY);
    end
    if issparse(blockVectorY) ~= issparse(blockVectorX)
        warning('BLOPEX:lobpcg:SparsityInconsistent',...
            strcat('The sparsity formats of the initial guess and ',...
            'the constraints are inconsistent'));
    end
end
% Set defaults
if isempty(residualTolerance)
    residualTolerance = sqrt(eps)*n;
end
if isempty(maxIterations)
    maxIterations = min(n,20);
end
if isempty(verbosityLevel)
    verbosityLevel = 0;
end
if verbosityLevel
    fprintf('Tolerance %e and maximum number of iterations %i \n',...
        residualTolerance,maxIterations)
end
%constraints preprocessing
if isempty(blockVectorY)
    constraintStyle = 0;
else
    %    constraintStyle = SYMMETRIC_CONSTRAINTS; % more accurate?
    constraintStyle = CONVENTIONAL_CONSTRAINTS;
end
if constraintStyle == CONVENTIONAL_CONSTRAINTS

    if isempty(operatorB)
        gramY = blockVectorY'*blockVectorY;
    else
        if isnumeric(operatorB)
            blockVectorBY = operatorB*blockVectorY;
        else
            blockVectorBY = feval(operatorB,blockVectorY);
        end
        gramY=blockVectorY'*blockVectorBY;
    end
    gramY=(gramY'+gramY)*0.5;
    if isempty(operatorB)
        blockVectorX = blockVectorX - ...
            blockVectorY*(gramY\(blockVectorY'*blockVectorX));
    else
        blockVectorX =blockVectorX - ...
            blockVectorY*(gramY\(blockVectorBY'*blockVectorX));
    end

elseif constraintStyle == SYMMETRIC_CONSTRAINTS

    if ~isempty(operatorB)
        if isnumeric(operatorB)
            blockVectorY = operatorB*blockVectorY;
        else
            blockVectorY = feval(operatorB,blockVectorY);
        end
    end
    if isempty(operatorT)
        gramY = blockVectorY'*blockVectorY;
    else
        blockVectorTY = feval(operatorT,blockVectorY);
        gramY = blockVectorY'*blockVectorTY;
    end
    gramY=(gramY'+gramY)*0.5;
    if isempty(operatorT)
        blockVectorX = blockVectorX - ...
            blockVectorY*(gramY\(blockVectorY'*blockVectorX));
    else
        blockVectorX = blockVectorX - ...
            blockVectorTY*(gramY\(blockVectorY'*blockVectorX));
    end

end
%Making the initial vectors (operatorB-) orthonormal
if isempty(operatorB)
    %[blockVectorX,gramXBX] = qr(blockVectorX,0);
    gramXBX=mygather(blockVectorX'*blockVectorX);
    if ~isreal(gramXBX)
        gramXBX=(gramXBX+gramXBX')*0.5;
    end
    [gramXBX,cholFlag]=chol(gramXBX);
    if  cholFlag ~= 0
        error('BLOPEX:lobpcg:ConstraintsTooTight',...
            'The initial approximation after constraints is not full rank');
    end
    blockVectorX = blockVectorX/gramXBX;
else
    %[blockVectorX,blockVectorBX] = orth(operatorB,blockVectorX);
    if isnumeric(operatorB)
        blockVectorBX = operatorB*blockVectorX;
    else
        blockVectorBX = feval(operatorB,blockVectorX);
    end
    gramXBX=blockVectorX'*blockVectorBX;
    if ~isreal(gramXBX)
        gramXBX=(gramXBX+gramXBX')*0.5;
    end
    [gramXBX,cholFlag]=chol(gramXBX);
    if  cholFlag ~= 0
        error('BLOPEX:lobpcg:InitialNotFullRank',...
            '%s\n%s', ...
            'The initial approximation after constraints is not ',...
            'full rank or/and operatorB is not positive definite');
    end
    blockVectorX = blockVectorX/gramXBX;
    blockVectorBX = blockVectorBX/gramXBX;
end
% Checking if the problem is big enough for the algorithm,
% i.e. n-sizeY > 5*blockSize
% Theoretically, the algorithm should be able to run if
% n-sizeY > 3*blockSize,
% but the extreme cases might be unstable, so we use 5 instead of 3 here.
if n-sizeY < 5*blockSize
    error('BLOPEX:lobpcg:MatrixTooSmall','%s\n%s', ...
        'The problem size is too small, relative to the block size.',...
        'Try using eig() or eigs() instead.');
end
% Preallocation
residualNormsHistory=zeros(blockSize,maxIterations);
lambdaHistory=zeros(blockSize,maxIterations+1);
condestGhistory=zeros(1,maxIterations+1);
blockVectorAR=zeros(n,blockSize, 'like', blockVectorX);
blockVectorP=zeros(n,blockSize, 'like', blockVectorX);
blockVectorAP=zeros(n,blockSize, 'like', blockVectorX);
if ~isempty(operatorB)
    blockVectorBR=zeros(n,blockSize, 'like', blockVectorX);
    blockVectorBP=zeros(n,blockSize, 'like', blockVectorX);
end
%Initial settings for the loop
if isnumeric(operatorA)
    blockVectorAX = operatorA*blockVectorX;
else
    blockVectorAX = feval(operatorA,blockVectorX);
end
gramXAX = full(blockVectorX'*blockVectorAX);
gramXAX = (gramXAX + gramXAX')*0.5;
% eig(...,'chol') uses only the diagonal and upper triangle -
% not true in MATLAB
% Octave v3.2.3-4, eig() does not support inputting 'chol'
[coordX,gramXAX]=eig(gramXAX,eye(blockSize));
lambda=diag(gramXAX); %eig returns non-ordered eigenvalues on the diagonal
if issparse(blockVectorX)
    coordX=sparse(coordX);
end
blockVectorX  =  blockVectorX*coordX;
blockVectorAX = blockVectorAX*coordX;
if ~isempty(operatorB)
    blockVectorBX = blockVectorBX*coordX;
end
clear coordX
condestGhistory(1)=-log10(eps)/2;  %if too small cause unnecessary restarts
lambdaHistory(1:blockSize,1) = mygather(lambda);
activeMask = true(blockSize,1);
% currentBlockSize = blockSize; %iterate all
%
% restart=1;%steepest descent
%The main part of the method is the loop of the CG method: begin
for iterationNumber=1:maxIterations

    %     %Computing the active residuals
    %     if isempty(operatorB)
    %%         if currentBlockSize > 1
    %%             blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ...
    %%                 blockVectorX(:,activeMask)*spdiags(lambda(activeMask),0,currentBlockSize,currentBlockSize);
    %%         else
    %%             blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ...
    %%                 blockVectorX(:,activeMask)*lambda(activeMask);
    %%         end
    %         blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ...
    %         bsxfun(@times,blockVectorX(:,activeMask),lambda(activeMask)');
    %     else
    %%         if currentBlockSize > 1
    %%             blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ...
    %%                 blockVectorBX(:,activeMask)*spdiags(lambda(activeMask),0,currentBlockSize,currentBlockSize);
    %%         else
    %%             blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ...
    %%                 blockVectorBX(:,activeMask)*lambda(activeMask);
    %%         end
    %         blockVectorR(:,activeMask)=blockVectorAX(:,activeMask) - ...
    %         bsxfun(@times,blockVectorBX(:,activeMask),lambda(activeMask)');
    %     end

    %Computing all residuals
    if isempty(operatorB)
        %         if blockSize > 1
        %             blockVectorR = blockVectorAX - ...
        %                 blockVectorX*spdiags(lambda,0,blockSize,blockSize);
        %         else
        %             blockVectorR = blockVectorAX - blockVectorX*lambda;
        %             %to make blockVectorR full when lambda is just a scalar
        %         end
        blockVectorR = blockVectorAX - ...
            bsxfun(@times,blockVectorX,lambda');
    else
        %        if blockSize > 1
        %             blockVectorR = blockVectorAX - ...
        %                 blockVectorBX*spdiags(lambda,0,blockSize,blockSize);
        %        else
        %             blockVectorR = blockVectorAX - blockVectorBX*lambda;
        %        end
        blockVectorR = blockVectorAX - ...
            bsxfun(@times,blockVectorBX,lambda');
    end

    %Satisfying the constraints for the active residulas
    if constraintStyle == SYMMETRIC_CONSTRAINTS
        if isempty(operatorT)
            blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
                blockVectorY*(gramY\(blockVectorY'*...
                blockVectorR(:,activeMask)));
        else
            blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
                blockVectorY*(gramY\(blockVectorTY'*...
                blockVectorR(:,activeMask)));
        end
    end

    residualNorms = full(sqrt(sum(conj(blockVectorR).*blockVectorR)'));
    residualNormsHistory(1:blockSize,iterationNumber) = ...
        mygather(residualNorms);

    %index antifreeze
    activeMask = full(residualNorms > residualTolerance) & activeMask;
    %activeMask = full(residualNorms > residualTolerance);
    %above allows vectors back into active, which causes problems with frosen Ps
    %activeMask = full(residualNorms > 0);      %iterate all, ignore freeze

    currentBlockSize = mygather(sum(activeMask));
    if  currentBlockSize == 0
        failureFlag=0; %all eigenpairs converged
        break
    end

    %Applying the preconditioner operatorT to the active residulas
    if ~isempty(operatorT)
        blockVectorR(:,activeMask) = ...
            feval(operatorT,blockVectorR(:,activeMask));
    end

    if constraintStyle == CONVENTIONAL_CONSTRAINTS
        if isempty(operatorB)
            blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
                blockVectorY*(gramY\(blockVectorY'*...
                blockVectorR(:,activeMask)));
        else
            blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
                blockVectorY*(gramY\(blockVectorBY'*...
                blockVectorR(:,activeMask)));
        end
    end

    %Making active (preconditioned) residuals orthogonal to blockVectorX
    if isempty(operatorB)
        blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
            blockVectorX*(blockVectorX'*blockVectorR(:,activeMask));
    else
        blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
            blockVectorX*(blockVectorBX'*blockVectorR(:,activeMask));
    end

    %Making active residuals orthonormal
    if isempty(operatorB)
        %[blockVectorR(:,activeMask),gramRBR]=...
        %qr(blockVectorR(:,activeMask),0); %to increase stability
        gramRBR=blockVectorR(:,activeMask)'*blockVectorR(:,activeMask);
        if ~isreal(gramRBR)
            gramRBR=(gramRBR+gramRBR')*0.5;
        end
        [gramRBR,cholFlag]=chol(gramRBR);
        if  cholFlag == 0
            blockVectorR(:,activeMask) = blockVectorR(:,activeMask)/gramRBR;
        else
            warning('BLOPEX:lobpcg:ResidualNotFullRank',...
                'The residual is not full rank.');
            break
        end
    else
        if isnumeric(operatorB)
            blockVectorBR(:,activeMask) = ...
                operatorB*blockVectorR(:,activeMask);
        else
            blockVectorBR(:,activeMask) = ...
                feval(operatorB,blockVectorR(:,activeMask));
        end
        gramRBR=blockVectorR(:,activeMask)'*blockVectorBR(:,activeMask);
        if ~isreal(gramRBR)
            gramRBR=(gramRBR+gramRBR')*0.5;
        end
        [gramRBR,cholFlag]=chol(gramRBR);
        if  cholFlag == 0
            blockVectorR(:,activeMask) = ...
                blockVectorR(:,activeMask)/gramRBR;
            blockVectorBR(:,activeMask) = ...
                blockVectorBR(:,activeMask)/gramRBR;
        else
            warning('BLOPEX:lobpcg:ResidualNotFullRankOrElse',...
                strcat('The residual is not full rank or/and operatorB ',...
                'is not positive definite.'));
            break
        end
    end
    clear gramRBR;

    if isnumeric(operatorA)
        blockVectorAR(:,activeMask) = ...
            mygather(operatorA*blockVectorR(:,activeMask));
    else
        blockVectorAR(:,activeMask) = ...
            feval(operatorA,blockVectorR(:,activeMask));
    end

    condestGmean = mean(condestGhistory(max(1,iterationNumber-10-...
        round(log(currentBlockSize))):iterationNumber));

    %  restart=1;

    % The Raileight-Ritz method for [blockVectorX blockVectorR blockVectorP]
    if isa(blockVectorAR, 'single')    % single initial
        myeps = 1; % play safe
    elseif isa(blockVectorR, 'single') %single somethings else
        myeps = eps(single(1));
    else                               % double everything
        myeps = eps;
    end
    if  mygather(residualNorms) > myeps^0.6
        explicitGramFlag = 0;
    else
        explicitGramFlag = 1;  %suggested by Garrett Moran, private
    end

    activeRSize=size(blockVectorR(:,activeMask),2);
    if iterationNumber == 1
        activePSize=0;
        restart=1;
    else
        activePSize=size(blockVectorP(:,activeMask),2);
        restart=0;
    end

    gramXAR=full(blockVectorAX'*blockVectorR(:,activeMask));
    gramRAR=full(blockVectorAR(:,activeMask)'*blockVectorR(:,activeMask));
    gramRAR=(gramRAR'+gramRAR)*0.5;

    if explicitGramFlag
        gramXAX=full(blockVectorAX'*blockVectorX);
        gramXAX=(gramXAX'+gramXAX)*0.5;
        if isempty(operatorB)
            gramXBX=full(blockVectorX'*blockVectorX);
            gramRBR=full(blockVectorR(:,activeMask)'*...
                blockVectorR(:,activeMask));
            gramXBR=full(blockVectorX'*blockVectorR(:,activeMask));
        else
            gramXBX=full(blockVectorBX'*blockVectorX);
            gramRBR=full(blockVectorBR(:,activeMask)'*...
                blockVectorR(:,activeMask));
            gramXBR=full(blockVectorBX'*blockVectorR(:,activeMask));
        end
        gramXBX=(gramXBX'+gramXBX)*0.5;
        gramRBR=(gramRBR'+gramRBR)*0.5;

    end

    if iterationNumber > 1
        %Making active conjugate directions orthonormal
        if isempty(operatorB)
            %[blockVectorP(:,activeMask),gramPBP] = qr(blockVectorP(:,activeMask),0);
            gramPBP=blockVectorP(:,activeMask)'*blockVectorP(:,activeMask);
            if ~isreal(gramPBP)
                gramPBP=(gramPBP+gramPBP')*0.5;
            end
            [gramPBP,cholFlag]=chol(gramPBP);
            if  cholFlag == 0
                blockVectorP(:,activeMask) = ...
                    blockVectorP(:,activeMask)/gramPBP;
                blockVectorAP(:,activeMask) = ...
                    blockVectorAP(:,activeMask)/gramPBP;
                restart = 0;
            else
                warning('BLOPEX:lobpcg:DirectionNotFullRank',...
                    'The direction matrix is not full rank.');
                restart = 1;
            end
        else
            gramPBP=blockVectorP(:,activeMask)'*blockVectorBP(:,activeMask);
            if ~isreal(gramPBP)
                gramPBP=(gramPBP+gramPBP')*0.5;
            end
            [gramPBP,cholFlag]=chol(gramPBP);
            if  cholFlag == 0
                blockVectorP(:,activeMask) = ...
                    blockVectorP(:,activeMask)/gramPBP;
                blockVectorAP(:,activeMask) = ...
                    blockVectorAP(:,activeMask)/gramPBP;
                blockVectorBP(:,activeMask) = ...
                    blockVectorBP(:,activeMask)/gramPBP;
                restart = 0;
            else
                warning('BLOPEX:lobpcg:DirectionNotFullRank',...
                    strcat('The direction matrix is not full rank ',...
                    'or/and operatorB is not positive definite.'));
                restart = 1;
            end
        end
        clear gramPBP
    end

    for cond_try=1:2           %cond_try == 2 when restart

        if ~restart
            gramXAP=full(blockVectorAX'*blockVectorP(:,activeMask));
            gramRAP=full(blockVectorAR(:,activeMask)'*...
                blockVectorP(:,activeMask));
            gramPAP=full(blockVectorAP(:,activeMask)'*...
                blockVectorP(:,activeMask));
            gramPAP=(gramPAP'+gramPAP)*0.5;

            if explicitGramFlag
                gramA = [ gramXAX     gramXAR     gramXAP
                    gramXAR'    gramRAR     gramRAP
                    gramXAP'     gramRAP'    gramPAP ];
            else
                gramA = [ diag(lambda)  gramXAR  gramXAP
                    gramXAR'      gramRAR  gramRAP
                    gramXAP'      gramRAP'  gramPAP ];
            end

            clear gramXAP  gramRAP gramPAP

            if isempty(operatorB)
                gramXBP=full(blockVectorX'*blockVectorP(:,activeMask));
                gramRBP=full(blockVectorR(:,activeMask)'*...
                    blockVectorP(:,activeMask));
            else
                gramXBP=full(blockVectorBX'*blockVectorP(:,activeMask));
                gramRBP=full(blockVectorBR(:,activeMask)'*...
                    blockVectorP(:,activeMask));
                %or blockVectorR(:,activeMask)'*blockVectorBP(:,activeMask);
            end

            if explicitGramFlag
                if isempty(operatorB)
                    gramPBP=full(blockVectorP(:,activeMask)'*...
                        blockVectorP(:,activeMask));
                else
                    gramPBP=full(blockVectorBP(:,activeMask)'*...
                        blockVectorP(:,activeMask));
                end
                gramPBP=(gramPBP'+gramPBP)*0.5;
                gramB = [ gramXBX  gramXBR  gramXBP
                    gramXBR' gramRBR  gramRBP
                    gramXBP' gramRBP' gramPBP ];
                clear   gramPBP
            else
                gramB=[eye(blockSize) zeros(blockSize,activeRSize) gramXBP
                    zeros(blockSize,activeRSize)' eye(activeRSize) gramRBP
                    gramXBP' gramRBP' eye(activePSize) ];
            end

            clear gramXBP  gramRBP;

        else

            if explicitGramFlag
                gramA = [ gramXAX   gramXAR
                    gramXAR'    gramRAR  ];
                gramB = [ gramXBX  gramXBR
                    gramXBR' eye(activeRSize)  ];
                clear gramXAX gramXBX gramXBR
            else
                gramA = [ diag(lambda)  gramXAR
                    gramXAR'        gramRAR  ];
                gramB = eye(blockSize+activeRSize);
            end

            clear gramXAR gramRAR;

        end

        condestG = log10(cond(gramB))+1;
        if (condestG/condestGmean > 2 && condestG > 2 )|| condestG > 8
            %black magic - need to guess the restart
            if verbosityLevel
                fprintf('Restart on step %i as condestG %5.4e \n',...
                    iterationNumber,condestG);
            end
            if cond_try == 1 && ~restart
                restart=1; %steepest descent restart for stability
            else
                warning('BLOPEX:lobpcg:IllConditioning',...
                    'Gramm matrix ill-conditioned: results unpredictable');
            end
        else
            break
        end

    end

    [gramA,gramB]=eig(gramA,gramB);
    lambda=diag(gramB(1:blockSize,1:blockSize));
    coordX=gramA(:,1:blockSize);

    clear gramA gramB

    if issparse(blockVectorX)
        coordX=sparse(coordX);
    end

    if ~restart
        blockVectorP =  blockVectorR(:,activeMask)*...
            coordX(blockSize+1:blockSize+activeRSize,:) + ...
            blockVectorP(:,activeMask)*...
            coordX(blockSize+activeRSize+1:blockSize + ...
            activeRSize+activePSize,:);
        blockVectorAP = blockVectorAR(:,activeMask)*...
            coordX(blockSize+1:blockSize+activeRSize,:) + ...
            blockVectorAP(:,activeMask)*...
            coordX(blockSize+activeRSize+1:blockSize + ...
            activeRSize+activePSize,:);
        if ~isempty(operatorB)
            blockVectorBP = blockVectorBR(:,activeMask)*...
                coordX(blockSize+1:blockSize+activeRSize,:) + ...
                blockVectorBP(:,activeMask)*...
                coordX(blockSize+activeRSize+1:blockSize+activeRSize+activePSize,:);
        end
    else %use block steepest descent
        blockVectorP =   blockVectorR(:,activeMask)*...
            coordX(blockSize+1:blockSize+activeRSize,:);
        blockVectorAP = blockVectorAR(:,activeMask)*...
            coordX(blockSize+1:blockSize+activeRSize,:);
        if ~isempty(operatorB)
            blockVectorBP = blockVectorBR(:,activeMask)*...
                coordX(blockSize+1:blockSize+activeRSize,:);
        end
    end

    blockVectorX = blockVectorX*coordX(1:blockSize,:) + blockVectorP;
    blockVectorAX = blockVectorAX*coordX(1:blockSize,:) + blockVectorAP;
    if ~isempty(operatorB)
        blockVectorBX = blockVectorBX*coordX(1:blockSize,:) + blockVectorBP;
    end
    clear coordX
    %%end RR

    lambdaHistory(1:blockSize,iterationNumber+1) = mygather(lambda);
    condestGhistory(iterationNumber+1) = mygather(condestG);

    if verbosityLevel
        fprintf('Iteration %i current block size %i \n',...
            iterationNumber,currentBlockSize);
        fprintf('Eigenvalues lambda %17.16e \n',mygather(lambda));
        fprintf('Residual Norms %e \n',mygather(residualNorms'));
    end
end
%The main step of the method was the CG cycle: end
%Postprocessing
%Making sure blockVectorX's "exactly" satisfy the blockVectorY constrains??
%Making sure blockVectorX's are "exactly" othonormalized by final "exact" RR
if isempty(operatorB)
    gramXBX=full(blockVectorX'*blockVectorX);
else
    if isnumeric(operatorB)
        blockVectorBX = operatorB*blockVectorX;
    else
        blockVectorBX = feval(operatorB,blockVectorX);
    end
    gramXBX = full(blockVectorX'*blockVectorBX);
end
gramXBX=(gramXBX'+gramXBX)*0.5;
if isnumeric(operatorA)
    blockVectorAX = operatorA*blockVectorX;
else
    blockVectorAX = feval(operatorA,blockVectorX);
end
gramXAX = full(blockVectorX'*blockVectorAX);
gramXAX = (gramXAX + gramXAX')*0.5;
%Raileigh-Ritz for blockVectorX, which is already operatorB-orthonormal
[coordX,gramXBX] = eig(gramXAX,gramXBX);
lambda=diag(gramXBX);
if issparse(blockVectorX)
    coordX=sparse(coordX);
end
blockVectorX  =   blockVectorX*coordX;
blockVectorAX  =  blockVectorAX*coordX;
if ~isempty(operatorB)
    blockVectorBX  =  blockVectorBX*coordX;
end
%Computing all residuals
if isempty(operatorB)
    %     if blockSize > 1
    %         blockVectorR = blockVectorAX - ...
    %             blockVectorX*spdiags(lambda,0,blockSize,blockSize);
    %     else
    %         blockVectorR = blockVectorAX - blockVectorX*lambda;
    %     end
    blockVectorR = blockVectorAX - ...
        bsxfun(@times,blockVectorX,lambda');
else
    %     if blockSize > 1
    %         blockVectorR=blockVectorAX - ...
    %             blockVectorBX*spdiags(lambda,0,blockSize,blockSize);
    %     else
    %         blockVectorR = blockVectorAX - blockVectorBX*lambda;
    %     end
    blockVectorR = blockVectorAX - ...
        bsxfun(@times,blockVectorBX,lambda');
end
residualNorms=full(sqrt(sum(conj(blockVectorR).*blockVectorR)'));
residualNormsHistory(1:blockSize,iterationNumber) = ...
    mygather(residualNorms);
if verbosityLevel
    fprintf('Final Eigenvalues lambda %17.16e \n',mygather(lambda));
    fprintf('Final Residual Norms %e \n',mygather(residualNorms'));
end
if verbosityLevel == 2
    whos
    figure(491)
    semilogy((abs(residualNormsHistory(1:blockSize,1:iterationNumber-1)))');
    title('Residuals for Different Eigenpairs','fontsize',16);
    ylabel('Eucledian norm of residuals','fontsize',16);
    xlabel('Iteration number','fontsize',16);
    %axis tight;
    %axis([0 maxIterations+1 1e-15 1e3])
    set(gca,'FontSize',14);
    figure(492);
    semilogy(abs((lambdaHistory(1:blockSize,1:iterationNumber)-...
        repmat(mygather(lambda),1,iterationNumber)))');
    title('Eigenvalue errors for Different Eigenpairs','fontsize',16);
    ylabel('Estimated eigenvalue errors','fontsize',16);
    xlabel('Iteration number','fontsize',16);
    %axis tight;
    %axis([0 maxIterations+1 1e-15 1e3])
    set(gca,'FontSize',14);
    drawnow;
end
varargout(1)={failureFlag};
varargout(2)={lambdaHistory(1:blockSize,1:iterationNumber)};
varargout(3)={residualNormsHistory(1:blockSize,1:iterationNumber-1)};
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function value = r4_mach ( i )

%*****************************************************************************80
%
%% R4_MACH returns single precision real machine constants.
%
%  Discussion:
%
%    Assume that single precision real numbers are stored with a mantissa
%    of T digits in base B, with an exponent whose value must lie
%    between EMIN and EMAX.  Then for values of I between 1 and 5,
%    R1MACH will return the following values:
%
%      R1MACH(1) = B^(EMIN-1), the smallest positive magnitude.
%      R1MACH(2) = B^EMAX*(1-B^(-T)), the largest magnitude.
%      R1MACH(3) = B^(-T), the smallest relative spacing.
%      R1MACH(4) = B^(1-T), the largest relative spacing.
%      R1MACH(5) = log10(B)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 April 2007
%
%  Author:
%
%    Original FORTRAN77 version by Phyllis Fox, Andrew Hall, Norman Schryer
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Phyllis Fox, Andrew Hall, Norman Schryer,
%    Algorithm 528,
%    Framework for a Portable Library,
%    ACM Transactions on Mathematical Software,
%    Volume 4, Number 2, June 1978, page 176-188.
%
%  Parameters:
%
%    Input, integer I, chooses the parameter to be returned.
%    1 <= I <= 5.
%
%    Output, real VALUE, the value of the chosen parameter.
%
  if ( i < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_MACH - Fatal error!\n' );
    fprintf ( 1, '  The input argument I is out of bounds.\n' );
    fprintf ( 1, '  Legal values satisfy 1 <= I <= 5.\n' );
    fprintf ( 1, '  I = %d\n', i );
    value = 0.0;
    error ( 'R4_MACH - Fatal error!' );
  elseif ( i == 1 )
    value = 1.1754944E-38;
  elseif ( i == 2 )
    value = 3.4028235E+38;
  elseif ( i == 3 )
    value = 5.9604645E-08;
  elseif ( i == 4 )
    value = 1.1920929E-07;
  elseif ( i == 5 )
    value = 0.3010300;
  elseif ( 5 < i )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_MACH - Fatal error!\n' );
    fprintf ( 1, '  The input argument I is out of bounds.\n' );
    fprintf ( 1, '  Legal values satisfy 1 <= I <= 5.\n' );
    fprintf ( 1, '  I = %d\n', i );
    value = 0.0;
    error ( 'R4_MACH - Fatal error!' );
  end

  return
end


text:
function prev_rng = seed_rand(seed)
% seed_rand - Set random number generator to a fixed seed.
%   prev_rng = seed_rand(seed)
%
%   Strategic use ensures that results are reproducible.
%
%   To restore the previous rng after calling this do:
%   rng(prev_rng);

% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

if nargin < 1
    % This value works best for me.
    seed = 3;
    % Just kidding, of course ;-).
end

prev_rng = rng;
rng(seed, 'twister')

author: ShaoqingRen


text:
function domi_factor = dominant_slope(data,fs)
% decide whch side has the dominant slope (up or down)
%
% return domi_factor = 1, up slope is dominant
%                     -1, down slope is dominant
% so the data may be multiplied by domi_factor

domi_factor=1;
up_slope=0;
down_slope=0;
% split data to every 2s and find the median amplitude of each segment
d=reshape(data(1:floor(length(data)/(fs*2))*fs*2),fs*2,[]);
amplitude=median(max(d)-min(d));
% use amplitude/3 as MinPeakProminence to find peaks
[pks locs]=findpeaks(data,'MinPeakProminence',amplitude/3);
% calculate x-coordinate of each slope
for i=1:length(pks)-1
    [mn mni]=min(data(locs(i):locs(i+1)));
    l=length(data(locs(i):locs(i+1)));
    up_slope(i)=l-mni;
    down_slope(i)=mni-1;
end
% sum(up_slope)
% sum(down_slope)
% median(up_slope)
% median(down_slope)
if median(up_slope)>median(down_slope)
    domi_factor=-1;
end
if median(up_slope)==median(down_slope)
    if sum(up_slope)>sum(down_slope)
        domi_factor=-1;
    end
end
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox