text:
function popt = partition_options(pdefault, varargin)
%PARTITION_OPTIONS Options for EQ partition
%
%Syntax
% popt = partition_options(pdefault,options);
%
%Description
% POPT = PARTITION_OPTIONS(PDEFAULT,options) collects partition options,
% specified as name, value pairs, and places these into the structure POPT.
% The structure PDEFAULT is used to define default option values.
%
% The structures pdefault and popt may contain the following fields:
% extra_offset:  boolean
%
% The following partition options are available.
%
% 'offset':      Control extra rotation offsets for S^2 and S^3 regions.
%     'extra':   Use extra rotation offsets for S^2 and S^3 regions, to try
%                to minimize energy.
%                Sets opt.extra_offset to true.
%     'normal':  Do not use extra offsets
%                Sets opt.extra_offset to false.
%
% Some shortcuts are also provided.
% POPT = PARTITION_OPTIONS(pdefault) just sets POPT to PDEFAULT.
%
% The following are equivalent to PARTITION_OPTIONS(PDEFAULT,'offset','extra'):
% PARTITION_OPTIONS(PDEFAULT,true)
% PARTITION_OPTIONS(PDEFAULT,'extra')
%
% The following are equivalent to PARTITION_OPTIONS(PDEFAULT,'offset','normal'):
% PARTITION_OPTIONS(PDEFAULT,false)
% PARTITION_OPTIONS(PDEFAULT,'normal')
%
%Examples
% > pdefault.extra_offset=false;
% > popt=partition_options(pdefault,'offset','extra')
% popt =
%     extra_offset: 1
%
% > popt=partition_options(pdefault,false)
% popt =
%     extra_offset: 0

% Copyright 2004-2005 Paul Leopardi for the University of New South Wales.
% $Revision 1.10 $ $Date 2005-06-01 $
% Documentation files renamed
% $Revision 1.00 $ $Date 2005-02-12 $
%
% For licensing, see COPYING.
% For references, see AUTHORS.
% For revision history, see CHANGELOG.

popt = pdefault;
nargs = length(varargin);

if nargs == 1
    %
    % Short circuit: single argument is value of extra_offset
    %
    value = varargin{1};
    switch value
    case true
        popt.extra_offset = true;
    case false
        popt.extra_offset = false;
    case 'extra'
        popt.extra_offset = true;
    case 'normal'
        popt.extra_offset = false;
    otherwise
        value_error(value,varargin{:});
    end
    return;
end

nopts = floor(nargs/2);
opt_args = {varargin{1:2:2*nopts-1}};
for k=1:nopts
    if ~ischar([opt_args{k}])
        fprintf('Option names must be character strings\n');
        option_error(varargin{:});
    end
end
opt_vals = {varargin{2:2:2*nopts}};

option_name = 'offset';
pos = strmatch(option_name,opt_args,'exact');
if ~isempty(pos)
    if (length(pos) == 1)
        value = opt_vals{pos};
    else
        duplicate_error(option_name,varargin{:});
    end
    switch value
    case 'extra'
        popt.extra_offset = true;
    case 'normal'
        popt.extra_offset = false;
    otherwise
        value_error(value,varargin{:});
    end
end

function duplicate_error(option_name,varargin)
fprintf('Duplicate option %s\n',option_name);
option_error(varargin{:});
%
% end function

function value_error(value,varargin)
fprintf('Invalid option value ');
disp(value);
option_error(varargin{:});
%
% end function

function option_error(varargin)
fprintf('Error in options:\n');
disp(varargin);
error('Please check "help partition_options" for options');
%
% end function

repo: gspbox


text:
%| argmin(x)
%| argmin(f, x)
%| akin to julia argmin
function out = argmin(f, x)
if isa(f, 'function_handle')
   [~, out] = min(f(x));
else
   [~, out] = min(f);
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function savesnirf(data, outfile,varargin)
%
%    savesnirf(snirfdata, fname)
%       or
%    savesnirf(snirfdata, fname, 'Param1',value1, 'Param2',value2,...)
%
%    Load an HDF5 based SNIRF file, and optionally convert it to a JSON
%    file based on the JSNIRF specification:
%    https://github.com/NeuroJSON/jsnirf
%
%    author: Qianqian Fang (q.fang <at> neu.edu)
%
%    input:
%        snirfdata: a raw SNIRF data, preprocessed SNIRF data or JSNIRF
%             data (root object must be SNIRFData)
%        fname: the output SNIRF (.snirf) or JSNIRF data file name (.jnirs, .bnirs)
%
%    output:
%        data: a MATLAB structure with the grouped data fields
%
%    example:
%        data=loadsnirf('test.snirf');
%        savesnirf(data,'newfile.snirf');
%
%    this file is part of JSNIRF specification: https://github.com/NeuroJSON/jsnirf
%
%    License: GPLv3 or Apache 2.0, see https://github.com/NeuroJSON/jsnirf for details
%

if(nargin<2 || ~ischar(outfile))
    error('you must provide data and a file name');
end

opt=varargin2struct(varargin{:});
if(~isfield(opt,'root'))
    opt.rootname='';
end

if(isfield(data,'SNIRFData'))
    data.nirs=data.SNIRFData;
    data.formatVersion=data.SNIRFData.formatVersion;
    data.nirs=rmfield(data.nirs,'formatVersion');
    data=rmfield(data,'SNIRFData');
end

if(~isempty(outfile))
    if(~isempty(regexp(outfile,'\.[Hh]5$', 'once')))
        saveh5(data,outfile,opt);
    elseif(~isempty(regexp(outfile,'\.[Ss][Nn][Ii][Rr][Ff]$', 'once')))
        data.nirs.data=forceindex(data.nirs.data,'measurementList');
        data.nirs=forceindex(data.nirs,'data');
        data.nirs=forceindex(data.nirs,'stim');
        data.nirs=forceindex(data.nirs,'aux');
        saveh5(data,outfile,opt);
    elseif(~isempty(regexp(outfile,'\.[Jj][Nn][Ii][Rr][Ss]$', 'once'))|| ~isempty(regexp(outfile,'\.[Jj][Ss][Oo][Nn]$', 'once')))
        savejson('SNIRFData',data,'FileName',outfile,opt);
    elseif(regexp(outfile,'\.[Mm][Aa][Tt]$'))
        save(outfile,'data');
    elseif(regexp(outfile,'\.[Bb][Nn][Ii][Rr][Ss]$'))
        savebj('SNIRFData',data,'FileName',outfile,opt);
    else
        error('only support .snirf, .h5, .jnirs, .bnirs and .mat files');
    end
end

% force adding index 1 to the group name for singular struct and cell
function newroot=forceindex(root,name)
newroot=root;
fields=fieldnames(newroot);
idx=find(ismember(fields,name));
if(~isempty(idx) && length(newroot.(name))==1)
    newroot.(sprintf('%s1',name))=newroot.(name);
    newroot=rmfield(newroot,name);
    fields{idx(1)}=sprintf('%s1',name);
    newroot=orderfields(newroot,fields);
end

author: brainstorm-tools


text:
function varargout = drawEllipseCylinder(varargin)
%DRAWELLIPSECYLINDER Draw a cylinder with ellipse cross-section.
%
%   drawEllipseCylinder(CYL)
%   draws the cylinder CYL on the current axis.
%   CYL is a cylinder defined by [x1 y1 z1 x2 y2 z2 r1 r2 roll], with:
%   * [x1 y2 z1] are coordinates of starting point,
%   * [x2 y2 z2] are coordinates of ending point,
%   * R1 and R2 are the lengths of the ellipse semi axes, and
%   * ROLL is the rotation of the cylinder around its main axis (in
%      degrees)
%
%   drawEllipseCylinder(CYL, N)
%   uses N points for discretisation of angle. Default value is 32.
%
%   drawEllipseCylinder(..., OPT)
%   with OPT = 'open' (default) or 'closed', specify if bases of the
%   cylinder should be drawn.
%
%   drawEllipseCylinder(..., 'FaceColor', COLOR)
%   Specifies the color of the cylinder. Any couple of parameters name and
%   value can be given as argument, and will be transfered to the 'surf'
%   matlab function
%
%   H = drawEllipseCylinder(...)
%   returns a handle to the patch representing the cylinder.
%
%
%   Example:
%     figure; drawEllipseCylinder([0 0 0 10 20 30 5 2]);
%
%     figure; drawEllipseCylinder([0 0 0 10 20 30 5 2], 'open');
%
%     figure; drawEllipseCylinder([0 0 0 10 20 30 5 2], 'FaceColor', 'r');
%
%     figure;
%     h = drawEllipseCylinder([0 0 0 10 20 30 5 2]);
%     set(h, 'facecolor', 'b');
%
%     % Draw three mutually intersecting elliptic cylinders
%     p1 = [30 0 0];
%     p2 = [0 30 0];
%     p3 = [0 0 30];
%     radii = [20 10];
%     figure;
%     drawEllipseCylinder([-p1 p1 radii 0], 'FaceColor', 'r');
%     hold on
%     drawEllipseCylinder([-p2 p2 radii 90], 'FaceColor', 'g');
%     drawEllipseCylinder([-p3 p3 radii 90], 'FaceColor', 'b');
%     axis equal
%     set(gcf, 'renderer', 'opengl')
%     view([60 30]); light;
%
%   See also
%     drawCylinder, drawSphere, cylinderMesh, drawLine3d, surf

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2014-02-27
% Copyright 2014-2022 INRA - TPV URPOI - BIA IMASTE

%% Input argument processing

% Check if axes handle is specified
if isAxisHandle(varargin{1})
    hAx = varargin{1};
    varargin(1) = [];
else
    hAx = gca;
end

cyl = varargin{1};
varargin(1) = [];

if iscell(cyl)
    res = zeros(length(cyl), 1);
    for i = 1:length(cyl)
        res(i) = drawEllipseCylinder(hAx, cyl{i}, varargin{:});
    end

    if nargout > 0
        varargout{1} = res;
    end
    return;
end

% default values
N = 32;
closed = true;

% check number of discretization steps
if ~isempty(varargin)
    var = varargin{1};
    if isnumeric(var)
        N = var;
        varargin = varargin(2:end);
    end
end

% check if cylinder must be closed or open
if ~isempty(varargin)
    var = varargin{1};
    if ischar(var)
        if strncmpi(var, 'open', 4)
            closed = false;
            varargin = varargin(2:end);
        elseif strncmpi(var, 'closed', 5)
            closed = true;
            varargin = varargin(2:end);
        end
    end
end

%% Computation of mesh coordinates

% extreme points of cylinder
p1 = cyl(1:3);
p2 = cyl(4:6);

% radius of cylinder
r1 = cyl(7);
r2 = cyl(8);
roll = 0;
if size(cyl, 2) > 8
    roll = cyl(9);
end

% compute orientation angle of cylinder (in degrees)
[theta, phi, rho] = cart2sph2d(p2 - p1);
dphi = linspace(0, 2*pi, N+1);

% generate a cylinder oriented upwards
x = repmat(cos(dphi) * r1, [2 1]);
y = repmat(sin(dphi) * r2, [2 1]);
z = repmat([0 ; rho], [1 length(dphi)]);

% transform points
trans   = localToGlobal3d(p1, theta, phi, roll);
pts     = transformPoint3d([x(:) y(:) z(:)], trans);

% reshape transformed points
x2 = reshape(pts(:,1), size(x));
y2 = reshape(pts(:,2), size(x));
z2 = reshape(pts(:,3), size(x));

%% Display cylinder mesh

% add default drawing options
varargin = [{'FaceColor', 'g', 'edgeColor', 'none'} varargin];

% plot the cylinder as a surface
hSurf = surf(hAx, x2, y2, z2, varargin{:});

% eventually plot the ends of the cylinder
if closed
    ind = find(strcmpi(varargin, 'facecolor'), 1, 'last');
    if isempty(ind)
        color = 'k';
    else
        color = varargin{ind+1};
    end

    patch(hAx, x2(1,:)', y2(1,:)', z2(1,:)', color, 'edgeColor', 'none');
    patch(hAx, x2(2,:)', y2(2,:)', z2(2,:)', color, 'edgeColor', 'none');
end

% format ouptut
if nargout == 1
    varargout{1} = hSurf;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function cmp_feat_caches(dir1, dir2)
% For comparing new feat cache files to known good
% feat cache files after some change

s = dir([dir1 '/*.mat']);

for i = 1:length(s)
  [a,b,c] = fileparts(s(i).name);
  s1 = [dir1 '/' s(i).name];
  s2 = [dir2 '/' b '.mat'];
  if ~exist(s2)
    fprintf('Skipping %s\n', s(i).name);
    continue;
  end

  d1 = load(s1);
  d2 = load(s2);

  %  gt: [2448x1 logical]
  %  overlap: [2448x20 single]
  %  boxes: [2448x4 single]
  %  feat: [2448x9216 single]
  %  class: [2448x1 uint8]

  assert(sum(abs(d1.gt - d2.gt)) == 0);
  assert(sum(sum(abs(d1.overlap - d2.overlap))) == 0);
  assert(sum(sum(abs(d1.boxes - d2.boxes))) == 0);
  assert(sum(sum(abs(d1.feat - d2.feat))) == 0);
  assert(sum(abs(d1.class - d2.class)) == 0);

  fprintf('%d is ok\n', i);
end

author: rbgirshick

repo: rcnn

save_path: github-repos/MATLAB/rbgirshick-rcnn


text:
function [prob] = predictRF(X,RF)

nInst = size(X,1);
trees = RF.Trees; numTrees = numel(trees);
predictions = zeros(nInst,1);
for t = 1:numTrees
    temp = trees{t}.predict(X); temp = str2double(temp);
    predictions = predictions + temp;
end
prob = predictions/numTrees;
end

author: mvallieres

repo: radiomics


text:
%% Importing Tensor Data
%
%%
% Single crystal tensor are imported using the command <tensor.load.html
% tensor.load>. This function automatically detect the file format and
% imports the data. In dependency of the specific format it might be
% necessary to specify the crystal symmetry seperately

% define crystal symmetry
CS = crystalSymmetry('32', [4.916 4.916 5.4054],...
  'X||a*', 'Z||c', 'mineral', 'Quartz');

% define the file name
fname = fullfile(mtexDataPath,'tensor', 'Single_RH_quartz_poly.P');

% import the single crystal tensor
P = tensor.load(fname,CS,'propertyname','piecoelectricity','unit','C/N','DoubleConvention')

%%
% For specific types of tensors, e.g. stiffness tensors there exist
% dedicated import functions that have the form *tensorName.load*

fname = fullfile(mtexDataPath,'tensor','Olivine1997PC.GPa');

cs = crystalSymmetry('mmm',[4.7646 10.2296 5.9942],'mineral','Olivin');

C = stiffnessTensor.load(fname,cs)

author: mtex-toolbox

repo: mtex

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


text:
function [ x, w ] = bashforth_set ( n )

%*****************************************************************************80
%
%% BASHFORTH_SET sets abscissas and weights for Adams-Bashforth quadrature.
%
%  Discussion:
%
%    Adams-Bashforth quadrature formulas are normally used in solving
%    ordinary differential equations, and are not really suitable for
%    general quadrature computations.  However, an Adams-Bashforth formula
%    is equivalent to approximating the integral of F(Y(X)) between X(M)
%    and X(M+1), using an explicit formula that relies only on known values
%    of F(Y(X)) at X(M-N+1) through X(M).  For this reason, the formulas
%    have been included here.
%
%    Suppose the unknown function is denoted by Y(X), with derivative
%    F(Y(X)), and that approximate values of the function are known at a
%    series of X values, which we write as X(1), X(2), ..., X(M).  We write
%    the value Y(X(1)) as Y(1) and so on.
%
%    Then the solution of the ODE Y'=F(X,Y) at the next point X(M+1) is
%    computed by:
%
%      Y(M+1) = Y(M) + Integral ( X(M) < X < X(M+1) ) F(Y(X)) dX
%             = Y(M) + H * Sum ( 1 <= I <= N ) W(I) * F(Y(M+1-I)) approximately.
%
%    In the documentation that follows, we replace F(Y(X)) by F(X).
%
%    The integral:
%
%      Integral ( 0 <= X <= 1 ) F(X) dX.
%
%    The quadrature rule:
%
%      Sum ( 1 <= I <= N ) W(I) * F ( 1 - I ),
%
%    The Adams-Bashforth formulas require equally spaced data.
%
%    Here is how the formula is applied in the case with non-unit spacing:
%
%      Integral ( A <= X <= A+H ) F(X) dX =
%      H * Sum ( 1 <= I <= N ) W(I) * F ( A - (I-1)*H ),
%      approximately.
%
%    The reference lists the second coefficient of the order 8 Adams-Bashforth
%    formula as
%      w(2) =  -1162169.0 / 120960.0
%    but this should be
%      w(2) =  -1152169.0 / 120960.0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964.
%
%    Jean Lapidus, John Seinfeld,
%    Numerical Solution of Ordinary Differential Equations,
%    Academic Press, 1971.
%
%    Daniel Zwillinger, editor,
%    Standard Mathematical Tables and Formulae,
%    30th Edition,
%    CRC Press, 1996.
%
%  Parameters:
%
%    Input, integer N, the order.
%    N should be between 1 and 10, or 12, 14, 16, 18 or 20.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%    W(1) is the weight at X = 0, W(2) the weight at X = -1,
%    and so on.  The weights are rational, and should sum to 1.  Some
%    weights may be negative.
%
  x = zeros ( n, 1 );
  w = zeros ( n, 1 );

  if ( n == 1 )

    w(1) = 1.0;

  elseif ( n == 2 )

    d = 2.0;

    w(1) =   3.0 / d;
    w(2) = - 1.0 / d;

  elseif ( n == 3 )

    d = 12.0;

    w(1) =   23.0 / d;
    w(2) = - 16.0 / d;
    w(3) =    5.0 / d;

  elseif ( n == 4 )

    d = 24.0;

    w(1) =   55.0 / d;
    w(2) = - 59.0 / d;
    w(3) =   37.0 / d;
    w(4) =  - 9.0 / d;

  elseif ( n == 5 )

    d = 720.0;

    w(1) =   1901.0 / d;
    w(2) = - 2774.0 / d;
    w(3) =   2616.0 / d;
    w(4) = - 1274.0 / d;
    w(5) =    251.0 / d;

  elseif ( n == 6 )

    d = 1440.0;

    w(1) =   4277.0 / d;
    w(2) = - 7923.0 / d;
    w(3) =   9982.0 / d;
    w(4) = - 7298.0 / d;
    w(5) =   2877.0 / d;
    w(6) =  - 475.0 / d;

  elseif ( n == 7 )

    d = 60480.0;

    w(1) =    198721.0 / d;
    w(2) =  - 447288.0 / d;
    w(3) =    705549.0 / d;
    w(4) =  - 688256.0 / d;
    w(5) =    407139.0 / d;
    w(6) =  - 134472.0 / d;
    w(7) =     19087.0 / d;

  elseif ( n == 8 )

    d = 120960.0;

    w(1) =     434241.0 / d;
    w(2) =  - 1152169.0 / d;
    w(3) =    2183877.0 / d;
    w(4) =  - 2664477.0 / d;
    w(5) =    2102243.0 / d;
    w(6) =  - 1041723.0 / d;
    w(7) =     295767.0 / d;
    w(8) =    - 36799.0 / d;

  elseif ( n == 9 )

    d = 3628800.0;

    w(1) =   14097247.0 / d;
    w(2) =  -43125206.0 / d;
    w(3) =   95476786.0 / d;
    w(4) = -139855262.0 / d;
    w(5) =  137968480.0 / d;
    w(6) =  -91172642.0 / d;
    w(7) =   38833486.0 / d;
    w(8) =   -9664106.0 / d;
    w(9) =    1070017.0 / d;

  elseif ( n == 10 )

    d = 7257600.0;

    w( 1) =   30277247.0 / d;
    w( 2) = -104995189.0 / d;
    w( 3) =  265932680.0 / d;
    w( 4) = -454661776.0 / d;
    w( 5) =  538363838.0 / d;
    w( 6) = -444772162.0 / d;
    w( 7) =  252618224.0 / d;
    w( 8) =  -94307320.0 / d;
    w( 9) =   20884811.0 / d;
    w(10) =   -2082753.0 / d;

  elseif ( n == 12 )

    d = 958003200.0;

    w( 1) =    4527766399.0 / d;
    w( 2) =  -19433810163.0 / d;
    w( 3) =   61633227185.0 / d;
    w( 4) = -135579356757.0 / d;
    w( 5) =  214139355366.0 / d;
    w( 6) = -247741639374.0 / d;
    w( 7) =  211103573298.0 / d;
    w( 8) = -131365867290.0 / d;
    w( 9) =   58189107627.0 / d;
    w(10) =  -17410248271.0 / d;
    w(11) =    3158642445.0 / d;
    w(12) =    -262747265.0 / d;

  elseif ( n == 14 )

    d = 5230697472000.0;

    w( 1) =    27511554976875.0 / d;
    w( 2) =  -140970750679621.0 / d;
    w( 3) =   537247052515662.0 / d;
    w( 4) = -1445313351681906.0 / d;
    w( 5) =  2854429571790805.0 / d;
    w( 6) = -4246767353305755.0 / d;
    w( 7) =  4825671323488452.0 / d;
    w( 8) = -4204551925534524.0 / d;
    w( 9) =  2793869602879077.0 / d;
    w(10) = -1393306307155755.0 / d;
    w(11) =   505586141196430.0 / d;
    w(12) =  -126174972681906.0 / d;
    w(13) =    19382853593787.0 / d;
    w(14) =    -1382741929621.0 / d;

  elseif ( n == 16 )

    d = 62768369664000.0;

    w( 1) =     362555126427073.0 / d;
    w( 2) =   -2161567671248849.0 / d;
    w( 3) =    9622096909515337.0 / d;
    w( 4) =  -30607373860520569.0 / d;
    w( 5) =   72558117072259733.0 / d;
    w( 6) = -131963191940828581.0 / d;
    w( 7) =  187463140112902893.0 / d;
    w( 8) = -210020588912321949.0 / d;
    w( 9) =  186087544263596643.0 / d;
    w(10) = -129930094104237331.0 / d;
    w(11) =   70724351582843483.0 / d;
    w(12) =  -29417910911251819.0 / d;
    w(13) =    9038571752734087.0 / d;
    w(14) =   -1934443196892599.0 / d;
    w(15) =     257650275915823.0 / d;
    w(16) =     -16088129229375.0 / d;

  elseif ( n == 18 )

    d = 64023737057280000.0;

    w( 1) =     401972381695456831.0 / d;
    w( 2) =   -2735437642844079789.0 / d;
    w( 3) =   13930159965811142228.0 / d;
    w( 4) =  -51150187791975812900.0 / d;
    w( 5) =  141500575026572531760.0 / d;
    w( 6) = -304188128232928718008.0 / d;
    w( 7) =  518600355541383671092.0 / d;
    w( 8) = -710171024091234303204.0 / d;
    w( 9) =  786600875277595877750.0 / d;
    w(10) = -706174326992944287370.0 / d;
    w(11) =  512538584122114046748.0 / d;
    w(12) = -298477260353977522892.0 / d;
    w(13) =  137563142659866897224.0 / d;
    w(14) =  -49070094880794267600.0 / d;
    w(15) =   13071639236569712860.0 / d;
    w(16) =   -2448689255584545196.0 / d;
    w(17) =     287848942064256339.0 / d;
    w(18) =     -15980174332775873.0 / d;

  elseif ( n == 20 )

    d = 102181884343418880000.0;

    w( 1) =      691668239157222107697.0 / d;
    w( 2) =    -5292843584961252933125.0 / d;
    w( 3) =    30349492858024727686755.0 / d;
    w( 4) =  -126346544855927856134295.0 / d;
    w( 5) =   399537307669842150996468.0 / d;
    w( 6) =  -991168450545135070835076.0 / d;
    w( 7) =  1971629028083798845750380.0 / d;
    w( 8) = -3191065388846318679544380.0 / d;
    w( 9) =  4241614331208149947151790.0 / d;
    w(10) = -4654326468801478894406214.0 / d;
    w(11) =  4222756879776354065593786.0 / d;
    w(12) = -3161821089800186539248210.0 / d;
    w(13) =  1943018818982002395655620.0 / d;
    w(14) =  -970350191086531368649620.0 / d;
    w(15) =   387739787034699092364924.0 / d;
    w(16) =  -121059601023985433003532.0 / d;
    w(17) =    28462032496476316665705.0 / d;
    w(18) =    -4740335757093710713245.0 / d;
    w(19) =      498669220956647866875.0 / d;
    w(20) =      -24919383499187492303.0 / d;

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'BASHFORTH_SET - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of N = %d\n', n );
    fprintf ( 1, '  Legal values are 1 through 10,\n' );
    fprintf ( 1, '  or 12, 14, 16, 18 or 20.\n' );
    error ( 'BASHFORTH_SET - Fatal error!' );

  end

  for i = 1 : n
    x(i) = 1 - i;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ n_data, x, fx ] = tran09_values ( n_data )

%*****************************************************************************80
%
%% TRAN09_VALUES returns some values of the order 9 transportation function.
%
%  Discussion:
%
%    The function is defined by:
%
%      TRAN09(x) = Integral ( 0 <= t <= x ) t^9 * exp(t) / ( exp(t) - 1 )^2 dt
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Allan McLeod,
%    Algorithm 757, MISCFUN: A software package to compute uncommon
%    special functions,
%    ACM Transactions on Mathematical Software,
%    Volume 22, Number 3, September 1996, pages 288-301.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 20;

  fx_vec = [ ...
     0.26469772870084897671E-22, ...
     0.11367943653594246210E-12, ...
     0.74428246255329800255E-08, ...
     0.48022728485415366194E-03, ...
     0.11700243014358676725E+00, ...
     0.27648973910899914391E+01, ...
     0.24716631405829192997E+02, ...
     0.12827119828849828583E+03, ...
     0.46842894800662208986E+03, ...
     0.31673967371627895718E+04, ...
     0.46140886546630195390E+04, ...
     0.11952718545392302185E+05, ...
     0.20001612666477027728E+05, ...
     0.31011073271851366554E+05, ...
     0.10352949905541130133E+06, ...
     0.19743173017140591390E+06, ...
     0.33826030414658460679E+06, ...
     0.36179607036750755227E+06, ...
     0.36360622124777561525E+06, ...
     0.36360880558827162725E+06 ];

  x_vec = [ ...
       0.0019531250E+00, ...
       0.0312500000E+00, ...
       0.1250000000E+00, ...
       0.5000000000E+00, ...
       1.0000000000E+00, ...
       1.5000000000E+00, ...
       2.0000000000E+00, ...
       2.5000000000E+00, ...
       3.0000000000E+00, ...
       4.0000000000E+00, ...
       4.2500000000E+00, ...
       5.0000000000E+00, ...
       5.5000000000E+00, ...
       6.0000000000E+00, ...
       8.0000000000E+00, ...
      10.0000000000E+00, ...
      15.0000000000E+00, ...
      20.0000000000E+00, ...
      30.0000000000E+00, ...
      50.0000000000E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
  function hs = jf_add_slider(varargin)
%|function hs = jf_add_slider([options])
%|
%| add a slider control to the current figure to enable interaction
%|
%| option
%|	'callback'	@ (value,data{:})	function to call with
%|						slider value [0,1] and data
%|	'pos'		[l b w h]	slider position
%|	'data'		{}		data to be passed to callback
%|	'value'		0 <= x <= 1	initial slider value, default 0
%|
%| Copyright 2010-1-11, Jeff Fessler, University of Michigan

if ~nargin && ~nargout, ir_usage, end
if nargin && streq(varargin{1}, 'test'), jf_add_slider_test, return, end

arg.callback = @jf_add_slider_call;
arg.pos = [0.1, 0.0, 0.8, 0.03];
arg.data = {};
arg.name = '';
arg.value = 0;
arg.min = [];
arg.max = [];
arg.style = 'slider';
arg.sliderstep = [0.01 0.10];
arg.string = 'button';
arg = vararg_pair(arg, varargin);

args = {};
if ~isempty(arg.min)
	args = {args{:}, 'min', arg.min};
end
if ~isempty(arg.max)
	args = {args{:}, 'max', arg.max};
end

if isempty(arg.callback)
	arg.callback = @jf_add_slider_call_null;
end

% construct gui components
if ir_is_octave % && exist('uicontrol') ~= 2
	warn 'no uicontrol in octave'
	% todo: draw box or something?
else
	hs = uicontrol('style', arg.style, 'string', arg.string, ...
		'units', 'normalized', 'position', arg.pos, ...
		'value', arg.value, ...
		args{:}, ...
		'sliderstep', arg.sliderstep, ...
		'callback', {@jf_add_slider_callback});
end

if ~isempty(arg.name)
	set(gcf, 'name', arg.name) % assign GUI a name for window title
end

if ~nargout
	clear hs
end

% callbacks: these automatically have access to component handles
% and initialized data because they are nested at a lower level.

% slider callback
function jf_add_slider_callback(source, eventdata)
	val = get(source, 'value');
	arg.callback(val, arg.data{:});
end % jf_add_slider_callback

end % jf_add_slider

function jf_add_slider_call(value)
pr value
end % jf_add_slider_call

function jf_add_slider_call_null(value)
end % jf_add_slider_call_null

function jf_add_slider_test
jf_add_slider_test_call(0.5);
if im
	jf_add_slider('callback', @jf_add_slider_test_call, 'data', {gca}, ...
		'value', 0.5)
%	jf_add_slider_test_call(0);
	jf_add_slider('callback', @jf_add_slider_test_call, 'data', {gca}, ...
		'pos', [0.1 0.04 0.8 0.03], ...
		'value', 0, 'style', 'togglebutton', 'string', 'push me')
end
end % jf_add_slider_test

function jf_add_slider_test_call(value, h)
if isvar('h')
	axes(h)
end
if im
	plot([0 value], [0 value], '-')
	axis([0 1 0 1])
end
end % jf_add_slider_test_call


text:
function [A, B] = kafbox_covProd(covfunc, logtheta, x, z);

% covProd - compose a covariance function as the product of other covariance
% functions. This function doesn't actually compute very much on its own, it
% merely does some bookkeeping, and calls other covariance functions to do the
% actual work.
%
% For more help on design of covariance functions, try "help covFunctions".
%
% (C) Copyright 2006 by Carl Edward Rasmussen, 2006-04-06.

for i = 1:length(covfunc)                   % iterate over covariance functions
  f = covfunc(i);
  if iscell(f{:}), f = f{:}; end          % dereference cell array if necessary
  j(i) = cellstr(feval(f{:}));
end

if nargin == 1,                                   % report number of parameters
  A = char(j(1)); for i=2:length(covfunc), A = [A, '+', char(j(i))]; end
  return
end

[n, D] = size(x);

v = [];              % v vector indicates to which covariance parameters belong
for i = 1:length(covfunc), v = [v repmat(i, 1, eval(char(j(i))))]; end

switch nargin
case 3                                              % compute covariance matrix
  A = ones(n, n);                       % allocate space for covariance matrix
  for i = 1:length(covfunc)                   % iteration over factor functions
    f = covfunc(i);
    if iscell(f{:}), f = f{:}; end        % dereference cell array if necessary
    A = A .* feval(f{:}, logtheta(v==i), x);             % multiply covariances
  end

case 4                      % compute derivative matrix or test set covariances
  if nargout == 2                                % compute test set cavariances
    A = ones(size(z,1),1); B = ones(size(x,1),size(z,1));      % allocate space
    for i = 1:length(covfunc)
      f = covfunc(i);
      if iscell(f{:}), f = f{:}; end      % dereference cell array if necessary
      [AA BB] = feval(f{:}, logtheta(v==i), x, z);   % compute test covariances
      A = A .* AA; B = B .* BB;                                % and accumulate
    end
  else                                            % compute derivative matrices
    A = ones(n, n);
    ii = v(z);                                      % which covariance function
    j = sum(v(1:z)==ii);                   % which parameter in that covariance
    for i = 1:length(covfunc)
      f = covfunc(i);
      if iscell(f{:}), f = f{:}; end      % dereference cell array if necessary
      if i == ii
        A = A .* feval(f{:}, logtheta(v==i), x, j);       % multiply derivative
      else
        A = A .* feval(f{:}, logtheta(v==i), x);          % multiply covariance
      end
    end
  end

end

repo: kafbox


text:
function [stlStruct] = import_STL_txt(fileName)

% function [stlStruct] = import_STL_txt(fileName)
% ------------------------------------------------------------------------
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2014/09/25
% 2016/02/24 Updated commenting and header
%------------------------------------------------------------------------

T=txtfile2cell(fileName);

logicSolids= ~cellfun(@isempty,strfind(T,'solid'));
numSolids=nnz(logicSolids)/2;
stlStruct.solidNames=cell(1,numSolids);
stlStruct.solidVertices=cell(1,numSolids);
stlStruct.solidFaces=cell(1,numSolids);
stlStruct.solidNormals=cell(1,numSolids);

for q=1:1:numSolids

    %Find current solid
    indStart=find(logicSolids,1);
    logicSolids(indStart)=0;
    indEnd=find(logicSolids,1);
    logicSolids(indEnd)=0;

    T_solid=T(indStart:indEnd); %Current solid text

    solidName = sscanf(T_solid{1}, 'solid %s'); %Solid name

    %Get vertices
    logicVertices= ~cellfun(@isempty,strfind(T_solid,'vertex'));
    numVertices=nnz(logicVertices);
    T_vertices=T_solid(logicVertices);
    V=cell2mat(cellfun(@(x) sscanf(x,'    vertex %f %f %f')',T_vertices,'UniformOutput',0));

    %Get face normals
    logicNormals= ~cellfun(@isempty,strfind(T_solid,'facet normal'));
    T_normals=T_solid(logicNormals);
    N=cell2mat(cellfun(@(x) sscanf(x,'  facet normal %f %f %f')',T_normals,'UniformOutput',0));

    %Create F
    F=reshape(1:numVertices,[3 numVertices/3])';

    stlStruct.solidNames{q}=solidName;
    stlStruct.solidVertices{q}=V;
    stlStruct.solidFaces{q}=F;
    stlStruct.solidNormals{q}=N;

end

%%
% _*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 outwarp_test(outfun)
%OUTWARP_TEST Test correct implementation of an output warping function.

% Generate random observations
N = randi(50);
y = rand(N,1)*10;

[Noutwarp,info] = outfun('info',y);

% Generate random hyperparameters from plausible box
PLB = info.PLB(:);
PUB = info.PUB(:);
hyp = rand(Noutwarp,1).*(PUB - PLB) + PLB;

hyp

fprintf('---------------------------------------------------------------------------------\n');
fprintf('Check error on inverse of output warping function...\n\n');

sum(abs(y - outfun(hyp,outfun(hyp,y),'inv')))

fprintf('---------------------------------------------------------------------------------\n');
fprintf('Check 1st-order derivative of output warping function...\n\n');

yy = y(randi(N));
derivcheck(@(t) f(t,hyp,outfun),yy);

fprintf('---------------------------------------------------------------------------------\n');
fprintf('Check gradient of output warping function wrt hyperparameters...\n\n');

derivcheck(@(hyp_) f2(yy,hyp_,outfun),hyp);

fprintf('---------------------------------------------------------------------------------\n');
fprintf('Check gradient of derivative of output warping function wrt hyperparameters...\n\n');

derivcheck(@(hyp_) f3(yy,hyp_,outfun),hyp);

end

function [y,dy] = f(t,hyp,outfun)
    [y,dy] = outfun(hyp,t);
end

function [y,dy] = f2(y,hyp,outfun)
    [y,~,dy] = outfun(hyp,y);
end

function [y,dy] = f3(y,hyp,outfun)
    [~,y,~,dy] = outfun(hyp,y);
end

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
function HTotal=uvHessian(xG,lRx,M,includeW)
%%UVHESSIAN Determine the Hessian matrix (second derivative matrix) of a
%           direction cosine measurement with respect to 3D position.
%           Relativity and atmospheric effects are not taken into account.
%
%INPUTS: xG The 3XN target position vectors in the global coordinate system
%           with [x;y;z] components for which gradients are desired.
%       lRx The 3X1  position vector of the receiver. If omitted, the
%           receiver is placed at the origin.
%         M A 3X3 rotation matrix from the global coordinate system to the
%           orientation of the coordinate system at the receiver. If
%           omitted, it is assumed to be the identity matrix.
%  includeW An optional boolean value indicating whether a third direction
%           cosine component should be included. The u and v direction
%           cosines are two parts of a 3D unit vector. Generally, one might
%           assume that the target is in front of the sensor, so the third
%           component would be positive and is not needed. However, the
%           third component can be included if ambiguity exists. The
%           default if this parameter is omitted or an empty matrix is
%           passed is false.
%
%OUTPUTS: HTotal A 3X3X(2+includeW)XN matrix such that for the ith measurement,
%          HTotal(:,:,1,i) is the Hessian matrix with respect to the u
%          component and HTotal(:,:,2,i) is the Hessian matrix with respect
%          to the v component and, if included, HTotal(:,:,3,i) is the
%          Hessian matrix with respect to the w component. The elements in
%          the matrices for each component/ point are
%          ordered [d^2/(dxdx), d^2/(dxdy), d^2/(dxdz);
%                   d^2/(dydx), d^2/(dydy), d^2/(dydz);
%                   d^2/(dzdx), d^2/(dzdy), d^2/(dzdz)];
%          note that each matrix is symmetric (i.e.
%           d^2/(dydx)=d^2/(dxdy) ).
%
%A derivation of the components of the Jacobian is given in [1]. The
%Hessian is just one derivative higher.
%
%EXAMPLE:
%Here, we verify that a numerically differentiated Hessian is consistent
%with the analytic one produced by this function.
% xG=[100;-1000;500];
% lRx=[500;20;-400];
% epsVal=1e-6;
% M=randRotMat(3);
% includeW=true;
%
% H=uvHessian(xG,lRx,M,includeW);
% J=uvGradient(xG,lRx,M,includeW);
% JdX=uvGradient(xG+[epsVal;0;0],lRx,M,includeW);
% JdY=uvGradient(xG+[0;epsVal;0],lRx,M,includeW);
% JdZ=uvGradient(xG+[0;0;epsVal],lRx,M,includeW);
% HNumDiff=zeros(3,3,2);
% for k=1:(2+includeW)
%     HNumDiff(1,1,k)=(JdX(k,1)-J(k,1))/epsVal;
%     HNumDiff(1,2,k)=(JdX(k,2)-J(k,2))/epsVal;
%     HNumDiff(2,1,k)=HNumDiff(1,2,k);
%     HNumDiff(1,3,k)=(JdX(k,3)-J(k,3))/epsVal;
%     HNumDiff(3,1,k)=HNumDiff(1,3,k);
%     HNumDiff(2,2,k)=(JdY(k,2)-J(k,2))/epsVal;
%     HNumDiff(2,3,k)=(JdY(k,3)-J(k,3))/epsVal;
%     HNumDiff(3,2,k)=HNumDiff(2,3,k);
%     HNumDiff(3,3,k)=(JdZ(k,3)-J(k,3))/epsVal;
% end
% max(abs((H(:)-HNumDiff(:))./H(:)))
%The relative error will be on the order of 1e-6 or better, indicating good
%agreement between the numerical Hessian matrix and the actual Hessian
%matrix.
%
%REFERENCES:
%[1] D. F. Crouse, "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%June 2017 David F.Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(includeW))
    includeW=false;
end

if(nargin<3||isempty(M))
    M=eye(3,3);
end

if(nargin<2||isempty(lRx))
    lRx=zeros(3,1);
end

N=size(xG,2);
HTotal=zeros(3,3,2+includeW,N);
for curPoint=1:N
    %Convert the state into the local coordinate system.
    xLocal=M*(xG(1:3)-lRx(1:3));

    x=xLocal(1);
    y=xLocal(2);
    z=xLocal(3);

    x2=x*x;
    y2=y*y;
    z2=z*z;

    r5=norm(xLocal)^5;

    H=zeros(3,3,2);

    %u Hessian values
    %du/(dxdx)
    H(1,1,1)=-((3*x*(y2+z2))/r5);
    %du/(dydy)
    H(2,2,1)=-((x*(x2-2*y2+z2))/r5);
    %du/(dzdz)
    H(3,3,1)=-((x*(x2+y2-2*z2))/r5);
    %du/(dxdy)
    H(1,2,1)=-((y*(-2*x2+y2+z2))/r5);
    %du/(dydx)
    H(2,1,1)=H(1,2,1);
    %du/(dxdz)
    H(1,3,1)=-((z*(-2*x2+y2+z2))/r5);
    %du/(dzdx)
    H(3,1,1)=H(1,3,1);
    %du/(dydz)
    H(2,3,1)=(3*x*y*z)/r5;
    %du/(dzdy)
    H(3,2,1)=H(2,3,1);

    %v Hessian values
    %dv/(dxdx)
    H(1,1,2)=-((y*(-2*x2+y2+z2))/r5);
    %dv/(dydy)
    H(2,2,2)=-((3*y*(x2+z2))/r5);
    %dv/(dzdz)
    H(3,3,2)=-((y*(x2+y2-2*z2))/r5);
    %dv/(dxdy)
    H(1,2,2)=-((x*(x2-2*y2+z2))/r5);
    %dv/(dydx)
    H(2,1,2)=H(1,2,2);
    %dv/(dxdz)
    H(1,3,2)=(3*x*y*z)/r5;
    %dv/(dzdx)
    H(3,1,2)=H(1,3,2);
    %dv/(dydz)
    H(2,3,2)=-((z*(x2-2*y2+z2))/r5);
    %dv/(dzdy)
    H(3,2,2)=H(2,3,2);

    %Rotate the values back into global coordinates
    HTotal(:,:,1,curPoint)=M'*H(:,:,1)*M;
    HTotal(:,:,2,curPoint)=M'*H(:,:,2)*M;

    %w Hessian values
    if(includeW)
        %dw/(dxdx)
        H(1,1,3)=-z*(-2*x2+y2+z2)/r5;
        %dw/(dydy)
        H(2,2,3)=-z*(x2-2*y2+z2)/r5;
        %dw/(dzdz)
        H(3,3,3)=-3*(x2+y2)*z/r5;
        %dw/(dxdy)
        H(1,2,3)=3*x*y*z/r5;
        %dw/(dydx)
        H(2,1,3)=H(1,2,3);
        %dw/(dxdz)
        H(1,3,3)=-x*(x2+y2-2*z2)/r5;
        %dw/(dzdx)
        H(3,1,3)=H(1,3,3);
        %dw/(dydz)
        H(2,3,3)=-y*(x2+y2-2*z2)/r5;
        %dw/(dzdy)
        H(3,2,3)=H(2,3,3);

        HTotal(:,:,3,curPoint)=M'*H(:,:,3)*M;
    end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function indx = r8vec2_sort_heap_index_a ( n, x, y )

%*****************************************************************************80
%
%% R8VEC2_SORT_HEAP_INDEX_A does an indexed heap ascending sort of an R8VEC2.
%
%  Discussion:
%
%    An R8VEC2 is two R8VEC's.
%
%    An R8VEC is a vector of R8 values.
%
%    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.
%
%    ( X(I), Y(I) ) < ( X(J), Y(J) ) if:
%
%    * X(I) < X(J), or
%
%    * X(I) = X(J), and Y(I) < Y(J).
%
%    Once the index array is computed, the sorting can be carried out
%    "implicitly:
%
%      ( X(INDX(1:N)), Y(INDX(1:N) ), is sorted,
%
%    or explicitly, by the call
%
%      call dvec_permute ( n, x, indx )
%      call dvec_permute ( n, y, indx )
%
%    after which ( X(1:N), Y(1:N) ), is sorted.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, real X(N),Y(N), pairs of X, Y coordinates of points.
%
%    Output, integer INDX(N), the sort index.  The
%    I-th element of the sorted array has coordinates ( X(INDX(I)), Y(INDX(I) ).
%
  if ( n < 1 )
    indx = [];
    return
  end

  indx(1:n) = 1 : n;

  if ( n == 1 )
    return
  end

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

  while ( 1 )

    if ( 1 < l )

      l = l - 1;
      indxt = indx(l);
      xval = x(indxt);
      yval = y(indxt);

    else

      indxt = indx(ir);
      xval = x(indxt);
      yval = y(indxt);
      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 ( x(indx(j)) < x(indx(j+1)) || ...
          ( x(indx(j)) == x(indx(j+1)) && y(indx(j)) < y(indx(j+1)) ) )
          j = j + 1;
        end

      end

      if ( xval < x(indx(j)) || ...
          ( xval == x(indx(j)) && yval < y(indx(j)) ) )
        indx(i) = indx(j);
        i = j;
        j = j + j;
      else
        j = ir + 1;
      end

    end

    indx(i) = indxt;

  end

  return
end


text:
function h=smooth_diff(n)
% A smoothed differentiation filter (digital differentiator).
%
% Such a filter has the following advantages:
%
% First, the filter involves both the smoothing operation and differentation operation.
% It can be regarded as a low-pass differention filter (digital differentiator).
% It is well known that the common differentiation operation amplifies the high-frequency noises.
% Therefore, the smoothded differentiation filter would be valuable in experimental (noisy) data processing.
%
% Secondly, the filter coefficients are all convenient integers (simple units) except for an integer scaling factor,
% as may be especially significant in some applications such as those in some single-chip microcomputers
% or digital signal processors.
%
% Usage:
% h=smooth_diff(n)
% n: filter length (positive integer larger no less than 2)
% h: filter coefficients (anti-symmetry)
%
% Examples:
% smooth_demo
%
% Author:
% Jianwen Luo <luojw@bme.tsinghua.edu.cn, luojw@ieee.org> 2004-11-02
% Department of Biomedical Engineering, Department of Electrical Engineering
% Tsinghua University, Beijing 100084, P. R. China
%
% References:
% Usui, S.; Amidror, I.,
% Digital Low-Pass Differentiation for Biological Signal-Processing.
% IEEE Transactions on Biomedical Engineering 1982, 29, (10), 686-693.
% Luo, J. W.; Bai, J.; He, P.; Ying, K.,
% Axial strain calculation using a low-pass digital differentiator in ultrasound elastography.
% IEEE Transactions on Ultrasonics Ferroelectrics and Frequency Control 2004, 51, (9), 1119-1127.

if n>=2 && floor(n)==ceil(n)
    if mod(n,2)==1%is odd
        m=fix((n-1)/2);
        h=[-ones(1,m) 0 ones(1,m)]/m/(m+1);
    else%is even
        m=fix(n/2);
        h=[-ones(1,m) ones(1,m)]/m^2;
    end
else
    error('The input parameter (n) should be a positive integer larger no less than 2.');
end

repo: REST


text:
function kf = gaussian_correlation(xf, yf, sigma)
%GAUSSIAN_CORRELATION Gaussian Kernel at all shifts, i.e. kernel correlation.
%   Evaluates a Gaussian kernel with bandwidth SIGMA for all relative
%   shifts between input images X and Y, which must both be MxN. They must
%   also be periodic (ie., pre-processed with a cosine window). The result
%   is an MxN map of responses.
%
%   Inputs and output are all in the Fourier domain.
%
%   Joao F. Henriques, 2014
%   http://www.isr.uc.pt/~henriques/

	N = size(xf,1) * size(xf,2);
	xx = xf(:)' * xf(:) / N;  %squared norm of x
	yy = yf(:)' * yf(:) / N;  %squared norm of y

	%cross-correlation term in Fourier domain
	xyf = xf .* conj(yf);
	xy = sum(real(ifft2(xyf)), 3);  %to spatial domain

	%calculate gaussian response for all positions, then go back to the
	%Fourier domain
	kf = fft2(exp(-1 / sigma^2 * max(0, (xx + yy - 2 * xy) / numel(xf))));

end

repo: Context-Aware-CF-Tracking


text:
%  smin smooth approximation of min(x,0)
%
%  ::
%
%
%  Args:
%
%  Returns:
%     :
%
%  Note:
%
%  Example:
%
%     See also:
%

author: jmaih


text:
function word_num = s_word_count ( s )

%*****************************************************************************80
%
%% S_WORD_COUNT counts the number of "words" in a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be examined.
%
%    Output, integer WORD_NUM, the number of "words" in the string.
%    Words are presumed to be separated by one or more blanks.
%
  word_num = 0;
  s_length = length ( s );

  if ( s_length <= 0 )
    return;
  end

  blank = 1;

  for i = 1 : s_length

    if ( s(i) == ' ' )
      blank = 1;
    elseif ( blank == 1 )
      word_num = word_num + 1;
      blank = 0;
    end

  end

  return
end


text:
function [col4row,row4col,gain]=linBottleneckAssign(C,solSel)
%%LINBOTTLENECKASSIGN Solve the linear bottleneck assignment problem for a
%          square cost matrix. The problem being solved can be formulated
%          as
%          gain=minimize_phi max_i C_{i,phi(i)} where phi is a permutation
%                                               vector.
%          It can also be formulated as:
%          gain=minimize_{x} max_{i,j} C_{i,j}*x_{i,j}
%          subject to
%          \sum_{j=1}^{n}x_{i,j}=1 for all i
%          \sum_{i=1}^{n}x_{i,j}=1 for all j
%          x_{i,j}=0 or 1.
%
%INPUTS: C An nXn cost matrix that does not contain any NaNs.
%   solSel Multiple optimal solutions typically exist to this problem. This
%          parameter specifies which solution to choose. Possible values
%          are:
%          0 (The default if omitted or an empty matrix is passed) Choose
%            a feasible solution without any particular property.
%          1 Choose the optimal solution that also minimizes
%            sum_i C_{i,phi(i)}.
%          2 Choose the optimal solution that also maximizes
%            sum_i C_{i,phi(i)}.
%
%OUTPUTS: col4row An nX1 vector where the entry in each element is an
%                 assignment of the element in that row to a column. This
%                 is equivalent to phi above.
%         row4col A nX1 vector where the entry in each element is an
%                 assignment of the element in that column to a row.
%            gain The value of the cost function of the linear bottleneck
%                 assignment problem at the optimal point.
%
%This function implements Algorithm 6.1 of Chapter 6.2 of [1]. This is an
%algorithm based on thresholding. The choice of which solution is obtained
%in the end depends on whether one uses maxCardBipartMatching in the final
%step, or whether one modifies the cost matrix in the final step according
%to G (forbidden assignments become either Inf or -Inf) and then uses
%assign2D to solve the problem.
%
%EXAMPLE:
%This is example 6.5 from Chapter 6.2 of [1].
% C=[8,2,3,3;
%    2,7,5,8;
%    0,9,8,4;
%    2,5,6,3];
% [col4row,row4col,gain]=linBottleneckAssign(C)
% %Verify the value of the cost function for the permutation.
% n=size(C,1);
% costVal=-Inf;
% for curRow=1:n
%     cij=C(curRow,col4row(curRow));
%     if(cij>costVal)
%         costVal=cij;
%     end
% end
% costVal
%One will see that gain=5 and costVal agrees with the gain. The cost value
%is five even though the col4row assignment is not identical to the in the
%book. Solution variations depend on variations to the algiruthm used for
%the perfect matching. Here, maxCardBipartMatching is used.
%
%REFERENCES:
%[1] R. Burkard, M. Dell'Amico, and S. Martello,Assignment Problems.
%    Philadelphia: Society for Industrial and Applied Mathematics, 2009.
%
%October 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

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

n=size(C,1);

c0=min(C(:));
c1=max(C(:));

if(c0==c1)
    col4row=(1:n).';
    row4col=(1:n).';
    gain=c0;
    return;
end

sel=C<c1&C>c0;
CSel=C(sel);
if(isempty(CSel))
    %This instance can occur, for example, if the matrix contains only two
    %values, such as if C is the identity matrix.
    cStar=c1;
else
    while(~isempty(CSel))
        cStar=medianHi(CSel(:));
        [c0,c1]=feasCheck(C,cStar,c0,c1);

        sel=C<c1&C>c0;
        CSel=C(sel);
    end
end

if(cStar~=c0)
    [~,c1]=feasCheck(C,c0,c0,c1);
end

G=(C<=c1);
switch(solSel)
    case 0
        [col4row,row4col]=maxCardBipartMatching(G);
    case 1
        CNew=Inf(n,n);
        CNew(G)=C(G);
        [col4row,row4col]=assign2D(CNew,false);
        if(isempty(col4row))%If the cost wouldn't have been finite.
            [col4row,row4col]=maxCardBipartMatching(G);
        end
    case 2
        CNew=-Inf(n,n);
        CNew(G)=C(G);
        [col4row,row4col]=assign2D(CNew,true);
        if(isempty(col4row))%If the cost wouldn't have been finite.
            [col4row,row4col]=maxCardBipartMatching(G);
        end
    otherwise
        error('Unknown solution selected.')
end
gain=c1;

end

function [c0,c1]=feasCheck(C,cStar,c0,c1)

G=(C<=cStar);
col4row=maxCardBipartMatching(G);

if(all(col4row~=0))
    %It is feasible.
    c1=cStar;
else
    %It is infeasible.
    c0=cStar;
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
function STResNet_stage1(varargin)

if ~isempty(gcp('nocreate'))
    delete(gcp)
end

opts.train.gpus = [ 1 ] ;
opts = cnn_setup_environment(opts);

addpath('network_surgery');
opts.dataSet = 'ucf101';
opts.inputdim  = [ 224,  224, 20] ;
initMethod = 'sumAB';
opts.dropOutRatio = 0;
opts.dataDir = fullfile(opts.dataPath, opts.dataSet) ;
opts.splitDir = [opts.dataSet '_splits'];
opts.train.fuseInto = 'spatial'; opts.train.fuseFrom = 'temporal';
opts.train.removeFuseFrom = 0;
opts.backpropFuseFrom = 1;
opts.nSplit = 1 ;
opts.addConv3D = 1 ;
opts.addPool3D  = 0 ;
opts.doSum = 1 ;
opts.doScale = 0 ;
opts.doBiDirectional = 0 ;
opts.nFrames = 5 ;
opts.singleSoftMax = 0;

opts.train.learningRate =  1*[ 1e-3*ones(1,3) 1e-4*ones(1,3)  1e-5*ones(1,3) 1e-6*ones(1,3)]  ;
opts.train.augmentation = 'multiScaleRegular';
opts.train.augmentation = 'f25noCtr';
opts.train.epochFactor = 10 ;
opts.train.batchSize = 128  ;
opts.train.numSubBatches = 32 ;
opts.train.cheapResize = 0 ;
opts.train.fusionLayer = {'res2a_relu', 'res2a'; 'res3a_relu', 'res3a'; 'res4a_relu', 'res4a'; ...
  'res5a_relu', 'res5a'; };

model = ['ST-ResNet50-split=' num2str(opts.nSplit)];
opts.train.numSubBatches =  ceil(opts.train.numSubBatches / max(numel(opts.train.gpus),1));
opts.train.memoryMapFile = fullfile(tempdir, 'ramdisk', ['matconvnet' num2str(opts.nSplit ) '.bin']) ;

opts.modelA = fullfile(opts.modelPath, [opts.dataSet '-img-resnet50-split' num2str(opts.nSplit) '-dr0.mat']) ;
opts.modelB = fullfile(opts.modelPath, [opts.dataSet '-resnet50-split' num2str(opts.nSplit) '-dr0.8.mat']) ;

opts.expDir = fullfile(opts.dataDir, [opts.dataSet '-' model]) ;

opts.train.startEpoch = 1 ;
opts.train.epochStep = 1 ;
opts.train.numEpochs = 8 ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.imdbPath = fullfile(opts.dataDir, [opts.dataSet '_split' num2str(opts.nSplit) 'imdb.mat']);

opts.train.saveAllPredScores = 1;
opts.train.denseEval = 1;
opts.train.plotDiagnostics = 0 ;
opts.train.continue = 1 ;
opts.train.prefetch = 1 ;
opts.train.expDir = opts.expDir ;
opts.train.numAugments = 1;
opts.train.frameSample = 'random';
opts.train.nFramesPerVid = 1;

opts = vl_argparse(opts, varargin) ;
% -------------------------------------------------------------------------
%                                                   Database initialization
% -------------------------------------------------------------------------
if exist(opts.imdbPath)
  imdb = load(opts.imdbPath) ;
  imdb.flowDir = opts.flowDir;
else
  imdb = cnn_setup_data(opts) ;
  save(opts.imdbPath, '-struct', 'imdb', '-v6') ;
end
nClasses = length(imdb.classes.name);

% -------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------
for model = {opts.modelA opts.modelB}
  model = model{:};
  if ~exist(model,'file')
    [~, baseModel] = fileparts(model);
    fprintf('Downloading base model file: %s ...\n', baseModel);
    mkdir(fileparts(model)) ;
    urlwrite(...
    ['http://ftp.tugraz.at/pub/feichtenhofer/st-res/ts-base/' baseModel '.mat'], ...
      model) ;
  end
end
netA = load(opts.modelA) ;
netB = load(opts.modelB) ;
if isfield(netA, 'net'), netA=netA.net;end
if isfield(netB, 'net'), netB=netB.net;end

if ~isfield(netA, 'meta')
  netA = update_legacy_net(netA);
  netA = dagnn.DagNN.fromSimpleNN(netA) ;
  netA = netA.saveobj() ;
end
if ~isfield(netB, 'meta'),    netB = dagnn.DagNN.fromSimpleNN(netB) ; netB = netB.saveobj() ; end;

if ~isfield(opts.train, 'fusionLayer');
  sum_layersA = find(strcmp({netA.layers(:).type},'dagnn.Sum')) + 1 ; % +1 for relu
  sum_layersB = find(strcmp({netB.layers(:).type},'dagnn.Sum')) + 1 ;
  opts.train.fusionLayer = {netA.layers(sum_layersA).name; netB.layers(sum_layersB).name};
  opts.train.fusionLayer = reshape({opts.train.fusionLayer{:}},2,[])';
end

f = find(strcmp({netA.layers(:).type}, 'dagnn.Loss'));
netA.layers(f(1)-1).name = 'prediction';
f = find(strcmp({netB.layers(:).type}, 'dagnn.Loss'));
netB.layers(f(1)-1).name = 'prediction';

fusionLayerA = []; fusionLayerB = [];
if ~isempty(opts.train.fusionLayer)
  for i=1:numel(netA.layers)
   if isfield(netA.layers(i),'name') && any(strcmp(netA.layers(i).name,opts.train.fusionLayer(:,1)))
     fusionLayerA = [fusionLayerA i];
   end
  end
  for i=1:numel(netB.layers)
   if  isfield(netB.layers(i),'name') && any(strcmp(netB.layers(i).name,opts.train.fusionLayer(:,2)))
     fusionLayerB = [fusionLayerB i];
   end
  end
end

netNormalization = netB.meta.normalization;
netA.meta.normalization.averageImage = mean(mean(netA.meta.normalization.averageImage, 1), 2);
netB.meta.normalization.averageImage = mean(mean(netB.meta.normalization.averageImage, 1), 2);

netB.meta.normalization.averageImage = gather(cat(3,netB.meta.normalization.averageImage, netA.meta.normalization.averageImage));

% rename layers, params and vars
for x=1:numel(netA.layers)
  if isfield(netA.layers(x), 'name'), netA.layers(x).name = [netA.layers(x).name '_spatial'] ;  end
end
for x=1:numel(netB.layers)
  if isfield(netB.layers(x), 'name'), netB.layers(x).name = [netB.layers(x).name '_temporal']; end
end

netA =  dagnn.DagNN.loadobj(netA);
for i = 1:numel(netA.vars),  if~strcmp(netA.vars(i).name,'label'), netA.renameVar(netA.vars(i).name, [netA.vars(i).name '_spatial']); end; end;
for i = 1:numel(netA.params),  netA.renameParam(netA.params(i).name, [netA.params(i).name '_spatial']); end;

netB =  dagnn.DagNN.loadobj(netB);
for i = 1:numel(netB.vars), if~strcmp(netB.vars(i).name,'label'), netB.renameVar(netB.vars(i).name, [netB.vars(i).name '_temporal']); end;end;
for i = 1:numel(netB.params),  netB.renameParam(netB.params(i).name, [netB.params(i).name '_temporal']); end;
if opts.addConv3D & any(~cellfun(@isempty,(strfind(opts.train.fusionLayer, 'prediction'))))
  if strcmp(opts.train.fuseInto,'temporal')
    [ netB ] = insert_conv_layers( netB, fusionLayerB(end), 'initMethod', initMethod, 'batchNormalization', bnormFusion, 'dropOutRatio', 0  );
  else
    [ netA ] = insert_conv_layers( netA, fusionLayerA(end), 'initMethod', initMethod, 'batchNormalization', bnormFusion, 'dropOutRatio', 0  );
  end
end
if ~opts.addConv3D && ~opts.doSum
  if strcmp(opts.train.fuseInto,'temporal')
    [ netB ] = insert_conv_layers( netB, fusionLayerB, 'initMethod', initMethod, 'batchNormalization', bnormFusion, 'dropOutRatio', 0  );
  else
    [ netA ] = insert_conv_layers( netA, fusionLayerA, 'initMethod', initMethod, 'batchNormalization', bnormFusion, 'dropOutRatio', 0  );
  end
end

if opts.train.removeFuseFrom,
  switch opts.train.fuseFrom
    case 'spatial'
      netA.layers = netA.layers(1:fusionLayerA(end)); netA.rebuild;
    case'temporal'
      netB.layers = netB.layers(1:fusionLayerB(end)); netB.rebuild;
  end
end

netA = netA.saveobj() ;
netB = netB.saveobj() ;

net.layers = [netA.layers netB.layers] ;
net.params =  [netA.params netB.params] ;
net.meta = netB.meta;
net = dagnn.DagNN.loadobj(net);

net = dagnn.DagNN.setLrWd(net, 'convFiltersLRWD', [1 1], 'convBiasesLRWD', [2 0], ...
  'fusionFiltersLRWD', [1 1], 'fusionBiasesLRWD', [2 0], ...
  'filtersLRWD' , [1 1], 'biasesLRWD' , [2 0] ) ;

clear netA netB;
if opts.doBiDirectional, nFuse = 2; else nFuse=1; end;
for s=1:nFuse
  if s==2
    [opts.train.fuseInto, opts.train.fuseFrom] = deal(opts.train.fuseFrom, opts.train.fuseInto);
  end
  for i = 1:size(opts.train.fusionLayer,1)
    if strcmp(opts.train.fuseInto,'spatial')
      i_fusion = find(~cellfun('isempty', strfind({net.layers.name}, ...
        [opts.train.fusionLayer{i,1} '_' opts.train.fuseInto])));
    else
        i_fusion = find(~cellfun('isempty', strfind({net.layers.name}, ...
        [opts.train.fusionLayer{i,2} '_' opts.train.fuseInto])));
    end

    if opts.doSum
      inputVars =  net.layers(strcmp({net.layers.name},[opts.train.fusionLayer{i,1}, ...
          '_' opts.train.fuseFrom])).outputs;

      if opts.doScale
        name = [net.layers(i_fusion(end)).name '_scale'];

        block = dagnn.Scale() ;
        for j = i_fusion(end):-1:1
          p_from = net.layers(j).params;
          if numel(p_from) == 1 ||  numel(p_from) == 2, break; end;
        end

        block.size = [1 1 size(net.params(net.getParamIndex(p_from{1})).value,4)];

        pars = {[name '_f'], [name '_b']} ;
        net.addLayerAt(i_fusion(end), name, ...
          block, ...
          inputVars, ...
          name, ...
          pars) ;
         params = block.initParams(0.1);
        [net.params(net.getParamIndex(pars)).value] = deal(params{:}) ;

        inputVars = {name};
        i_fusion(end) = i_fusion(end)+1;
      end

      name = [net.layers(i_fusion(end)).name '_sum'];

      block = dagnn.Sum() ;

      inputVars = [inputVars{:} net.layers(strcmp({net.layers.name},[opts.train.fusionLayer{i,2}, ...
          '_' opts.train.fuseInto])).outputs];
      net.addLayerAt(i_fusion(end), name, block, ...
                 inputVars, ...
                  name) ;
        % next line is for disabling ReLUshortcut
        net.layers(cellfun(@(x) any(isequal(x,net.getVarIndex(inputVars(end)))), {net.layers.inputIndexes} )).block.useShortCircuit = 0;

    else
      name = [net.layers(i_fusion(end)).name '_concat'];
      block = dagnn.Concat() ;
      net.addLayerAt(i_fusion(end), name, block, ...
                 [net.layers(strcmp({net.layers.name},[opts.train.fusionLayer{i,1} '_spatial'])).outputs ...
                  net.layers(strcmp({net.layers.name},[opts.train.fusionLayer{i,2} '_temporal'])).outputs], ...
                  name) ;
    end

    net.layers(i_fusion(end)+2).inputs{1} = name;

  end
end

net.addVar('input_flow')
net.vars(net.getVarIndex('input_flow')).fanout = net.vars(net.getVarIndex('input_flow')).fanout + 1 ;

% set input for conv layers
i_conv1= find(~cellfun('isempty', strfind({net.layers.name},'temporal')));
net.layers(i_conv1(1)).inputs = {'input_flow'};

net.renameVar(net.vars(1).name, 'input');

if opts.addConv3D
    conv_layers = find(arrayfun(@(x) isa(x.block,'dagnn.Conv') && isequal(x.block.size(1),1) ...
      && (~isempty(strfind(x.name,'b_branch2a')) || ~isempty(strfind(x.name,'b_branch2c'))) && isequal(x.block.stride(1),1), net.layers ))

    for l=conv_layers
       disp(['converting ' net.layers(l).name ' to ConvTime'])
       block = dagnn.ConvTime() ;   block.net = net ;

       kernel = net.params(net.getParamIndex(net.layers(l).params{1})).value;
       sz = size(kernel);
       kernel = cat(2, kernel/3, kernel/3, kernel/3  );

       net.params(net.getParamIndex(net.layers(l).params{1})).value = kernel;
       pads = size(kernel); pads = ceil(pads(1:2) / 2) - 1;

       block.pad = [pads(1),pads(1), pads(2),pads(2)] ;
       block.size = size(kernel);
       block.hasBias = net.layers(l).block.hasBias;

       block.stride(1) =  net.layers(l).block.stride(1)*net.layers(l).block.stride(2);
       net.layers(l).block = block;
    end
end
if opts.addPool3D

  poolLayers = {'pool5'; };

  for j=1:numel(poolLayers)
    for s = {'spatial', 'temporal'}
      i_pool = find(strcmp({net.layers.name},[poolLayers{j} '_' char(s)]));
      block = dagnn.PoolTime() ;
      block.poolSize = [1 opts.nFrames ];
      block.pad = [0 0 0 0];
      block.stride = [1 1];
      block.method = 'avg';
      nPool = 1 ;
      if opts.addPool3D == 2
        block.method = 'max';
      end

      name = {};
      outputVars = {};
      inputVars = net.layers(i_pool).outputs;
      for k=1:nPool
        name{end+1} = [poolLayers{j} '_pool_' char(s) num2str(k)];
        disp(['injecting ' name{end} ' as PoolTime'])
        net.addLayerAt(i_pool, name{end}, block, ...
                      [inputVars], {name{end}}) ;
      end

      % chain input of l that has layer as input
      for l = 1:numel(net.layers)
          if ~strcmp(net.layers(l).name, name)
            sel = find(strcmp(net.layers(l).inputs, net.layers(i_pool).outputs{1})) ;
            if any(sel)
                net.layers(l).inputs{sel} = inputVars;
            end
          end
      end

    end
  end
end % add pool3d

net = dagnn.DagNN.insertLossLayers(net, 'numClasses', nClasses) ;

opts.train.augmentation = 'f25noCtr';
opts.train.frameSample = 'temporalStrideRandom';
opts.train.nFramesPerVid = opts.nFrames * 1;
opts.train.temporalStride = 5:15;

opts.train.valmode = 'temporalStrideRandom';
opts.train.numValFrames = 25 ;
opts.train.saveAllPredScores = 1 ;
opts.train.denseEval = 1;
opts.train.opts.nFramestack = opts.nFrames;

net.meta.normalization.rgbVariance = [];
opts.train.train = find(ismember(imdb.images.set, [1])) ;
opts.train.train = repmat(opts.train.train,1,opts.train.epochFactor);

% opts.train.train = NaN; % for testing only
opts.train.denseEval = 1;

%%
zero_drs =  find(arrayfun(@(x) isa(x.block,'dagnn.DropOut') && x.block.rate == 0, net.layers)) ;
net.removeLayer({net.layers(zero_drs).name}, true);

if ~isnan(opts.dropOutRatio)
  dr_layers = find(arrayfun(@(x) isa(x.block,'dagnn.DropOut'), net.layers)) ;
  if opts.dropOutRatio > 0
    net.layers(dr_layers).block.rate = opts.dropOutRatio;
  else
    net.removeLayer({net.layers(dr_layers).name}, true);
  end
end

if opts.singleSoftMax
  pred_layers = [];
  for l=1:numel(net.layers)
    if isempty( net.layers(l).params ), continue; end;
    if size(net.params(net.getParamIndex(net.layers(l).params{1})).value,4) == nClasses || ...
        size(net.params(net.getParamIndex(net.layers(l).params{1})).value,5) == nClasses % 3D FC layer
          pred_layers = [pred_layers l];
          net.vars(net.layers(l).outputIndexes).precious = 1;
    end
  end
  pred_layers = fliplr(pred_layers) ; % remove the spatial layer
  paramsIdx1 = net.getParamIndex(net.layers(pred_layers(1)).params);
  paramsIdx2 = net.getParamIndex(net.layers(pred_layers(2)).params);
  for p = 1:numel(paramsIdx1)
    sz = size(net.params(paramsIdx1(p)).value);
    if numel(sz) > 2
      net.params(paramsIdx1(p)).value = cat(3,net.params(paramsIdx1(p)).value, net.params(paramsIdx2(p)).value);
    else
      net.params(paramsIdx1(p)).value = net.params(paramsIdx1(p)).value + net.params(paramsIdx2(p)).value;
    end
  end
   block = dagnn.Concat() ;
   newName = ['opts.singleSoftMaxConcat'];
   net.addLayer(newName, block, ...
                 [net.layers(pred_layers).inputs], ...
                 newName) ;
  net.layers(pred_layers(1)).inputs =    newName ;
  % remove layers of the second prediction path
  for l = numel(net.layers):-1:1
    for f = net.layers(pred_layers(2)).outputs
       sel = find(strcmp(f, net.layers(l).inputs )) ;
       if ~isempty(sel)
         fprintf('removing ayer %s \n', net.layers(l).name);
         net.removeLayer({net.layers(l).name});
       end
    end
  end
  net.removeLayer({net.layers(pred_layers(2)).name});
end

net.layers(~cellfun('isempty', strfind({net.layers(:).name}, 'err'))) = [] ;
net.rebuild() ;
opts.train.derOutputs = {} ;
for l=1:numel(net.layers)
  if isa(net.layers(l).block, 'dagnn.Loss') && isempty(strfind(net.layers(l).block.loss, 'err'))
    if opts.backpropFuseFrom || ~isempty(strfind(net.layers(l).name, opts.train.fuseInto ))
      fprintf('setting derivative for layer %s \n', net.layers(l).name);
      opts.train.derOutputs = [opts.train.derOutputs, net.layers(l).outputs, {1}] ;
    end
    net.addLayer(['err1_' net.layers(l).name(end-7:end) ], dagnn.Loss('loss', 'classerror'), ...
             net.layers(l).inputs, 'error') ;
  end
end

net.conserveMemory = 1 ;
fn = getBatchWrapper_rgbflow(net.meta.normalization, opts.numFetchThreads, opts.train) ;
[info] = cnn_train_dag(net, imdb, fn, opts.train) ;

author: feichtenhofer

repo: st-resnet


text:
%PLANESTRAIN Plane strain physics mode.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function fem1d_heat_explicit_test ( )

%*****************************************************************************80
%
%% FEM1D_HEAT_EXPLICIT_TEST tests the FEM1D_HEAT_EXPLICIT library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_HEAT_EXPLICIT_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the FEM1D_HEAT_EXPLICIT library.\n' );

  fem1d_heat_explicit_test01 ( );
  fem1d_heat_explicit_test02 ( );
  fem1d_heat_explicit_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_HEAT_EXPLICIT_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function y=fshift(x,f,fs)

% y=fshift(x,f,fs)
%
% Shift time domain signal x, sampled at a sampling rate of fs by
% f. fs defaults to 2.

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% 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/>.

error(nargchk(2,3,nargin));
error(chk_param(x,'x','vector','horizontal'));
error(chk_param(f,'f','scalar','real'));
if (nargin<3)
  fs=2;
end
error(chk_param(fs,'fs','scalar','real','>',0));

y=x.*exp((j*pi*f/(fs/2))*(0:length(x)-1));

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
addpath(genpath('D:\GitHub\KiloSort2')) % path to kilosort folder
addpath('D:\GitHub\npy-matlab')

pathToYourConfigFile = 'D:\GitHub\KiloSort2\configFiles'; % take from Github folder and put it somewhere else (together with the master_file)
run(fullfile(pathToYourConfigFile, 'configFileBench384.m'))

% common options for every probe
% ops.chanMap     = 'D:\GitHub\KiloSort2\configFiles\neuropixPhase3A_kilosortChanMap_385.mat';
ops.chanMap = 'H:\DATA\Spikes\Bekesy\neuropixPhase3B1_kilosortChanMap.mat';
ops.trange      = [0 Inf]; % TIME RANGE IN SECONDS TO PROCESS

 % these settings overwrite any settings from config
ops.Th       = [10 10];     % threshold on projections (like in Kilosort1)
ops.lam      = 40^2;   % weighting on the amplitude penalty (like in Kilosort1, but it has to be much larger)

ops.ThS      = [8 8];  % lower bound on acceptable single spike quality
ops.momentum = [20 400]; % number of samples to average over
ops.minFR    = 1/50; % minimum spike rate (Hz)

ops.sigmaMask  = 30;

ops.Nfilt       = 1024; % max number of clusters
ops.nfullpasses = 3; % how many forward backward passes to do
ops.nPCs        = 3; % how many PCs to project the spikes into

ops.useRAM = 0;
ops.ccsplit = .99;
ops.NchanTOT = 385;
ops.nPCs = 3;

% rootZ = 'D:\DATA\ALLEN\mouse366119\probeC_2018-03-02_15-18-32_SN619041624\experiment1\recording1\continuous\Neuropix-120.0\';
% rootZ = 'H:\DATA\Spikes\Robbins\ZNP1';
rootZ = 'H:\DATA\Spikes\Bekesy';
rootH = 'H:\DATA\Spikes\temp\';

fs = dir(fullfile(rootZ, '*.bin'));

fname = fs(1).name;

ops.fbinary     = fullfile(rootZ,  fname);
ops.fproc       = fullfile(rootH, 'temp_wh.dat'); % residual from RAM of preprocessed data

%%
% preprocess data
rez = preprocessDataSub(ops);

rez.ops.NchanTOT

rez.iorig = 1:rez.temp.Nbatch;

% fname = fullfile(rootZ, 'rez.mat');
% save(fname, 'rez');
%%

% clusterSingleBatches;
clusterSingleBatches;

% figure(191);
% imagesc(rez.ccb(rez.iorig, rez.iorig), [20 100])
%%
rez = learnAndSolve8(rez);

rez2    = splitAllClusters(rez);

rezToPhy(rez2, rootZ);

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function z = ne(x,y)
%NE Not equal (~=) for sptensors.
%
%   A ~= B compares the elements of A and B for equality. The arguments can
%   be a pair of sptensors, an sptensor and a tensor, or an sptensor and a
%   scalar.  Regardless, the result is always returned as a sparse tensor.
%
%   See also SPTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

%% Observations for sparse matrix case.
% The result of a ~= 5 is sparse.
% The result of a ~= 0 is sparse.
% The result of a ~= full(a) is sparse.

%% Case 1: One argument is a scalar
if isscalar(y)
    if y == 0
        z = sptensor(x.subs,true,size(x));
    else
        subs1 = x.subs(x.vals ~= y,:);
        subs2 = setdiff(allsubs(x),x.subs,'rows');
        z = sptensor([subs1;subs2],true,size(x));
    end
    return;
end

% Call back with the arguments reversed.
if isscalar(x)
    z = ne(y,x);
    return;
end

%% Case 2: Both x and y are tensors or some sort
% Check that the sizes match
if ~isequal(x.size,y.size)
    error('Size mismatch');
end

% Case 2a: Two sparse tensors
if isa(x,'sptensor') && isa(y,'sptensor')

    % find entries where either x *or* y is nonzero, but not both
    subs1 = setxor(x.subs,y.subs,'rows');
    % find entries where both are nonzero, but inequal
    subs2 = intersect(x.subs,y.subs,'rows');
    subs2 = subs2(extract(x,subs2) ~= extract(y,subs2),:);
    % put it all together
    z = sptensor([subs1;subs2],true,size(x));
    return;

end

% Case 2b: y is a dense tensor
if isa(y,'tensor')

    % find entries where x is zero but y is nonzero
    subs1 = setdiff(allsubs(x),union(x.subs,find(y == 0),'rows'),'rows');

    % find entries where x is nonzero but not equal to y
    subs2 = x.subs(x.vals ~= y(x.subs,'extract'),:);

    % put it all together
    z = sptensor([subs1;subs2],true,size(x));
    return;

end

%% Otherwise
error('The arguments must be two sptensors or an sptensor and a scalar.');

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
% Gauss-Seidel 迭代法（分量形式）求解方程组
clear;
% 输入值
A = [10, -1, -2; -1, 10, -2; -1, -1, 5];
b = [7.2; 8.3; 4.2];
tol = 1e-5;
N = 100;
x = [0; 0; 0];
x_backup = [0; 0; 0];
y = [0; 0; 0];
%
A_ = A;
for i = 1 : length(A)
    A_(i,i) = 0;
end

for i = 0 : N
    for j = 1 : length(A)
        y(j,1) = (b(j) - A_(j,:) * x) / A(j,j);
        x_backup(j) = x(j);   % 备份“老值”
        x(j) = y(j);     % “新值” 替换 “老值”
    end
    if (max(abs(x_backup - y)) < tol)
        fprintf('迭代次数: %d\n', i);
        fprintf('方程组的根: %10.8f\n', y);
        break;
    end
end
if i == N
    fprintf('迭代方法失败\n');
end

repo: NumericalAnalysis


text:
% function fval = f11(x)
%
% Bound = [-65.536 65.536];
%
% if nargin==0
%     fval = Bound;
% else
%     data = [-32 -16   0  16  32 ...
%             -32 -16   0  16  32 ...
%             -32 -16   0  16  32 ...
%             -32 -16   0  16  32 ...
%             -32 -16   0  16  32;
%             -32 -32 -32 -32 -32 ...
%             -16 -16 -16 -16 -16 ...
%               0   0   0   0   0 ...
%              16  16  16  16  16 ...
%              32  32  32  32  32];
%     fval = 1/(1/500+sum(1./((x(1)-data(1,:)).^6+(x(2)-data(2,:)).^6+(1:25))));
% end
function fval =f11(x)

Bound=[-600 600];

if nargin==0
   fval = Bound;
else
    [Dim, PopSize] = size(x);
    indices = repmat(1:Dim, PopSize, 1);
    fval  = .00025*sum((x-100).^2) - prod( cos((x-100)./sqrt(indices')) ) +1;
end

author: sfvsfv

repo: Mathematical-modeling


text:
function [g,B,mu]=BaylissTapering(sidelobedB,N,xyPoints,a)
%%BAYLISSTAPERING The Bayliss tapering is a set of complex amplitude
%          weights for a continuous circular (narrowband) aperture antenna
%          that will form a difference beam (odd symmetry about an axis)
%          and hold the four closest sidelobes to a desired level. Such a
%          tapering can be discretized and applied to the elements in a
%          circular phased array (An array of antenna elements can be
%          viewed as a discrete approximation to a continuous aperture).
%          This function will provide the tapering values at a set of
%          discrete points given by xyPoints (the origin is taken to be the
%          center of the aperture). The radius of the aperture can either
%          be provided or is taken as value of the farthest point provided.
%          This function also returns coefficients for one to efficiently
%          compute the tapering at points on their own. If xyPoints is
%          empty, only the coefficients are returned. The difference axis
%          generated by this function is the y-axis.
%
%INPUTS: sidelobedB The number of decibels of the ratio of the close-in
%             sidelobe voltages to the main lobe voltage. This must be a
%             negative number. A typical value is -30.
%           N The Bayliss tapering is computed using a certain number of
%             terms. Using too many terms can be undesirable as edge
%             illumination increases, as noted in [1]. If this parameter is
%             omitted or an empty matrix is passed, then the default of 17
%             is used. In [1], it is suggested that N be chosen to be
%             <2*a/lambda, where a is the radius of the aperture and
%             lambda the wavelength.
%    xyPoints A 2XnumPoints set of numPoints points in the aperture plane
%             at which the tapering values should be evaluated. If this
%             parameter is omitted or an empty matrix is passed, then an
%             empty matrix is returned for the output g. The center of the
%             aperture is taken to be the origin.
%           a The radius of the aperture. Tapering weights for points in
%             xyPoints outside of the aperture are taken to be zero. If
%             this parameter is omitted or an empty matrix is passed, then
%             the radius is taken to be the distance of the farthest point
%             from the origin in xyPoints.
%
%OUTPUTS: g The NX1 set of discretized Bayliss tapering values evaluated at
%           the points given in xyPoints. If xyPoints is omitted, then this
%           is an empty matrix. All Bayliss tapering values are imaginary.
%           The values are not normalized.
%     B, mu These two outputs can be used to evaluate the bayliss tapering
%           values at arbitrary points. B is an NX1 vector and mu is an
%           (N+1)X1 vector. Given a 2X1 point xy and a radius of the
%           aperture, set the normalized radius to p=pi*norm(xy)/a and
%           the Bayliss tapering weight g at the point is
%           g=(xy(1)/norm(xy))*sum(B.*besselj(1,mu(1:N)*p));
%
%This function implements the algorithm of [1] using the polynomial
%interpolation values in the table below Figure 4. This approximation means
%that low sidelobe patterns (-45 dB and below) will not have good fidelity
%sidelobes.
%
%EXAMPLE 1:
%Here, we evaluate the tapering values for 30dB down on a fine grid of
%points to plot what the imaginary part of the tapering weights looks
%like. The real part is all zero.
% numPoints=300;
% points1D=linspace(-1,1,numPoints);
% [X,Y]=meshgrid(points1D,points1D);
% %The Bayliss tapering weights, evaluated across the aperture. Points
% %outside the aperture are assigned a weight of 0. All Bayliss weights are
% %imaginary.
% xyPoints=[X(:)';Y(:)'];
% a=1;%Aperture radius=1.
% gBayliss=BaylissTapering(-30,17,xyPoints,a);
% gBayliss=reshape(gBayliss,numPoints,numPoints);
%
% figure(1)
% clf
% surface(X,Y,imag(gBayliss),'EdgeColor','None')
% colormap(jet(256))
% colorbar()
% view(45,45)
% light()
% axis square
% h1=xlabel('x');
% h2=ylabel('y');
% title('Bayliss Tapering Weight')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%EXAMPLE 2:
%Here, we consider the array response when using tapering values for 30dB
%sidelobes on a circular array with lambda/2 spacing between elements.
%First, we create a circular array. The element locations are given in
%terms of the wavelength lambda, so lambda will not appear in the
%equations for the sum beam.
% %First, we create a circular array. The element locations are given in
% %terms of the wavelength lambda, so lambda will not appear in the
% %equations for the sum beam.
% xyVals=getShaped2DLattice([25;25],'circular');
% %Get the tapering. It is -30dB and  nBar=17;
% N=17;
% sidelobedB=-30;
% g=BaylissTapering(sidelobedB,N,xyVals);
%
% %The tapering matrix
% T=diag(g);
%
% %Now, display the response with the tapering. We normalize it with
% %respect to the peak value and plot the result in decibels (power).
% [Rsp,U,V]=standardUVBeamPattern(T,xyVals,'NormPowGain');
%
% figure(1)
% clf
% surface(U,V,10*log10(Rsp),'EdgeColor','None')
% colormap(jet(256));
% caxis([-40,0])
% colorbar()
% view(45,30)
% light()
% axis square
% h1=xlabel('u');
% h2=ylabel('v');
% h3=zlabel('Response, Decibels');
% title('Bayliss Weighted Array Response')
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h3,'FontSize',14,'FontWeight','bold','FontName','Times')
%
%Note that the difference pattern produced has a line of symmetry about
%the y axis. To flip the symmetry axis (e.g. for a vertical difference
%beam), then simply use g=BaylissTapering(sidelobedB,N,flipud(xyVals));
%
%REFERENCES:
%[1] E. T. Bayliss, "Design of monopulse antenna difference patterns with
%    low sidelobes," The Bell System Technical Journal, vol. 47, no. 5, pp.
%    623-650, May-Jun. 1968.
%
%August 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(N))
   N=17;
end

%The definition og the mu_m terms in Equation 6 in [1]. The first zero is
%the one indexed by 0 in the paper. Thus, this goes from mu_0 to mu_N.
mu=BesselJDerivZeros(1,N+1)/pi;

%This holds the coefficients for the interpolating polynomials given below
%Figure 4 in [1]. The polynomials all take the desired sidelobe level in
%decibels as an input parameter. The first row is for the term A, which
%is a translation of the SNR parameter to a parameter in the paper. The
%next four rows are xi_1 to xi_4, which are the locations of the first four
%zeroes in the modified pattern. The final row is for p_0, which is related
%to the point at which the peak of the asymptotic difference pattern=1.
polyCoeffTable=[0.30387530,-0.05042922,-0.00027989,-0.00000343,-0.00000002;
                0.98583020,-0.03338850, 0.00014064, 0.00000190, 0.00000001;
                2.00337487,-0.01141548, 0.00041590, 0.00000373, 0.00000001;
                3.00636321,-0.00683394, 0.00029281, 0.00000161, 0;
                4.00518423,-0.00501795, 0.00021735, 0.00000088, 0;
                0.47972120,-0.01456692,-0.00018739,-0.00000218,-0.00000001];

A  =polyCoeffTable(1,1)+sidelobedB*(polyCoeffTable(1,2)+sidelobedB*(polyCoeffTable(1,3)+sidelobedB*(polyCoeffTable(1,4)+sidelobedB*polyCoeffTable(1,5))));
xi1=polyCoeffTable(2,1)+sidelobedB*(polyCoeffTable(2,2)+sidelobedB*(polyCoeffTable(2,3)+sidelobedB*(polyCoeffTable(2,4)+sidelobedB*polyCoeffTable(2,5))));
xi2=polyCoeffTable(3,1)+sidelobedB*(polyCoeffTable(3,2)+sidelobedB*(polyCoeffTable(3,3)+sidelobedB*(polyCoeffTable(3,4)+sidelobedB*polyCoeffTable(3,5))));
xi3=polyCoeffTable(4,1)+sidelobedB*(polyCoeffTable(4,2)+sidelobedB*(polyCoeffTable(4,3)+sidelobedB*(polyCoeffTable(4,4)+sidelobedB*polyCoeffTable(4,5))));
xi4=polyCoeffTable(5,1)+sidelobedB*(polyCoeffTable(5,2)+sidelobedB*(polyCoeffTable(5,3)+sidelobedB*(polyCoeffTable(5,4)+sidelobedB*polyCoeffTable(5,5))));
%p0 =polyCoeffTable(6,1)+sidelobedB*(polyCoeffTable(6,2)+sidelobedB*(polyCoeffTable(6,3)+sidelobedB*(polyCoeffTable(6,4)+sidelobedB*polyCoeffTable(6,5))));

Z=zeros(N+1,1);
%Equation 15
Z(1)=0;%The Z(0) term
%Now, the moved zeros
Z(2)=xi1;
Z(3)=xi2;
Z(4)=xi3;
Z(5)=xi4;
for k=5:N
    %The location of the non-moved zeros as given by Equation 13.
    Z(k+1)=sqrt(A^2+k^2);
end
%Equation 16 in [1].
sigma=mu(N+1)/Z(N+1);

%No normalization is performed. We just use C=1.
C=1;

B=zeros(N,1);%The first entry is B_0
%Equation 24
for m=0:(N-1)
    num=prod(1-(mu(m+1)./(sigma*Z(2:N))).^2);
    denom=prod(1-(mu(m+1)./mu([1:m,(m+2):N])).^2);

    B(m+1)=-(C*1j*2*mu(m+1)^2/besselj(1,pi*mu(m+1)))*num/denom;
end

%If discretized tapering values are desired.
if(nargin>2&&~isempty(xyPoints))
    numPoints=size(xyPoints,2);

    if(nargin<4||isempty(a))
        %The maximum distance from the origin to a point is taken to be the
        %radius of the aperture.
        a2=max(sum(xyPoints.^2,1));
        a=sqrt(a2);
    end

    g=zeros(numPoints,1);
    for curPoint=1:numPoints
        rho2=sum(xyPoints(:,curPoint).^2,1);

        if(rho2<=a2)
            rho=sqrt(rho2);
            %The normalized radius at this point.
            p=pi*rho/a;

            x=xyPoints(1,curPoint);
            cosVal=x/rho;

            %Equation 7 in [1].
            g(curPoint)=cosVal*sum(B.*besselj(1,mu(1:N)*p));
        end
    end
else%If no tapering values are requested, return an empty matrix.
    g=[];
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.


text:
%
% state = [x, y, yaw, delta]
% input = [v_des, delta_des]
% ref = [x_ref, y_ref, yaw_ref, v_ref]
%
%
%

clear variables;
close all;

set(0, 'defaultAxesFontSize', 12);
set(0, 'defaultTextFontSize', 20);
set(0, 'DefaultAxesLineWidth', 1.0, 'DefaultLineLineWidth', 1.0);

addpath ../path_design

control_mode_option = ["pure_pursuit", "pid", "mpc", "mpc_no_constraints"];
control_mode = control_mode_option(3);

save_video = 0; %1:save, 0:no

%% preliminaries
rad2deg = 180 / pi;
deg2rad = pi / 180;
kmh2ms = 1000 / 3600;

simulation_time = 35;
simulation_rk4_time_step = 0.002; % simulation time step

vel_ref = 30 * kmh2ms;

% for dynamics model
param.tau = 0.27; % steering dynamics: 1d-approximated time constant
param.wheelbase = 2.69;
param.steer_lim = 30 * deg2rad;
param.vel_max = 10;
param.vel_min = -5;

param.input_delay = 0.24; % [s]
param.control_dt = 0.03; % [s]
param.measurement_noise_stddev = [0.1, 0.1, 1.0*deg2rad, 0.5*deg2rad]; % measurement noise
% param.measurement_noise_stddev = [0,0,0,0]; % measurement noise
param.steering_steady_state_error_deg = 1;

% for pure pursuit only
param.pure_pursuit_lookahead = 8.0; % [m]

% for mpc only
param.mpc_dt = 0.1;
param.mpc_n = 30;
param.mpc_constraint_steering_deg = 30;
param.mpc_constraint_steer_rate_deg = 280;
param.mpc_model_dim = 3;
param.mpc_Q = diag([1,2]);
param.mpc_R = 0.5;
param.mpc_delay_comp_step = round(param.input_delay / param.control_dt);
% param.mpc_delay_comp_step = 0.0;

% use the input ahead of the delay time
param.mpc_sensor_delay = param.input_delay;

% for mpc2
param.mpc2_dt = 0.2;
param.mpc2_n = 10;
param.mpc2_steering_lim_deg = 40;
param.mpc2_model_dim = 4;
param.mpc2_Q = diag([1,1,0]);
param.mpc2_R = 0.05;

%% simulation parameters

% initial position (x, y, yaw, delta)
x0 = [0, 0.5, 0, 0];

ts = 0;
dt = simulation_rk4_time_step;
tf = simulation_time;
t = ts:dt:tf;

%% reference trajectory design

path_design; % using spline
load path; % x, y, yaw

ref = zeros(length(path), 6);
IDX_X = 1;
IDX_Y = 2;
IDX_XY = 1:2;
IDX_XYYAW = 1:3;
IDX_YAW = 3;
IDX_VEL = 4;
IDX_CURVATURE = 5;
IDX_TIME = 6;

IDX_STEER = 4;

path_size_scale = 15;
path(:,IDX_XY) = path(:,IDX_XY) * path_size_scale;
ref(:,IDX_XYYAW) = path(:,IDX_XYYAW);

ref(:,IDX_VEL) = ones(length(path),1)*vel_ref;

% insert curvature into path
for i = 2:length(ref)-1
    p1_ = ref(i-1,IDX_XY);
    p2_ = ref(i, IDX_XY);
    p3_ = ref(i+1, IDX_XY);
    A_ = ((p2_(1)-p1_(1))*(p3_(2)-p1_(2)) - (p2_(2)-p1_(2))*(p3_(1)-p1_(1))) / 2;
    ref(i, IDX_CURVATURE) = 4 * A_ / (norm(p1_-p2_) * norm(p2_-p3_) * norm(p3_-p1_));
end

% insert relative time into path
for i = 2:length(ref)
    v_ = ref(i,IDX_VEL);
    d_ = norm(ref(i,IDX_XY)-ref(i-1,IDX_XY));
    dt_ = d_ / v_;
    ref(i, IDX_TIME) = ref(i-1, IDX_TIME) + dt_;
end

%% simulation

if control_mode == "pure_pursuit"
    [X, U, debug] = simulate_rk4(@kinematics_model, @pure_pursuit, x0, ref, ts, dt, tf, param);
    lat_error_vec = debug(:,end);
elseif control_mode == "pid"
   [X, U, debug] = simulate_rk4(@kinematics_model, @pid_controller, x0, ref, ts, dt, tf, param);
   lat_error_vec = debug(:,end);
elseif control_mode == "mpc"
    param.mpc_solve_without_constraint = false;
    [X, U, debug] = simulate_rk4(@kinematics_model, @model_predictive_controller, x0, ref, ts, dt, tf, param);
    lat_error_vec = debug(:,end);
elseif control_mode == "mpc_no_constraints"
    param.mpc_solve_without_constraint = true;
    [X, U, debug] = simulate_rk4(@kinematics_model, @model_predictive_controller, x0, ref, ts, dt, tf, param);
    lat_error_vec = debug(:,end);
elseif control_mode == "mpc2"
    [X, U, debug] = simulate_rk4(@kinematics_model, @model_predictive_controller2, x0, ref, ts, dt, tf, param);
    lat_error_vec = debug(:,end);
end
fprintf("lattitude error: mean square = %f, max = %f", norm(lat_error_vec)/simulation_time, max(lat_error_vec));

%% movie plot

sp_num = 18;
subpl1 = 'subplot(sp_num,sp_num, sp_num+1:sp_num*12);';
subpl2 = 'subplot(sp_num,sp_num, sp_num*13+1:sp_num*15);';
subpl3 = 'subplot(sp_num,sp_num, sp_num*16+1:sp_num*18);';

% tire2steer = 12.5;

fig_trajectory_result = figure(1);
set(fig_trajectory_result, 'Position', [716 735 1026 1146]);
eval(subpl1);
plot(ref(:,1), ref(:,2),'k-.'); hold on; grid on;
% plot(X(:,1), X(:,2));
% legend('ref','tracked');
xlabel('x [m]'); ylabel('y [m]');
% img_orig = imread('handle.jpg');
% img = imrotate(img_orig, 10);
% handle_plt = image([150, 170],[110, 90], img);
% handle_plt_point = [160, 100, 0];

eval(subpl2);
plot(t, lat_error_vec, 'b'); grid on; hold on;
xlabel('t [s]'); ylabel('latitude error [m]');
ulim = ceil(2*max(lat_error_vec))/2;
dlim = floor(2*min(lat_error_vec))/2;
ylim([dlim, ulim]);

eval(subpl3);
p1 = plot(t, X(:,IDX_STEER)*rad2deg, 'b'); grid on; hold on;
p2 = plot(t, U(:,2)*rad2deg, 'Color', [0.7 0. 1]); hold on;
legend([p1,p2], {'measured','command'})
xlabel('t [s]'); ylabel('steering angle [deg]');
ulim = round(2*max(X(:,IDX_STEER)*rad2deg))/2;
dlim = round(2*min(X(:,IDX_STEER)*rad2deg))/2;
ylim([dlim, ulim]);

z_axis = [0 0 1];
setpoint = []; rear_tire = []; front_tire = []; body = []; tracked = [];
setpoint_ideal = []; error_point = []; steer_point = []; time_bar_laterror = []; time_bar_steer = [];
L = param.wheelbase;
rear_length = 1;
front_length = 1;
side_width = 0.9;
fig_draw_i = 1:round(1/dt/20):length(t);

% for movie
clear frame_vec;
frame_vec(length(fig_draw_i)) = struct('cdata', [], 'colormap',[]);

j = 1;
fig_trajectory_result; hold on;
for i = fig_draw_i
    eval(subpl1);
    disp(t(i))
    rear_x = X(i,1);
    rear_y = X(i,2);
    yaw = X(i,3);
    delta = X(i,IDX_STEER);
    front_x = rear_x + L;
    front_y = rear_y;
    delete([setpoint, rear_tire, front_tire, body, tracked, setpoint_ideal, error_point, steer_point, time_bar_laterror, time_bar_steer]);
    tracked = plot(X(1:i,1), X(1:i,2),'r');
    title_draw = "t = "+num2str(t(i),'%5.1f') + "[s], steer = " + num2str(delta*rad2deg,'%+3.1f') + "[deg], v = " + ...
        num2str(vel_ref*3600/1000,'%3.1f') + "[km/h], lat error = "+num2str(lat_error_vec(i),'%+2.2f') + "[m]";
    title_draw = [title_draw; "Simulation: solver = rk4, sensor-delay = "  + num2str(param.input_delay*1000, '%d') + "[ms], control freq=" + ...
        num2str(1/param.control_dt, '%d') + "[hz]"];
    title_draw = [title_draw; "noise-sigma = " + num2str(param.measurement_noise_stddev(1),'%2.2f')+"(pos), "+ ...
        num2str(param.measurement_noise_stddev(3),'%2.2f')+"(yaw), "+num2str(param.measurement_noise_stddev(4),'%2.2f')+"(steer)"];
    if control_mode == "mpc" || control_mode == "mpc_no_constraints"
        title_draw = [title_draw; "MPC: dt = " + num2str(param.mpc_dt, '%3.3f') + "[s], horizon step = " + num2str(param.mpc_n, '%d')];
%         pred_states = debug(i, param.mpc_n+1:param.mpc_n*(4+1));
%         pred_states = reshape(pred_states, param.mpc_n, length(pred_states)/param.mpc_n);
%         setpoint = plot(pred_states(:,1), pred_states(:,2), 'bo'); % include linealize error
        pred_error = debug(i, param.mpc_n*(4+1)+1:param.mpc_n*(2+4+1));
        pred_error = reshape(pred_error, param.mpc_n, length(pred_error)/param.mpc_n);
        setpoint_ideal = plot(pred_error(:,1), pred_error(:,2), 'mx'); % without linealize error
    elseif control_mode == "mpc2"
        title_draw = [title_draw; "MPC2: dt = " + num2str(param.mpc_dt, '%3.3f') + "[s], horizon step = " + num2str(param.mpc_n, '%d')];
        pred_states = debug(i, param.mpc_n+1:param.mpc_n*(4+1));
        pred_states = transpose(reshape(pred_states, length(pred_states)/param.mpc_n, param.mpc_n));
        setpoint = plot(pred_states(:,1), pred_states(:,2), 'bo'); % include linealize error
        pred_error = debug(i, param.mpc_n*(4+1)+1:param.mpc_n*(4+4+1));
        pred_error = transpose(reshape(pred_error, length(pred_error)/param.mpc_n, param.mpc_n));
        setpoint_ideal = plot(pred_error(:,1), pred_error(:,2), 'mx'); % without linealize error
    elseif control_mode == "pure_pursuit"
        title_draw = [title_draw; "pure-pursuit: lookahead dist="+num2str(param.pure_pursuit_lookahead, '%1.1f')+"[m]"];
        sp = debug(i,:);
        setpoint = plot(sp(1), sp(2), 'ro');
    elseif control_mode == "pid"
        title_draw = [title_draw; "PID: kp = " + num2str(0.3, '%3.3f') + ", ki = " + num2str(0, '%3.3f') + ", kd = " + num2str(1.5, '%3.3f')];
        sp = debug(i,:);
        setpoint = plot(sp(1), sp(2), 'ro');
    end
    rear_tire = plot([rear_x-0.3, rear_x+0.3],[rear_y, rear_y], 'k', 'LineWidth', 2.0);
    front_tire = plot([front_x-0.3, front_x+0.3],[front_y, front_y], 'k', 'LineWidth', 2.0);
    body = plot([rear_x-rear_length, front_x+front_length, front_x+front_length, rear_x-rear_length, rear_x-rear_length], ...
        [rear_y-side_width, front_y-side_width, front_y+side_width, rear_y+side_width, rear_y-side_width],'k');
    rear_origin = [rear_x, rear_y, 0];
    front_origin = [rear_x + L*cos(yaw), rear_y + L*sin(yaw), 0];
    rotate(body, z_axis, yaw * rad2deg, rear_origin);
    rotate(rear_tire, z_axis, yaw * rad2deg, rear_origin);
    rotate(front_tire, z_axis, yaw * rad2deg, rear_origin);
    rotate(front_tire, z_axis, delta * rad2deg, front_origin);
%     rotate(handle_plt, [0,0,1], delta*tire2steer*rad2deg, handle_plt_point)
    title(title_draw);
    xlim([0 120]);

    % lat error
    eval(subpl2);
    error_point = plot(t(i), lat_error_vec(i), 'ko');
    time_bar_laterror = plot([t(i), t(i)], [100, -100], 'k');

    % steering
    eval(subpl3);
    steer_point = plot(t(i), X(i, IDX_STEER)*rad2deg, 'ko');
    time_bar_steer = plot([t(i), t(i)], [100, -100], 'k');
    legend([p1,p2], {'measured','command'})
    ylim([-40 40]);

    drawnow;
    frame_vec(j) = getframe(fig_trajectory_result);

    j = j + 1;
end

% for video save
if (save_video == 1)
    cd ./movie
    frame_vec(1) = [];
    vidObj = VideoWriter('result.avi');
    vidObj.FrameRate = 25;
    open(vidObj);
    writeVideo(vidObj, frame_vec);
    close(vidObj);
    cd ../
end

% plot_pid;

author: TakaHoribe


text:
classdef (InferiorClasses = {?rotation,?quaternion}) SO3Grid < orientation
% represent orientations in a gridded structure to allow quick access
%
% Syntax
%   S3G = regularSO3Grid(CS)
%   S3G = regularSO3Grid(CS,SS,'resolution',2.5*degree)     % specify the resolution
%   S3G = regularSO3Grid(CS,SS,'resolution',5*degree,'ZYZ') % use ZYZ convention
%   S3G = regularSO3Grid(CS,SS,'phi2','sections',10) % 10 phi2 sections
%   S3G = regularSO3Grid
%   S3G = equispacedSO3Grid(CS,SS,'points',n)
%   S3G = equispacedSO3Grid(CS,'resolution',res)
%
%   % fill only a ball with radius of 20 degree
%   S3G = equispacedSO3Grid(CS,'maxAngle',20*degree)
%
% Input
%  CS  - @crystalSymmetry
%  SS  - @specimenSymmetry
%   n  - approximate number of points
%  res - resolution in radiant
%
% Output
%  S3G - @SO3Grid
%
% Options
%  maxAngle - radius of the ball to be filles
%  center - center of the ball
%  sections - number of sections
%
% Flags
%  phi1 | Phi | phi2 - sections along which variable
%

  properties
    alphabeta = [];
    gamma    = [];
    resolution = 2*pi;
    center  = [];
  end

  methods
    function S3G = SO3Grid(ori,alphabeta,gamma,varargin)

      % call superclass method
      S3G = S3G@orientation(ori);

      S3G.alphabeta = alphabeta;
      S3G.gamma = gamma;
      S3G.center = get_option(varargin,'center',quaternion.id);
      S3G.resolution = get_option(varargin,'resolution',2*pi);
    end
  end
end