text:
%     Steffen Urban email: steffen.urban@kit.edu
%     Copyright (C) 2014  Steffen Urban
%
%     This program is free software; you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation; either version 2 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License along
%     with this program; if not, write to the Free Software Foundation, Inc.,
%     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

clear all
close all

% !!! change those if you changed Step1_perform_test_calibrations.m !!!
nr_cams    = 7;
nr_tests   = 5;

%% load scaramuzza data
for i=1:nr_tests
    path = sprintf('CalibData%i.mat', i);
    calibMeth{i} = load(path);
end

for cam=1:nr_cams
    for meth = 1:nr_tests
        rms(cam,meth) = calibMeth{meth}.calib_data{cam}.rms;
    end
end

figure
bar(rms)
title('root mean square error')
xlabel('data set')
ylabel('root mean square error [pixel]')
legend('ocam standard','ocam subpixel','urban', 'urban subpixel' ,'urban subpixel robust')

for cam=1:nr_cams
    for meth = 1:nr_tests
        runtime(cam,meth) = calibMeth{meth}.calib_data{cam}.runtime;
    end
end

figure
bar(runtime)
title('runtime')
xlabel('data set')
ylabel('time [s]')
legend('ocam standard','ocam subpixel','urban', 'urban subpixel' ,'urban subpixel robust')

badIdx = [];
for cam=1:nr_cams
    for meth = 1:2
        stdEOangle = [];
        stdEOpos = [];
        for imgs = calibMeth{meth}.calib_data{cam}.ima_proc
            try % if an image was omitted during corner extraction
            if (~isempty(calibMeth{meth}.calib_data{cam}.statEO{imgs}.stdEO))
                stdEOangle = [stdEOangle
                              1e3*[calibMeth{meth}.calib_data{cam}.statEO{imgs}.stdEO(1)
                                   calibMeth{meth}.calib_data{cam}.statEO{imgs}.stdEO(2)
                                   calibMeth{meth}.calib_data{cam}.statEO{imgs}.stdEO(3)]];
                stdEOpos   = [stdEOpos
                              [calibMeth{meth}.calib_data{cam}.statEO{imgs}.stdEO(4)
                              calibMeth{meth}.calib_data{cam}.statEO{imgs}.stdEO(5)
                              calibMeth{meth}.calib_data{cam}.statEO{imgs}.stdEO(6)]];
               badIdx = [badIdx, imgs];
            end
            catch

            end
        end
        stdEOangleA(cam,meth) = mean(stdEOangle);
        stdEOposA(cam,meth) = mean(stdEOpos);
        badIdx = [];
    end
end

lauf = 0;
for cam=1:nr_cams
    for meth = 3:nr_tests
        stdEOangle = [];
        stdEOpos = [];
        for imgs = 1:length(calibMeth{meth}.calib_data{cam}.ima_proc)
            stdEOangle = [stdEOangle
                          1e3*[calibMeth{meth}.calib_data{cam}.statEO.stdEO(1+lauf)
                               calibMeth{meth}.calib_data{cam}.statEO.stdEO(2+lauf)
                               calibMeth{meth}.calib_data{cam}.statEO.stdEO(3+lauf)]];
            stdEOpos   = [stdEOpos
                          [calibMeth{meth}.calib_data{cam}.statEO.stdEO(4+lauf)
                          calibMeth{meth}.calib_data{cam}.statEO.stdEO(5+lauf)
                          calibMeth{meth}.calib_data{cam}.statEO.stdEO(6+lauf)]];
            lauf = lauf+6;
        end
        lauf = 0;
        stdEOangleA(cam,meth) = median(stdEOangle);
        stdEOposA(cam,meth) = median(stdEOpos);
    end
end

figure
bar(stdEOangleA)
title('mean standard deviation of camera orientations')
xlabel('data set')
ylabel('mean standard deviation[mrad]')
legend('ocam standard','ocam subpixel','urban', 'urban subpixel' ,'urban subpixel robust')

figure
bar(stdEOposA)
title('mean standard deviation of camera positions')
xlabel('data set')
ylabel('standard deviation [mm]')
legend('ocam standard','ocam subpixel','urban', 'urban subpixel' ,'urban subpixel robust')

%% IO, c and a0
lauf = 0;
for cam=1:nr_cams
    for meth = 1:nr_tests
        std_cde(cam,meth) = calibMeth{meth}.calib_data{cam}.statIO.stdIO(3);
        std_a0(cam,meth) = calibMeth{meth}.calib_data{cam}.statIO.stdIO(6);
    end
end

author: urbste

repo: ImprovedOcamCalib


text:
function [value, isterminal, direction, causes] = getSoITransitionOdeEvents(ut, rVect, vVect, bodyInfo, celBodyData)
        value = [];
        isterminal = [];
        direction = [];
%         causes = AbstractIntegrationTerminationCause.empty(0,1);

        %Max Radius (SoI Radius) Constraint (Leave SOI Upwards)
        parentBodyInfo = bodyInfo.getParBodyInfo(celBodyData);
%         rSOI = getSOIRadius(bodyInfo, parentBodyInfo);
        rSOI = bodyInfo.getCachedSoIRadius();
        radius = norm(rVect);

%         if(isempty(parentBodyInfo))
%             parentBodyInfo = KSPTOT_BodyInfo.empty(0,1);
%         end

        value(end+1) = rSOI - radius;
        isterminal(end+1) = 1;
        direction(end+1) = -1;
        causes(1) = SoITransitionUpIntTermCause(bodyInfo, parentBodyInfo, celBodyData);

        %Leave SoI Downwards
        [sma, ecc, ~, ~, ~, ~] = getKeplerFromState(rVect, vVect, bodyInfo.gm, false);
        [rApSC, rPeSC] = computeApogeePerigee(sma, ecc);

        if(ecc >= 1)
            rApSC = Inf;
        end

        children = bodyInfo.getChildrenBodyInfo(celBodyData);
        if(~isempty(children))
%             soiDownCauses(length(children)) = SoITransitionDownIntTermCause(bodyInfo, children(end), celBodyData);
            for(i=1:length(children)) %#ok<*NO4LP>
                childBodyInfo = children(i);
                rSOI = childBodyInfo.getCachedSoIRadius();
                [rApCB, rPeCB] = computeApogeePerigee(childBodyInfo.sma, childBodyInfo.ecc);

                if((rApSC < (rPeCB - rSOI)) || ...
                    rPeSC > (rApCB + rSOI))
                    val = realmax;

                else
                    dVect = getAbsPositBetweenSpacecraftAndBody(ut, rVect, bodyInfo, childBodyInfo, celBodyData);
                    distToChild = norm(dVect);

                    val = distToChild - rSOI;
                end

                value(end+1) = val; %#ok<AGROW>
                direction(end+1) = -1; %#ok<AGROW>
                isterminal(end+1) = 1; %#ok<AGROW>
                soiDownCauses(i) = SoITransitionDownIntTermCause(bodyInfo, childBodyInfo, celBodyData); %#ok<AGROW>
            end
            causes = [causes, soiDownCauses];
        end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function exact = p04_exact ( )

%*****************************************************************************80
%
%% P04_EXACT returns the estimated integral for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the estimated value of the integral.
%
  exact = -0.0046848541335080643181;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function lvec_print ( n, a, title )

%*****************************************************************************80
%
%% LVEC_PRINT prints a logical vector.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%
%    Input, logical A(N), the vector to be printed.
%
%    Input, string TITLE, a title to be printed first.
%    TITLE may be blank.
%
  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  fprintf ( 1, '\n' );
  for i = 1 : n
    value = ( a(i) ~= 0 );
    fprintf ( 1, '%6d  %1d\n', i, value );
  end

  return
end


text:
%DRAWCHARUCOBOARD  Draw a ChArUco board
%
%     img = cv.drawCharucoBoard(board, outSize)
%     img = cv.drawCharucoBoard(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __board__ ChArUco board that will be drawn.
% * __outSize__ size of the output image in pixels `[w,h]`.
%
% ## Output
% * __img__ output image with the board. The size of this image will be
%   `outSize` and the board will be on the center, keeping the board
%   proportions.
%
% ## Options
% * __MarginSize__ minimum margins (in pixels) of the board in the output
%   image. default 0
% * __BorderBits__ width of the marker borders. default 1
%
% This function return the image of the ChArUco board, ready to be printed.
%
% See also: cv.estimatePoseCharucoBoard, cv.drawPlanarBoard, cv.drawMarkerAruco
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%% Naive Bayes with independent Bernoulli
close all; clear;
d = 10;
k = 2;
n = 2000;
[X,t,mu] = mixBernRnd(d,k,n);
m = floor(n/2);
X1 = X(:,1:m);
X2 = X(:,(m+1):end);
t1 = t(1:m);
t2 = t((m+1):end);
model = nbBern(X1,t1);
y2 = nbBernPred(model,X2);
err = sum(t2~=y2)/numel(t2);

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function order = lyness_order ( rule )

%*****************************************************************************80
%
%% LYNESS_ORDER returns the order of a Lyness quadrature rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    James Lyness, Dennis Jespersen,
%    Moderate Degree Symmetric Quadrature Rules for the Triangle,
%    Journal of the Institute of Mathematics and its Applications,
%    Volume 15, Number 1, February 1975, pages 19-32.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%
%    Output, integer ORDER, the order of the rule.
%
  if ( rule == 0 )
    order = 1;
  elseif ( rule == 1 )
    order = 3;
  elseif ( rule == 2 )
    order = 4;
  elseif ( rule == 3 )
    order = 4;
  elseif ( rule == 4 )
    order = 7;
  elseif ( rule == 5 )
    order = 6;
  elseif ( rule == 6 )
    order = 10;
  elseif ( rule == 7 )
    order = 9;
  elseif ( rule == 8 )
    order = 7;
  elseif ( rule == 9 )
    order = 10;
  elseif ( rule == 10 )
    order = 12;
  elseif ( rule == 11 )
    order = 16;
  elseif ( rule == 12 )
    order = 13;
  elseif ( rule == 13 )
    order = 13;
  elseif ( rule == 14 )
    order = 16;
  elseif ( rule == 15 )
    order = 16;
  elseif ( rule == 16 )
    order = 21;
  elseif ( rule == 17 )
    order = 16;
  elseif ( rule == 18 )
    order = 19;
  elseif ( rule == 19 )
    order = 22;
  elseif ( rule == 20 )
    order = 27;
  elseif ( rule == 21 )
    order = 28;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LYNESS_ORDER - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized rule index.\n' );
    error ( 'LYNESS - Fatal error!' );
  end

  return
end


text:
function option_num = p16_option_num ( )

%*****************************************************************************80
%
%% P16_OPTION_NUM returns the number of options for problem 16.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer OPTION_NUM, the number of options.
%
  option_num = 1;

  return
end


text:
function sphere_stereograph_test03 ( )

%*****************************************************************************80
%
%% SPHERE_STEREOGRAPH_TEST03 checks that the two functions are inverses.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_STEREOGRAPH_TEST03\n' );
  fprintf ( 1, '  SPHERE_STEREOGRAPH2 maps from sphere to plane.\n' );
  fprintf ( 1, '  SPHERE_STEREOGRAPH2_INVERSE is the inverse map.\n' );
  fprintf ( 1, '  Check that these two functions are inverses.\n' );

  m = 3;
  n = 100;
  seed = 123456789;

  [ focus, seed ] = r8vec_uniform_01 ( m, seed );
  [ center, seed ] = r8vec_uniform_01 ( m, seed );
  r = norm ( focus - center );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Using radius = %f\n', r );
  fprintf ( 1, '  center = (%f,%f,%f)\n', center );
  fprintf ( 1, '  focus =  (%f,%f,%f)\n', focus );
%
%  Check #1.
%
  p1 = uniform_on_sphere01_map ( m, n, seed );
  p1 = repmat ( center, 1, n ) + r * p1;

  q = sphere_stereograph2 ( p1, focus, center );

  p2 = sphere_stereograph2_inverse ( q, focus, center );

  dif = sqrt ( sum ( ( p1 - p2 ).^2, 1 ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Map points from sphere to plane to sphere.\n' );
  fprintf ( 1, '  Maximum difference for %d points was %f\n', n, max ( dif ) );
%
%  Check #2.
%  We have to work hard to get random points on the plane, since
%  all we know to begin with is the point of tangency and the normal.
%
  tang = 2 * center - focus;
  normal = center - focus;
  [ pr, pq ] = plane_normal_basis_3d ( tang, normal );
  q1 = zeros(m,n);
  alpha = rand(1,n);
  beta = rand(1,n);
  tt = repmat ( tang, 1, n );

  q1(1:m,1:n) = tt + pr(1:m,1) * alpha(1,1:n) + pq(1:m,1) * beta(1,1:n);
  p = sphere_stereograph2_inverse ( q1, focus, center );
  q2 = sphere_stereograph2 ( p, focus, center );

  dif = sqrt ( sum ( ( q1 - q2 ).^2, 1 ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Map points from plane to sphere to plane.\n' );
  fprintf ( 1, '  Maximum difference for %d points was %f\n', n, max ( dif ) );

  return
end


text:
function parameters = convert_ma_roots(parameters, q)
% Checks an MA parameterization for invertability and inverts it if possible
%
%
%
% COMMENTS:
% It is not always possible to invert irregular MAs since the indices of the non-zero coefficients
% in the invertable MA may differ.  If this is the case PARAMETERS will be equal to the input.

MAparameters = zeros(1,max(q));
MAparameters(q) = parameters;
lambda = roots([1 MAparameters]);
lambda(abs(lambda)>1) = 1./lambda(abs(lambda)>1);
parameters = poly(lambda);
parameters(abs(parameters)<1e-10) = 0;

if length(q)<max(q)
    % Check for index change
    newq = find(abs(parameters)>1e-10);
    newq = newq(2:length(newq))-1;
    if length(q)~=length(newq) || any(q~=newq)
        warning('MFEToolbox:Invertability',['The irregular MA cannot be inverted without changing the lag indices. \nThe require lag indices are:' num2str(newq)])
        parameters = MAparameters(q);
    else
        parameters = parameters(q+1);
    end
else
    parameters = parameters(q+1);
end

if size(parameters,2)>size(parameters,1)
    parameters = parameters';
end

author: bashtage

repo: mfe-toolbox


text:
function value = legendre_determinant ( n )

%*****************************************************************************80
%
%% LEGENDRE_DETERMINANT returns the determinant of the LEGENDRE matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real VALUE, the determinant.
%
  value = 1.0;
  t = 1.0;

  for i = 3 : n
    t = t * ( 2 * i - 3 ) / (  i - 1 );
    value = value * t;
  end

  return
end


text:
function ti = treeinfo(X, y, treeopts, tstopts)
%TREEINFO Collects decision tree information based on tree walking.
%
%   ti = treeinfo(X, y, treeopts, tstopts)
%
%   It creates CART according to data set [X, y] and
%   optional parameters in treeopts, performs
%   tree pruning (if required) according to optional
%   parameters in tstopts, and implements tree walking
%   to collect all necessary information for genfis4.
%   treeopts and tstopts should be either empty matrix or
%   cell array containing one or more optional
%   parameter name/value pairs.
%
%   The result is a structure ti containing the following fields:
%   - structure array "node" with fields "variable" (integer-valued scalar),
%       "cutpoint" (float scalar), and "sample" (index-valued vector);
%   - structure array "leaf" with field "sample" (index-valued vector);
%   - structure array "branch" with fields "nodes" (integer-valued vector),
%       and "ineqs" (vector of "1" or "-1" values).
%
%   It requires Statistics Toolbox realizing CART algorithm.
%
%   Example:
%   load fisheriris;
%   treeopts = {'minparent', 5, 'prune', 'off'};
%   tstopts = {'crossvalidation', meas, species};
%   ti = treeinfo(meas, species, treeopts, tstopts);

%   Per Konstantin A. Sidelnikov, 2009.

% Check inputs
error(nargchk(2, 4, nargin));

if (nargin < 4) || isempty(tstopts)
    tstopts = {};
end
if (nargin < 3) || isempty(treeopts)
    treeopts = {};
end

% Create a decision tree using CART algorithm
t = classregtree(X, y, treeopts{:});
% If necessary, perfom optimal tree pruning
if ~isempty(tstopts)
    [~, ~, ~, bestlevel] = t.test(tstopts{:});
    t = t.prune('level', bestlevel);
end
% Get indices of nodes and leaves (terminal nodes)
indn = find(t.isbranch());
indl = find(~t.isbranch());
% Get overall numbers of nodes and leaves
numn = length(indn);
numl = length(indl);
% Allocate memory for structure array,
% which stores information about each branch nodes:
% cut variable number, cut point, and observations from
% the original data that satisfy the conditions for the node
node(numn).variable = [];
node(numn).cutpoint = [];
node(numn).sample = [];
% Allocate memory for structure array,
% which stores information about each leaves:
% observations from the original data that
% satisfy the conditions for the leaf
leaf(numl).sample = [];
% Allocate memory for structure array,
% which stores information about each branches:
% group of nodes that form a path to the leaf, and
% types of inequalities that are satisfied to reach the leaf
branch(numl).nodes = [];
branch(numl).ineqs = [];
% Implement tree walking
treewalk(1, [], [], 1 : length(y));
% Collect tree information
ti = struct('node', node, 'leaf', leaf, 'branch', branch);

    function treewalk(n, nodes, ineqs, sample)
        %   This nested function implements tree walking
        %   based on recursive algorithm in reference and
        %   collects tree information using data sharing.

        %   Reference: Cormen T.H., Leiserson C.E., Rivest R.L.,
        %   Introduction to Algorithms (1st ed.),
        %   MIT Press and McGraw-Hill, 1990.

        % Get numbers of the child nodes
        kids = t.children(n);
        if all(kids) % branch node
            % Get cut variable and cut point
            var = t.var(n);
            cut = t.cutpoint(n);
            % Save cut variable, cut point, and observations
            k = find(n == indn, 1, 'first');
            node(k).variable = var;
            node(k).cutpoint = cut;
            node(k).sample = sample;
            % Form new subsample and continue walking to the left child
            s = sample(X(sample, var) < cut);
            treewalk(kids(1), [nodes, k], [ineqs, -1], s);
            % Form new subsample and continue walking to the right child
            s = sample(X(sample, var) >= cut);
            treewalk(kids(2), [nodes, k], [ineqs, 1], s);
        else % leaf
            % Save leaf and branch information
            k = find(n == indl, 1, 'first');
            leaf(k).sample = sample;
            branch(k).nodes = nodes;
            branch(k).ineqs = ineqs;
        end
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function [rsamobjects, ah]=plotrsam(sta, chan, snum, enum, DATAPATH)
% [rsamobjects, ah]=plotrsam_wrapper(sta, chan, snum, enum, DATAPATH)
%
%   Inputs:
%       sta - station code
%       chan - channel code
%       snum - start datenum
%       enum - end datenum
%       DATAPATH - path to data, including pattern
%
%   Outputs:
%       rsamobjects - vector of rsam objects
%       ah - vector of axes handles
%
%   Examples:
%       1. Data from the digital seismic network, Montserrat
%           DP = fullfile('/raid','data','antelope','mvo','SSSS_CCC_YYYY.DAT');
%           [rsamobjects, ah] = rsam.plotrsam('MBWH','SHZ',datenum(2001,2,24), datenum(2001,3,3), DP);
%       2. Data from the analog seismic network, Montserrat
%           DP = fullfile(DROPBOX, 'DOME', 'SEISMICDATA', 'RSAM_1', 'SSSSYYYY.DAT');
%           [rsamobjects, ah] = rsam.plotrsam('MWHZ','',datenum(1996,7,1), datenum(1996,8,13), DP);
%   Could use the following logic in a wrapper to decide DP:
%       strfind(sta{i},'MB') & ~strcmp(sta{i},'MBET')

    % validate
    if nargin ~= 5
        help rsam>plotrsam()
        return
    end

    % initialise
    if ~iscell(sta)
        sta={sta};
    end
    if ~iscell(chan)
        chan={chan};
    end
    numsta = length(sta);
    numrsams = 0;
    rsamobjects = [];
    ah = [];

    % load data
    for i=1:numsta
        s = rsam('file', DATAPATH, 'snum', snum, 'enum', enum, 'sta', sta{i},'chan',chan{i});
        if ~isempty(s.data)
            numrsams = numrsams + 1;
            rsamobjects = [rsamobjects resample(s.despike(100))];
            %ah(i)=subplot(numsta,1,i),plot(resample(s.despike(10)));
        end
    end

    % plot data
    if numrsams > 0
        figure
        for i=1:numrsams
            ah(i) = subplot(numrsams, 1, i), plot(rsamobjects(i))
        end
        linkaxes(ah, 'x')
        %datetick('x','keeplimits')
    end
end

author: geoscience-community-codes

repo: GISMO


text:
% Driver script for solving the 3D advection equations
Globals3D;

% Order of polymomials used for approximation
N = 8;

% Generate simple mesh
filename = 'cubeK6.neu'
[Nv, VX, VY, VZ, K, EToV] = MeshReaderGambit3D(filename);

% Initialize solver and construct grid and metric
StartUp3D;

% set initial conditions
u = exp(-1*(x.^2 + y.^2 + z.^2));

% Solve Problem
FinalTime = 0.8;
[u] = Advec3D(u,FinalTime);

repo: Aero-matlab


text:
% MANIPULATOR TRAJECTORY GENERATION
% Generates Cartesian trajectories with independent
% linearly interpolated rotation trajectories.
%
% Copyright 2019 The MathWorks, Inc.

%% Setup
clear, clc, close all

% Define waypoint information
createWaypointData;

% Define IK
ik = inverseKinematics('RigidBodyTree',gen3);
ikWeights = [1 1 1 1 1 1];
ikInitGuess = gen3.homeConfiguration;

% Set up plot
plotMode = 2; % 0 = None, 1 = Trajectory, 2 = Coordinate Frames
show(gen3,gen3.homeConfiguration,'Frames','off','PreservePlot',false);
xlim([-1 1]), ylim([-1 1]), zlim([0 1.2])
hold on
if plotMode == 1
    hTraj = plot3(waypoints(1,1),waypoints(2,1),waypoints(3,1),'b.-');
end
plot3(waypoints(1,:),waypoints(2,:),waypoints(3,:),'ro','LineWidth',2);

%% Generate and follow trajectory
% Loop through segments one at a time
trajType = 'trap';
numWaypoints = size(waypoints,2);
for w = 1:numWaypoints-1
    % Get the initial and final rotations and times for the segment
    R0 = eul2quat(orientations(:,w)');
    Rf = eul2quat(orientations(:,w+1)');
    timeInterval = waypointTimes(w:w+1);
    trajTimes = timeInterval(1):ts:timeInterval(2);

    % Cartesian Motion only
    switch trajType
        case 'trap'
            [q,qd,qdd] = trapveltraj(waypoints(:,w:w+1),numel(trajTimes), ...
                'AccelTime',waypointAccelTimes(w), ...
                'EndTime',diff(waypointTimes(w:w+1)));

        case 'cubic'
            [q,qd,qdd] = cubicpolytraj(waypoints(:,w:w+1),waypointTimes(w:w+1),trajTimes, ...
                'VelocityBoundaryCondition',waypointVels(:,w:w+1));

        case 'quintic'
            [q,qd,qdd] = quinticpolytraj(waypoints(:,w:w+1),waypointTimes(w:w+1),trajTimes, ...
                'VelocityBoundaryCondition',waypointVels(:,w:w+1), ...
                'AccelerationBoundaryCondition',waypointAccels(:,w:w+1));

        case 'bspline'
            ctrlpoints = waypoints(:,idx:idx+1); % Can adapt this as needed
            [q,qd,qdd] = bsplinepolytraj(ctrlpoints,timeInterval,trajTimes);

        otherwise
            error('Invalid trajectory type! Use ''trap'', ''cubic'', ''quintic'', or ''bspline''');
    end

    % Find the quaternions from trajectory generation
    [R, omega, alpha] = rottraj(R0, Rf, timeInterval, trajTimes);

    % Plot trajectory
    if plotMode == 1
        set(hTraj,'xdata',q(1,:),'ydata',q(2,:),'zdata',q(3,:));
    elseif plotMode == 2
        plotTransforms(q',R','FrameSize',0.05)
    end

    % Trajectory following loop
    for idx = 1:numel(trajTimes)
        % Solve IK
        tgtPose = trvec2tform(q(:,idx)') * quat2tform(R(:,idx)');
        [config,info] = ik(eeName,tgtPose,ikWeights,ikInitGuess);
        ikInitGuess = config;

        % Show the robot
        show(gen3,config,'Frames','off','PreservePlot',false);
        title(['Trajectory at t = ' num2str(trajTimes(idx))])
        drawnow
    end

end

author: mathworks-robotics

repo: trajectory-planning-robot-manipulators


text:
function lA = shortAxis(grains,varargin)
% short axis of a grain
%
% the long axis is the direction of the smallest
% <grain2d.principalComponents.html,principal component> of a grain
%
% Syntax
%   sA = grains.shortAxis
%
% Input
%  grains - @grain2d
%
% Output
%  sA - @vector3d direction of the shortest elongation
%
% See also
% grain2d/principalComponents

omega = principalComponents(grains);

lA = vector3d.byPolar(pi/2,omega+pi/2,'antipodal');

author: mtex-toolbox

repo: mtex

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


text:
%DRAWCHESSBOARDCORNERS  Renders the detected chessboard corners
%
%     im = cv.drawChessboardCorners(im, patternSize, corners)
%     im = cv.drawChessboardCorners(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __im__ Source image. It must be an 8-bit color image.
% * __patternSize__ Number of inner corners per a chessboard row and column
%   (`patternSize = [points_per_row, points_per_column]`).
% * __corners__ Array of detected corners, the output of
%   cv.findChessboardCorners.
%
% ## Output
% * __im__ Destination image.
%
% ## Options
% * __PatternWasFound__ Parameter indicating whether the complete board was
%   found or not. The return value of cv.findChessboardCorners should be
%   passed here. default true
%
% The function draws individual chessboard corners detected either as red
% circles if the board was not found, or as colored corners connected with
% lines if the board was found.
%
% See also: cv.findChessboardCorners, cv.findCirclesGrid
%


text:
function pass = test_partitionCombine( )
% Test the partition and combine methods.

tol = 1e3*chebfunpref().cheb2Prefs.chebfun2eps;

% Check that partitioning an empty spherefun gives two empty spherefuns
f = spherefun;
[feven,fodd] = partition(f);
pass(1) = isempty(feven) & isempty(fodd);

feven = spherefun;
fodd = spherefun;
f = combine(feven,fodd);
pass(2) = isempty(f);

% Check that an even spherefun is partitioned correctly.
f = spherefun(@(x,y,z) sin(pi*x.*y));  % Strictly even/pi-periodic
[feven,fodd] = partition(f);
pass(3) = isequal(feven,f);
pass(4) = isempty(fodd);

% Check that odd spherefun is partitioned correctly.
f = spherefun(@(x,y,z) sin(pi*x.*z));  % Strictly odd/anti-periodic
[feven,fodd] = partition(f);
pass(5) = isequal(fodd,f);
pass(6) = isempty(feven);

% Check that the even and odd terms are partitioned correctly
fe = @(x,y,z) sin(pi*x.*y);  % Strictly even/pi-periodic
fo = @(x,y,z) sin(pi*x.*z);  % Strictly odd/anti-periodic
f = spherefun(@(x,y,z) fe(x,y,z) + fo(x,y,z));
[feven,fodd] = partition(f);
pass(7) = norm(spherefun(fe)-feven) < tol;
pass(8) = norm(spherefun(fo)-fodd) < tol;

% Check that combine puts an even spherefun and empty spherfun back
% together.
feven = spherefun(@(x,y,z) sin(pi*x.*y));  % Strictly even/pi-periodic
fodd = spherefun;
f = combine(feven,fodd);
pass(9) = norm(feven-f) < tol;

% Check that combine puts an even spherefun and empty spherfun back
% together.
feven = spherefun;
fodd = spherefun(@(x,y,z) sin(pi*x.*z));  % Strictly odd/anti-periodic
f = combine(feven,fodd);
pass(10) = norm(fodd-f) < tol;

% Check that combine puts non-empty odd and even spherefuns back together.
fe = @(x,y,z) sin(pi*x.*y);  % Strictly even/pi-periodic
fo = @(x,y,z) sin(pi*x.*z);  % Strictly odd/anti-periodic
fcombine = spherefun(@(x,y,z) fe(x,y,z) + fo(x,y,z));
f = combine(spherefun(fe),spherefun(fo));
pass(11) = norm(fcombine-f) < tol;

% Check that combine cannot put together spherefuns that are not strictly
% odd or even.
try
    f = spherefun(@(x,y,z) sin(pi*x.*y) + sin(pi*x.*z));
    g = combine(f,f);
    pass(12) = false;
catch ME
    pass(12) = strcmp(ME.identifier, 'CHEBFUN:SPHEREFUN:combine:parity');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [handles, hFig] = openPlotFigure(hObject, handles, Name, Title)
% openPlotFigure Opens a window for plots

hFig = figure('Position', getPlotWindowSize(), ...
              'Name', Name, ...
              'NumberTitle', 'off');
set(gca, 'Color', 'none');

hold on;
title(Title)

guidata(hObject, handles);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
% Copyright 2017 Google Inc.
%
% Licensed under the Apache License, Version 2.0 (the "License");
% you may not use this file except in compliance with the License.
% You may obtain a copy of the License at
%
%     https://www.apache.org/licenses/LICENSE-2.0
%
% Unless required by applicable law or agreed to in writing, software
% distributed under the License is distributed on an "AS IS" BASIS,
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing permissions and
% limitations under the License.

% Calls a denoiser and return the cleaned image f_est.
% In this example the chosen denoising method is TNRD.
% Note: The devision and multiplication by the constant 5 is done to handle
% a different noise-level than a fixed one (which is equal to 5)

function f_est = Denoiser(x_est,effective_sigma)

f_est = ReactionDiffusion(5/effective_sigma*x_est);
f_est = f_est*effective_sigma/5;

return

author: google

repo: RED

save_path: github-repos/MATLAB/google-RED


text:
% SPM5 UPDATE 17/12/06
% UPDATE 27/01/05
% Sets the default values for the FieldMap toolbox
%
% FORMAT pm_defaults_Allegra_eFoV
%_______________________________________________________________________
%
% This file is intended for use with the Siemens fieldmap sequence
% on the Allegra scanner at the FIL and the new EPI sequence with
% extended  FOV, and PE blips=-1:
% nw_mepi_v3d_efov, nw_mepi_v3d_Amy,
% nw_mepi_v3d_Hippoc, nw_mepi_v3d_OFC_Amy.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Chloe Hutton and Jesper Andersson
% $Id: pm_defaults_Allegra_eFoV.m 5015 2012-10-24 13:40:07Z guillaume $

global pm_def;

% Defaults for creating field map. (See pm_make_fieldmap.m and
%                                   FieldMap.man for more info.)
%=======================================================================
pm_def.INPUT_DATA_FORMAT = 'PM';      % 'RI' = load two real and
                                      % imaginary image pairs
                                      % 'PM' = load one or two
                                      % phase and magnitude image
                                      % pairs.
pm_def.SHORT_ECHO_TIME = 10.0;        % Short echo time in ms for Allegra
pm_def.LONG_ECHO_TIME = 12.46;        % Long echo time in ms for Allegra
pm_def.MASKBRAIN = 1;                 % Do brain masking (1 or 0,
                      % 0 for EPI fieldmaps)

% Defaults for unwrapping options. (See pm_make_fieldmap.m and
%                                   FieldMap.man for more info.)
%=======================================================================
pm_def.UNWRAPPING_METHOD = 'Mark3D';  % Unwrapping options are:
                                      % 'Huttonish', 'Mark3D' or 'Mark2D'
pm_def.FWHM = 10;                     % FWHM of Gaussian filter used to
                                      % implement weighted smoothing of
                                      % unwrapped maps.
pm_def.PAD = 0;                       % Size of padding kernel if required.
pm_def.WS = 1;                        % Weighted or normal smoothing.

% Flags for brain extraction
%=======================================================================
pm_def.MFLAGS.TEMPLATE = fullfile(spm('Dir'),'toolbox','FieldMap','T1.nii');
pm_def.MFLAGS.FWHM = 5;     % In mm
pm_def.MFLAGS.NERODE = 2;   % In voxels
pm_def.MFLAGS.NDILATE = 4;  % In voxels
pm_def.MFLAGS.THRESH = 0.5;
pm_def.MFLAGS.REG = 0.02;   % A larger value helps segmentation to converge
pm_def.MFLAGS.GRAPHICS = 0; % A larger value helps segmentation to converge

% Defaults for converting field map to voxel displacement map.
%=======================================================================
pm_def.EPI_BASED_FIELDMAPS = 0;         % EPI=1, other=0.
pm_def.K_SPACE_TRAVERSAL_BLIP_DIR = -1; % +ve k-space = 1, -ve = -1.
pm_def.TOTAL_EPI_READOUT_TIME = 23.76;   % Allegra EPI RO time (330E-6*72)

% Defaults for Unwarping.
%=======================================================================
pm_def.DO_JACOBIAN_MODULATION = 0;    % Do jacobian modulation to adjust
                                      % for compression or stretching
                                      % No = 0, Yes = 1

author: spm


text:
%% Solving LQR via SDP (LMIs)
% http://www1.se.cuhk.edu.hk/~zhang/Courses/Seg5660/Lecture%206.pdf

%% Continuous LQR
clc
Gs = tf(2,[0.5 0.2 0.1]);
[A,B] = ssdata(ss(Gs));
Q = diag([10;10]); R = 1;

[K,S] = lqr(A,B,Q,R)

%LQR LMI - Vars are P11 P12 P22
P = sym('[P11 P12; P12 P22]')
Lp = vpa([R B'*P; P*B Q + A'*P + P*A],5)

% Objective (max)
f = -[1 0 1]; %just diagonal elements

% Semidefinite Constraint (cheap way to pick off elements)
C = -double(subs(Lp,P,zeros(2)));
A1 = double(diff(Lp,P(1)));
A2 = double(diff(Lp,P(2)));
A3 = double(diff(Lp,P(end)));
sdcone = sparse([C(:) A1(:) A2(:) A3(:)]);

% Create OPTI Object
Opt = opti('f',f,'sdcone',sdcone)
% Solve the SDP
x = solve(Opt)

%Reform P and solve K
P = [x(1) x(2); x(2) x(3)];
Ksdp = R\(B'*P) %see help lqr

%% Discrete LQR
clc
Gd = c2d(Gs,0.1);
[A,B] = ssdata(ss(Gd));
Q = diag([10;10]); R = 1;

[K,S] = dlqr(A,B,Q,R)

%LQR LMI - Vars are P11 P12 P22
P = sym('P',[2 2]);
Lp = vpa([B'*P*B+R B'*P*A; A'*P*B Q + A'*P*A - P],5)

% Objective (max)
f = -[1 0 1]; %just diagonal elements

% Semidefinite Constraint (cheap way to pick off elements)
C = -double(subs(Lp,P,zeros(2)));
A1 = double(diff(Lp,P(1)));
A2 = double(diff(Lp,P(2))) + double(diff(Lp,P(3)));
A3 = double(diff(Lp,P(end)));

sdcone = sparse([C(:) A1(:) A2(:) A3(:)]);

% Create OPTI Object
Opt = opti('f',f,'sdcone',sdcone)
% Solve the SDP
x = solve(Opt)

%Reform P and solve K
P = [x(1) x(2); x(2) x(3)];
Ksdp = (B'*P*B + R)\(B'*P*A) %see help dlqr

%% Continuous Lyapunov
clc
[A,B] = ssdata(ss(Gs));
X = lyap(A,Q)

%LYAP LMI - Vars are P11 P12 P22
P = sym('P',[2 2]);
Lp = vpa(A*P + P*A' + Q)

% Objective (max)
f = -[1 1 1]; %just diagonal elements

% Semidefinite Constraint (cheap way to pick off elements)
C = -double(subs(Lp,P,zeros(2)));
A1 = double(diff(Lp,P(1)));
A2 = double(diff(Lp,P(2))) + double(diff(Lp,P(3)));
A3 = double(diff(Lp,P(end)));

sdcone = sparse([C(:) A1(:) A2(:) A3(:)]);

% Create OPTI Object
Opt = opti('f',f,'sdcone',sdcone,'solver','dsdp')
% Solve the SDP
[x,f,e,i] = solve(Opt)
P = [x(1) x(2); x(2) x(3)];

A*P + P*A' + Q
A*X + X*A' + Q

%% Discrete Lyapunov
clc
[A,B] = ssdata(ss(Gd));
X = dlyap(A,Q)

%LYAP LMI - Vars are P11 P12 P22
P = sym('P',[2 2]);
Lp = vpa(A*P*A' - P + Q)

% Objective (max)
f = -[1 1 1]; %just diagonal elements

% Semidefinite Constraint (cheap way to pick off elements)
C = -double(subs(Lp,P,zeros(2)));
A1 = double(diff(Lp,P(1)));
A2 = double(diff(Lp,P(2))) + double(diff(Lp,P(3)));
A3 = double(diff(Lp,P(end)));

sdcone = sparse([C(:) A1(:) A2(:) A3(:)]);

% Create OPTI Object
Opt = opti('f',f,'sdcone',sdcone,'solver','dsdp')
% Solve the SDP
[x,f,e,i] = solve(Opt)
P = [x(1) x(2); x(2) x(3)];

A*P*A' - P + Q
A*X*A' - X + Q


text:
function [tt]=tt_qshift_p(d,m)

% for integer m
% returns the periodic downward m-position shift matrix of size 2^d x 2^d
% in the QTT format
% The shift is (-m)-position upward for m < 0
% m=0 corresponds to the identity matrix
% m=2^d corresponds to the identity matrix
% m=-2^d corresponds to the identity matrix
%
% January 28, 2011
% Vladimir Kazeev
% vladimir.kazeev@gmail.com
% INM RAS
% Moscow, Russia
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For details please see the preprint
% http://www.mis.mpg.de/publications/preprints/2011/prepr2011-36.html
% Vladimir A. Kazeev, Boris N. Khoromskij and Eugene E. Tyrtyshnikov
% Multilevel Toeplitz matrices generated by QTT tensor-structured vectors and convolution with logarithmic complexity
% January 12, 2012
% Vladimir Kazeev,
% Seminar for Applied Mathematics, ETH Zurich
% vladimir.kazeev@sam.math.ethz.ch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I=eye(2);
J=[0,1;0,0];
O=zeros(2);
P=J+J';

m=mod(m,2^d);
m=dec2bin(m,d);

tt=cell(d,1);

if (m(1) == '1')
	tt{1}=[P,I];
else
	tt{1}=[I,P];
end
tt{1}=reshape(tt{1},[2,2,2]);

for k=2:d-1
	if (m(k) == '1')
		tt{k}=[J',O;J,I];
	else
		tt{k}=[I,J';O,J];
	end
	tt{k}=permute(reshape(tt{k},[2,2,2,2]),[1,3,2,4]);
end

if (m(d) == '1')
	tt{d}=[J',J];
else
	tt{d}=[I,O];
end
tt{d}=reshape(tt{d},[2,2,2]);

tt=tt_reverse(tt,2);

return
end

author: oseledets

repo: TT-Toolbox


text:
function problem = getDefaultOptions(problem)
% problem = getDefaultOptions(problem)
%
% This function fills in any blank entries in the problem.options struct.
% It is designed to be called from inside of optimTraj.m, and not by the
% user.
%

%%%% Top-level default options:
OPT.method = 'trapezoid';
OPT.verbose = 2;
OPT.defaultAccuracy = 'medium';

%%%% Basic setup

% ensure that options is not empty
if ~isfield(problem,'options')
    problem.options.method = OPT.method;
end
opt = problem.options;

% Loop over each options struct and fill in top-level options
for i=1:length(opt)
    if ~isfield(opt(i),'method')
        opt(i).method = OPT.method;
    elseif isempty(opt(i).method)
        opt(i).method = OPT.method;
    end
    if ~isfield(opt(i),'verbose')
        opt(i).verbose = OPT.verbose;
    elseif isempty(opt(i).verbose)
        opt(i).verbose = OPT.verbose;
    end
    if ~isfield(opt(i),'defaultAccuracy')
        opt(i).defaultAccuracy = OPT.defaultAccuracy;
    elseif isempty(opt(i).defaultAccuracy)
        opt(i).defaultAccuracy = OPT.defaultAccuracy;
    end
end

% Figure out basic problem size:
nState = size(problem.guess.state,1);
nControl = size(problem.guess.control,1);

% Loop over opt and fill in nlpOpt struct:
for i=1:length(opt)
    switch opt(i).verbose
        case 0
            NLP_display = 'notify';
        case 1
            NLP_display = 'final-detailed';
        case 2
            NLP_display = 'iter';
        case 3
            NLP_display = 'iter-detailed';
        otherwise
            error('Invalid value for options.verbose');
    end
    switch opt(i).defaultAccuracy
        case 'low'
            OPT.nlpOpt = optimset(...
                'Display',NLP_display,...
                'TolFun',1e-4,...
                'MaxIter',200,...
                'MaxFunEvals',1e4*(nState+nControl));
        case 'medium'
            OPT.nlpOpt = optimset(...
                'Display',NLP_display,...
                'TolFun',1e-6,...
                'MaxIter',400,...
                'MaxFunEvals',5e4*(nState+nControl));
        case 'high'
            OPT.nlpOpt = optimset(...
                'Display',NLP_display,...
                'TolFun',1e-8,...
                'MaxIter',800,...
                'MaxFunEvals',1e5*(nState+nControl));
        otherwise
            error('Invalid value for options.defaultAccuracy')
    end
    if isfield(opt(i),'nlpOpt')
        if isstruct(opt(i).nlpOpt) && ~isempty(opt(i).nlpOpt)
            names = fieldnames(opt(i).nlpOpt);
            for j=1:length(names)
                if ~isfield(OPT.nlpOpt,names{j})
                    disp(['WARNING: options.nlpOpt.' names{j} ' is not a valid option']);
                else
                    OPT.nlpOpt.(names{j}) = opt(i).nlpOpt.(names{j});
                end
            end
        end
    end
    opt(i).nlpOpt = OPT.nlpOpt;
end

% Check ChebFun dependency:
missingChebFun = false;
for i=1:length(opt)
if strcmp(opt(i).method,'chebyshev')
try
    chebpts(3);  %Test call to chebfun
catch ME %#ok<NASGU>
    missingChebFun = true;
    opt(i).method = 'trapezoid';  %Force default method
end
end
end
if missingChebFun
   warning('''chebyshev'' method requires the Chebfun toolbox');
   disp('   --> Install Chebfun toolbox:  (http://www.chebfun.org/)');
   disp('   --> Running with default method instead (''trapezoid'')');
end

% Fill in method-specific paramters:
for i=1:length(opt)
    OPT_method = opt(i).method;
    switch OPT_method
        case 'trapezoid'
            OPT.trapezoid = defaults_trapezoid(opt(i).defaultAccuracy);
        case 'hermiteSimpson'
            OPT.hermiteSimpson = defaults_hermiteSimpson(opt(i).defaultAccuracy);
        case 'chebyshev'
            OPT.chebyshev = defaults_chebyshev(opt(i).defaultAccuracy);
        case 'multiCheb'
            OPT.multiCheb = defaults_multiCheb(opt(i).defaultAccuracy);
        case 'rungeKutta'
            OPT.rungeKutta = defaults_rungeKutta(opt(i).defaultAccuracy);
        case 'gpops'
            OPT.gpops = defaults_gpops(opt(i).defaultAccuracy);
        otherwise
            error('Invalid value for options.method');
    end
    if isfield(opt(i),OPT_method)
        if isstruct(opt(i).(OPT_method)) && ~isempty(opt(i).(OPT_method))
            names = fieldnames(opt(i).(OPT_method));
            for j=1:length(names)
                if ~isfield(OPT.(OPT_method),names{j})
                    disp(['WARNING: options.' OPT_method '.' names{j} ' is not a valid option']);
                else
                    OPT.(OPT_method).(names{j}) = opt(i).(OPT_method).(names{j});
                end
            end
        end
    end
    opt(i).(OPT_method) = OPT.(OPT_method);
end

problem.options = opt;
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                    Method-specific parameters                           %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

function OPT_trapezoid = defaults_trapezoid(accuracy)

switch accuracy
    case 'low'
        OPT_trapezoid.nGrid = 12;
    case 'medium'
        OPT_trapezoid.nGrid = 30;
    case 'high'
        OPT_trapezoid.nGrid = 60;
    otherwise
        error('Invalid value for options.defaultAccuracy')
end

OPT_trapezoid.adaptiveDerivativeCheck = 'off';

end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

function OPT_hermiteSimpson = defaults_hermiteSimpson(accuracy)

switch accuracy
    case 'low'
        OPT_hermiteSimpson.nSegment = 10;
    case 'medium'
        OPT_hermiteSimpson.nSegment = 20;
    case 'high'
        OPT_hermiteSimpson.nSegment = 40;
    otherwise
        error('Invalid value for options.defaultAccuracy')
end

OPT_hermiteSimpson.adaptiveDerivativeCheck = 'off';

end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

function OPT_chebyshev = defaults_chebyshev(accuracy)

switch accuracy
    case 'low'
        OPT_chebyshev.nColPts = 9;
    case 'medium'
        OPT_chebyshev.nColPts = 13;
    case 'high'
        OPT_chebyshev.nColPts = 23;
    otherwise
        error('Invalid value for options.defaultAccuracy')
end

end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

function OPT_multiCheb = defaults_multiCheb(accuracy)

switch accuracy
    case 'low'
        OPT_multiCheb.nColPts = 6;
        OPT_multiCheb.nSegment = 3;
    case 'medium'
        OPT_multiCheb.nColPts = 8;
        OPT_multiCheb.nSegment = 6;
    case 'high'
        OPT_multiCheb.nColPts = 8;
        OPT_multiCheb.nSegment = 12;
    otherwise
        error('Invalid value for options.defaultAccuracy')
end

end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

function OPT_rungeKutta = defaults_rungeKutta(accuracy)

switch accuracy
    case 'low'
        OPT_rungeKutta.nSegment = 10;
        OPT_rungeKutta.nSubStep = 2;
    case 'medium'
        OPT_rungeKutta.nSegment = 20;
        OPT_rungeKutta.nSubStep = 2;
    case 'high'
        OPT_rungeKutta.nSegment = 20;
        OPT_rungeKutta.nSubStep = 4;
    otherwise
        error('Invalid value for options.defaultAccuracy')
end

OPT_rungeKutta.adaptiveDerivativeCheck = 'off';

end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

function OPT_gpops = defaults_gpops(accuracy)

OPT_gpops.bounds.phase.integral.lower = -inf;
OPT_gpops.bounds.phase.integral.upper = inf;
OPT_gpops.guess.phase.integral = 0;

OPT_gpops.name = 'OptimTraj_GPOPS';
OPT_gpops.auxdata = [];
OPT_gpops.nlp.solver = 'ipopt'; % {'ipopt','snopt'}
OPT_gpops.derivatives.dependencies = 'full';  %�full�, �sparse� or �sparseNaN�
OPT_gpops.derivatives.supplier = 'sparseCD'; %'sparseCD';  %'adigator'
OPT_gpops.derivatives.derivativelevel = 'first'; %'second';
OPT_gpops.mesh.method = 'hp-PattersonRao';
OPT_gpops.method = 'RPM-Integration';
OPT_gpops.mesh.phase.colpoints = 10*ones(1,10);
OPT_gpops.mesh.phase.fraction = ones(1,10)/10;
OPT_gpops.scales.method = 'none'; % { 'none' , automatic-hybridUpdate' , 'automatic-bounds';

switch accuracy
    case 'low'
        OPT_gpops.mesh.tolerance = 1e-2;
        OPT_gpops.mesh.maxiterations = 0;
    case 'medium'
        OPT_gpops.mesh.tolerance = 1e-3;
        OPT_gpops.mesh.maxiterations = 1;
    case 'high'
        OPT_gpops.mesh.tolerance = 1e-4;
        OPT_gpops.mesh.maxiterations = 3;
    otherwise
        error('Invalid value for options.defaultAccuracy')
end

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function x = emailFeatures(word_indices)
%EMAILFEATURES takes in a word_indices vector and produces a feature vector
%from the word indices
%   x = EMAILFEATURES(word_indices) takes in a word_indices vector and
%   produces a feature vector from the word indices.

% Total number of words in the dictionary
n = 1899;

% You need to return the following variables correctly.
x = zeros(n, 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return a feature vector for the
%               given email (word_indices). To help make it easier to
%               process the emails, we have have already pre-processed each
%               email and converted each word in the email into an index in
%               a fixed dictionary (of 1899 words). The variable
%               word_indices contains the list of indices of the words
%               which occur in one email.
%
%               Concretely, if an email has the text:
%
%                  The quick brown fox jumped over the lazy dog.
%
%               Then, the word_indices vector for this text might look
%               like:
%
%                   60  100   33   44   10     53  60  58   5
%
%               where, we have mapped each word onto a number, for example:
%
%                   the   -- 60
%                   quick -- 100
%                   ...
%
%              (note: the above numbers are just an example and are not the
%               actual mappings).
%
%              Your task is take one such word_indices vector and construct
%              a binary feature vector that indicates whether a particular
%              word occurs in the email. That is, x(i) = 1 when word i
%              is present in the email. Concretely, if the word 'the' (say,
%              index 60) appears in the email, then x(60) = 1. The feature
%              vector should look like:
%
%              x = [ 0 0 0 0 1 0 0 0 ... 0 0 0 0 1 ... 0 0 0 1 0 ..];
%
%

for loc = word_indices,
  x(loc) = 1;
endfor

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

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
function p09_story ( )

%*****************************************************************************80
%
%% P09_STORY prints the "story" for problem 09
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This example is due to Max Waldmeier.\n' );
  fprintf ( 1, '  This data represents a measure of sunspot activity over the\n' );
  fprintf ( 1, '  years 1700 to 1960.  The X value is the year, and the Y value\n' );
  fprintf ( 1, '  is a measure of the sunspot activity, which is usually, but\n' );
  fprintf ( 1, '  not always, an integer.\n' );

  return
end


text:
function A = shiftAnglesFromMinusPIToPI(InA)
    %dont do anything to angles that are already in range
    %if we leave out this step, we get +180 shifted to -180
    idx = find(InA>pi | InA <-pi);
    if(isempty(idx))
        A = InA;
    else
        A = InA;
        A(idx) = mod(InA(idx)+pi,2*pi)-pi;  % input is in radians
    end
end


text:
function a = c8mat_exp_a ( test, n )

%*****************************************************************************80
%
%% C8MAT_EXP_A returns the matrix for a given test.
%
%  Discussion:
%
%    1) Diagonal, real;
%    2) Diagonal, pure imaginary.
%    3) Diagonal, complex.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TEST, the index of the test case.
%
%    Input, integer N, the order of the matrix.
%
%    Output, real A(N,N), the matrix.
%
  if ( test == 1 )
    a = [ 1,  0;
          0,  2 ];
  elseif ( test == 2 )
    a = [  3 * i,  0;
           0,   - 4 * i ];
  elseif ( test == 3 )
    a = [  5 + 6 * i,    0;
          0,  7 - 8 * i ];
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'C8MAT_EXP_A - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of TEST = %d\n', test );
    error ( 'C8MAT_EXP_A - Fatal error!' );
  end

  return
end


text:
function y = daub4_scale ( n, x )

%*****************************************************************************80
%
%% DAUB4_SCALE recursively evaluates the DAUB4 scaling function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the recursion level.
%
%    Input, real X, the point at which the function is to be evaluated.
%
%    Output, real Y, the estimated value of the function.
%
  c = [  0.4829629131445341E+00; ...
         0.8365163037378079E+00; ...
         0.2241438680420133E+00; ...
        -0.1294095225512603E+00 ];

  c = sqrt ( 2.0 ) * c;

  if ( 0 < n )
    y = ( c(1) * daub4_scale ( n - 1, 2 * x     ) ...
        + c(2) * daub4_scale ( n - 1, 2 * x - 1 ) ...
        + c(3) * daub4_scale ( n - 1, 2 * x - 2 ) ...
        + c(4) * daub4_scale ( n - 1, 2 * x - 3 ) );
  else
    y = ( 0 <= x & x < 1 );
  end

  return
end


text:
%% Part of discussion on ADC quadratic form calculations
% Probably old between Bob and Brian years ago.
%
% Could be fixed up or deprecated.
%

% Here is a standard ellipsoid.
[X,Y,Z] = ellipsoid(0,0,0,3,1,1);
figure(1), subplot(1,2,1), surf(X,Y,Z)
axis equal

sz = size(X);

% Here is how an ellipse looks like a peanut
%
% Assign a length equal to the square of their length on the ellipsoid.
% This corresponds to what I believe are the ADC (variance) values.  Notice
% the peanut shape.

% The length of each ellipsoid vector
tmp = [X(:),Y(:),Z(:)];
l = sqrt(diag(tmp*tmp'));
% Multiply each vector by its own length
tmp = diag(l)*tmp;

X2 = reshape(tmp(:,1),sz);
Y2 = reshape(tmp(:,2),sz);
Z2 = reshape(tmp(:,3),sz);
figure(1), subplot(1,2,2), surf(X2,Y2,Z2)

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ a_lu, info ] = r83_np_fa ( n, a )

%*****************************************************************************80
%
%% R83_NP_FA factors a R83 matrix without pivoting.
%
%  Discussion:
%
%    The R83 storage format is used for a tridiagonal matrix.
%    The superdiagonal is stored in entries (1,2:N), the diagonal in
%    entries (2,1:N), and the subdiagonal in (3,1:N-1).  Thus, the
%    original matrix is "collapsed" vertically into the array.
%
%    Because this routine does not use pivoting, it can fail even when
%    the matrix is not singular, and it is liable to make larger
%    errors.
%
%    R83_NP_FA and R83_NP_SL may be preferable to the corresponding
%    LINPACK routine SGTSL for tridiagonal systems, which factors and solves
%    in one step, and does not save the factorization.
%
%  Example:
%
%    Here is how a R83 matrix of order 5 would be stored:
%
%       *  A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54  *
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be at least 2.
%
%    Input, real A(3,N), the tridiagonal matrix.
%
%    Output, integer INFO, singularity flag.
%    0, no singularity detected.
%    nonzero, the factorization failed on the INFO-th step.
%
%    Output, real A_LU(3,N), factorization information.
%
  info = 0;

  a_lu(1:3,1:n) = a(1:3,1:n);

  for i = 1 : n-1

    if ( a_lu(2,i) == 0.0E+00 )
      info = i;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R83_NP_FA - Fatal error!\n' );
      fprintf ( 1, '  Zero pivot on step %d\n', info );
      return;
    end
%
%  Store the multiplier in L.
%
    a_lu(3,i) = a_lu(3,i) / a_lu(2,i);
%
%  Modify the diagonal entry in the next column.
%
    a_lu(2,i+1) = a_lu(2,i+1) - a_lu(3,i) * a_lu(1,i+1);

  end

  if ( a_lu(2,n) == 0.0E+00 )
    info = n;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R83_NP_FA - Fatal error!\n' );
    fprintf ( 1, '  Zero pivot on step %d\n', info );
    return;
  end

  return
end


text:
function ht = fix_xticklabels(handle,margins,textopts)

%FIX_XTICKLABELS will determine maximum allowed width
%   of long XTickLabels and convert them into multiple line
%   format to fit maximum allowed width
%
%
%   IN:
%       handle      [optional] - axes handle (defaults to gca)
%       margins     [optional] - relative margins width (defaults to 0.1)
%       texopts     [optional] - cell array of addtitional text formating
%                                options (defaults to {} )
%
%   OUT:
%       ht - column vector of handles for created text objects
%
% This code was derrived from MY_XTICKLABELS by Pekka Kumpulainen
% http://www.mathworks.com/matlabcentral/fileexchange/19059-myxticklabels
%   EXAMPLE:
%       figure;
%
%       bar(rand(1,4));
%       set(gca,'XTickLabel',{'Really very very long text', ...
%                             'One more long long long text', ...
%                             'Short one','Long long long text again'});
%       fix_xticklabels();
%
%       figure;
%
%       bar(rand(1,4));
%       set(gca,'XTickLabel',{'Really very very long text', ...
%                             'One more long long long text', ...
%                             'Short one','Long long long text again'});
%       fix_xticklabels(gca,0.1,{'FontSize',16});
%

% Mikhail Erofeev 07.06.2013
% Pekka Kumpulainen 12.2.2008

if(~exist('handle','var'))
    handle = gca;
end

if(~exist('textopts','var'))
   textopts={};
end
if(~exist('margins','var'))
    margins = 0.1;
end

xtickpos = get(handle,'XTick');
xtickstring = get(handle,'XTickLabel');

set(handle, 'XTickLabel','')
h_olds = findobj(handle, 'Tag', 'MUXTL');
if ~isempty(h_olds)
    delete(h_olds)
end

%% Make XTickLabels
NTick = length(xtickpos);
Ybot = min(get(handle,'YLim'));

Xbot = get(handle,'XLim');
max_w = ((max(Xbot) - min(Xbot) )/NTick)*(1-margins);

ht = zeros(NTick,1);

NLabels = length(xtickstring);

for ii = 1:NTick
    str_id = mod( ii - 1 , NLabels) + 1;
    ht(ii) = text('String',xtickstring{str_id}, ...
        'Units','data', ...
        'VerticalAlignment', 'top', ...
        'HorizontalAlignment', 'center ', ...
        'Position',[xtickpos(ii) Ybot], ...
        'Tag','MUXTL');
end
if ~isempty(textopts)
    set(ht,textopts{:})
end

for ii = 1:NTick
    auto_width(ht(ii),max_w);
end

%% squeeze axis if needed

set(handle,'Units','pixels')
Axpos = get(handle,'Position');
% set(Hfig,'Units','pixels')
% Figpos = get(Hfig,'Position');

set(ht,'Units','pixels')
TickExt = zeros(NTick,4);
for ii = 1:NTick
    TickExt(ii,:) = get(ht(ii),'Extent');
end

needmove = -(Axpos(2) + min(TickExt(:,2)));

if needmove>0;
    Axpos(2) = Axpos(2)+needmove+2;
    Axpos(4) = Axpos(4)-needmove+2;
    set(handle,'Position',Axpos);
end

set(handle,'Units','normalized')
set(ht,'Units','normalized')
end

function []=auto_width(h,maxw)
    str = get(h,'string');
    set(h,'string',{});
    words = regexp(str,' ','split');

    tail = words(2:end);
    res = words(1);
    ln = 1;

    while(~isempty(tail))
        lin = res{ln};
        lin_e = [lin ' ' tail{1}];
        set(h,'string',lin_e);
        ext = get(h,'Extent');
        width=ext(3);
        if(width>maxw)
            ln =ln+1;
            res{ln} = tail{1};
        else
            res{ln} = lin_e;
        end
        tail = tail(2:end);
    end
    set(h,'string',res);
    ext = get(h,'Extent');
    width=ext(3);
    if(width>maxw)
        warning('FIX_XTICKLABELS:over','Overflown box');
    end
end


text:
function [A] = cm2A(cm)
% Convert length from centimeters to angstroms.
% Chad A. Greene 2012
A = cm*100000000;


text:
% eeg_mergechan() - merge channel structure while preserving channel
%                   order
%
% >> mergelocs = eeg_mergechan(locs1, locs2);
%
% Inputs:
%     locs1     - EEGLAB channel location structure
%     locs2     - second EEGLAB channel location structure
%
% Output:
%     mergelocs - merged channel location structure
%
% Author: Arnaud Delorme, August 2006

% Copyright (C) Arnaud Delorme, CERCO, 2006, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% union of two channel location structure
% without loosing the order information
% ---------------------------------------
function alllocs = myunion(locs1, locs2)

   labs1 = { locs1.labels };
   labs2 = { locs2.labels };

   count1 = 1;
   count2 = 1;
   count3 = 1;
   alllocs = locs1; alllocs(:) = [];
   while count1 <= length(locs1) | count2 <= length(locs2)

       if count1 > length(locs1)
           alllocs(count3) = locs2(count2);
           count2 = count2 + 1;
           count3 = count3 + 1;
       elseif count2 > length(locs2)
           alllocs(count3) = locs1(count1);
           count1 = count1 + 1;
           count3 = count3 + 1;
       elseif strcmpi(labs1{count1}, labs2{count2})
           alllocs(count3) = locs1(count1);
           count1 = count1 + 1;
           count2 = count2 + 1;
           count3 = count3 + 1;
       elseif isempty(strmatch(labs1{count1}, labs2))
           alllocs(count3) = locs1(count1);
           count1 = count1 + 1;
           count3 = count3 + 1;
       else
           alllocs(count3) = locs2(count2);
           count2 = count2 + 1;
           count3 = count3 + 1;
       end;

   end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [nodes, edges, grayPath] = loadGrayNodes(hemisphere, grayPath, numLayers)
%
% [nodes, edges, grayPath] = loadGrayNodes(hemisphere, [grayPath], [numLayers])
%
% Loads gray nodes from specified gray class file.
% If the file is not fully specified, opens dialog box to get it from the user.
%
% hemisphere: 'right' or 'left'
% grayPath: path string to file
%    If not specified, opens dialog box to have user find the file below:
%       /usr/local/mri/anatomy/<subject>/<hemisphere>
%
% djh, 8/98
%
% 022699 - Get anatomy path from getAnatomyPath - WAP
% 061899	- Browse with uigetfile - arw
% 080299 - fixed bug with returned restofpath,
%          change into correct directory before browsing - arw
% 080800 - allowed returning null nodes and edges if restofpath is sent in null -- djh, bw
% 120100 - fixed bug in which when restofPath was sent in it was never
%          appended to grayPath.  BW
% 120100 - AW reorganised the logic a little. Exits immediately if restOfPath
%          is null.
% 2/5/2001, djh - cleaned this up (major rewrite)
%          wrote getPathStrDialog.m to get the code for that outta this function
%          dumped restOfPath, replaced it with full grayPath
% 2008.01.?? RFD: added code to automatically grow gray matter in given a class file.
% 2008.02.04 RFD: when growing from a class file, we now prompt for the number of layers
%            if it isn't provided, rather than using a hard-coded default.

global mrSESSION

% Open dialog box to find gray class file
if notDefined('grayPath')

    % default initial directory: anatomy dir / hemisphere
    Hemi = hemisphere; Hemi(1) = upper(Hemi(1)); % e.g., left -> Left
    startDir = fullfile(getAnatomyPath(mrSESSION.subject), Hemi);
    % if 'Left' dir does not exist go back to top dir
    if ~exist(startDir,'dir')
	  startDir = getAnatomyPath(mrSESSION.subject);
    end

    filterspec = {'*.nii.gz;*.nii', 'NIFTI class file'; '*.*', 'All files'};
    grayPath = getPathStrDialog(startDir, 'Select class file', filterspec);
end

if exist(grayPath, 'file')
    disp(['Building gray graph from: ', grayPath]);
    class = readClassFile(grayPath, 0, 0, hemisphere);
    if(~exist('numLayers','var')||isempty(numLayers))
        a = inputdlg('Number of gray layers:','Gray Layers',1,{'5'});
        if(isempty(a)), error('user canceled.'); end
        numLayers = str2double(a{1});
    end
    [nodes,edges] = mrgGrowGray(class,numLayers);
elseif isempty(grayPath)
    fprintf('No file selected, so no coords loaded. \n');
    nodes = [];
    edges = [];

else
    myWarnDlg(['File ',grayPath,' does not exist. No gray nodes loaded']);
    nodes=[];
    edges=[];
end

return;


text:
%% DEMO 03: Generate sample data and add realistic CT noise to it.
%
% This demo will show how to generate sample data for image reconstruction
%
%
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------
%% Initialize
clear;
close all;
%% Geometry
geo=defaultGeometry();
%% Define angles of projection and load phatom image

% define projection angles (in radians)
angles=linspace(0,2*pi,100);
% load phatnom image
head=headPhantom(geo.nVoxel);

% Simulate forward projection.
% Strongly suggested to use 'iterpolated' option for more accurate
% projections. reduce geo.accuracy for better results
projections=Ax(head,geo,angles,'interpolated');

% Add realistic noise. Adds photon scattering noise ('Poisson') and
% electronic noise of the detector ('Gaussian').
%
% 'Poisson' is related to the maximum photon count in the detector. 1e5 is
% a standard clinical nuber, reduce it for more noise
% 'Gaussian' is related to possible electronic noise in the detector. mean
% of 0 and std of 10 is common in clinical scenario. Increase std for more
% noise.
noise_projections=addCTnoise(projections,'Poisson',1e5,'Gaussian',[0 10]);

% Plot Projections
plotProj(projections,angles)
% plot noise
plotProj(projections-noise_projections,angles)

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
%% (Internal) Guess PPG/ABP signals indexes in a multimodal recording
%
%   [PPG_ABP_idx, ECG_header ] = get_PPG_ABP_idx_from_header(ECG_header)
%
% Arguments:
%
%      + ECG_header: header info
%
% Output:
%
%      + PPG_ABP_idx: indexes of the ECG signals in the multimodal matrix
%
%      + ECG_header: header of a ECG only recording
%
% Example:
%
% See also get_ECG_idx_from_header, ECGwrapper
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function [PPG_ABP_idx, ECG_header ] = get_PPG_ABP_idx_from_header(ECG_header)

    PPG_ABP_idx = [];

    str_12_leads_desc = { 'PPG' 'PLETH' 'ABP' 'ART' 'BP' 'PULSEWAVE' 'SPO2' };
    % add ECG patterns to be matched in ECG_header.desc
    PPG_ABP_pattern_desc = {'PLETH' 'ABP' 'ART' 'BP' 'PPG' };

    [~, ~, aux_idx ] = intersect(str_12_leads_desc, upper(strtrim(cellstr(ECG_header.desc))) );

    if( ECG_header.nsig ~= length(aux_idx) )

        PPG_ABP_idx = colvec(aux_idx);

        aux_idx = find(any(cell2mat(cellfun(@(b)(cell2mat(cellfun(@(a)(~isempty(strfind(a, b))), rowvec(upper(cellstr(ECG_header.desc))), 'UniformOutput', false))), colvec(PPG_ABP_pattern_desc), 'UniformOutput', false)),1));

        if( isempty(aux_idx) && isempty(PPG_ABP_idx) )
            if( isempty(PPG_ABP_idx) )
                cprintf('[1,0.5,0]', disp_option_enumeration('Could not find any PPG/ABP signal, check the lead description of the recording:', cellstr(ECG_header.desc) ) );
                fprintf(1, '\n')
            end
        else
            PPG_ABP_idx = unique([PPG_ABP_idx; colvec(aux_idx)]);
            Not_ECG_idx = setdiff(1:ECG_header.nsig, PPG_ABP_idx);
            if( ~isempty(Not_ECG_idx) )
%                 warning('get_ECG_idx_from_header:SignalsDiscarded', disp_option_enumeration('Some signal/s present are not ECG:', cellstr(ECG_header.desc(Not_ECG_idx,:)) ))
            end
        end

    else
        % all leads are standard ECG leads
        PPG_ABP_idx = 1:ECG_header.nsig;
    end

    if(length(PPG_ABP_idx) ~= ECG_header.nsig && nargout > 1)
       %trim the heasig if needed

        ECG_header = trim_ECG_header(ECG_header, PPG_ABP_idx);

    end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function y = modu(x,dzdy)
%ModU activation function
%Ye, C., Yang, Y., Fermuller, C., & Aloimonos, Y. (2017).
%On the Importance of Consistency in Training Deep Neural Networks. arXiv preprint arXiv:1708.00631.

  if nargin <= 1 || isempty(dzdy)
    y = abs(x) ;
  else
    y = dzdy .* sign(x) ;
  end

end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
% -*- INTERNAL UNDOCUMENTED FUNCTION -*-
function [geometry, msh, sp, u, gnum] = ...
                    mp_solve_linear_elasticity_3d (problem_data, method_data)

warning ('geopdes:obsolete', 'Function MP_SOLVE_LINEAR_ELASTICITY_3D is obsolete. Using MP_SOLVE_LINEAR_ELASTICITY instead')

[geometry, msh, sp, u, gnum] = mp_solve_linear_elasticity (problem_data, method_data);

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function llh=local2llh(xy,origin)
%local2llh     llh=local2llh(xy,origin)
%
%Converts from local coorindates to longitude and latitude
%given the [lon, lat] of an origin. 'origin' should be in
%decimal degrees. Note that heights are ignored and that
%xy is in km.  Output is [lon, lat, height] in decimal
%degrees. This is an iterative solution for the inverse of
%a polyconic projection.

%-------------------------------------------------------------
%   Record of revisions:
%
%   Date          Programmer            Description of Change
%   ====          ==========            =====================
%
%   Aug 23, 2001  Jessica Murray        Clarification to help.
%
%   Apr 4, 2001   Peter Cervelli        Added failsafe to avoid
%                                       infinite loop because of
%                                       covergence failure.
%   Sep 7, 2000   Peter Cervelli		Original Code
%
%-------------------------------------------------------------

%Set ellipsoid constants (WGS84)

   a=6378137.0;
   e=0.08209443794970;

%Convert to radians / meters

%   xy=xy*1000;
%   origin=origin*pi/180;
   xy=double(xy)*1000;
   origin=double(origin)*pi/180;

%Iterate to perform inverse projection

   M0=a*((1-e^2/4-3*e^4/64-5*e^6/256)*origin(2) - ...
        (3*e^2/8+3*e^4/32+45*e^6/1024)*sin(2*origin(2)) + ...
        (15*e^4/256 +45*e^6/1024)*sin(4*origin(2)) - ...
        (35*e^6/3072)*sin(6*origin(2)));

   z=xy(2,:)~=-M0;

   A=(M0+xy(2,z))/a;
   B=xy(1,z).^2./a^2+A.^2;

   llh(2,z)=A;

   delta=Inf;

   c=0;

   while max(abs(delta))>1e-8

      C=sqrt((1-e^2*sin(llh(2,z)).^2)).*tan(llh(2,z));

      M=a*((1-e^2/4-3*e^4/64-5*e^6/256)*llh(2,z) - ...
           (3*e^2/8+3*e^4/32+45*e^6/1024)*sin(2*llh(2,z)) + ...
           (15*e^4/256 +45*e^6/1024)*sin(4*llh(2,z)) - ...
           (35*e^6/3072)*sin(6*llh(2,z)));

      Mn=1-e^2/4-3*e^4/64-5*e^6/256 - ...
         -2*(3*e^2/8+3*e^4/32+45*e^6/1024)*cos(2*llh(2,z)) + ...
         4*(15*e^4/256 +45*e^6/1024)*cos(4*llh(2,z)) + ...
         -6*(35*e^6/3072)*cos(6*llh(2,z));

      Ma=M/a;

      delta=-(A.*(C.*Ma+1)-Ma-0.5*(Ma.^2+B).*C)./ ...
           (e^2*sin(2*llh(2,z)).*(Ma.^2+B-2*A.*Ma)./(4*C)+(A-Ma).*(C.*Mn-2./sin(2*llh(2,z)))-Mn);

      llh(2,z)=llh(2,z)+delta;

      c=c+1;
      if c>100
          error('Convergence failure.')
      end
   end

   llh(1,z)=(asin(xy(1,z).*C/a))./sin(llh(2,z))+origin(1);

%Handle special case of latitude = 0

   llh(1,~z)=xy(1,~z)/a+origin(1);
   llh(2,~z)=0;

%Convert back to decimal degrees

   llh=llh*180/pi;

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function [volume, permutevec] = volumepermute(volume, permutevec)

% VOLUMEPERMUTE
%
% See also VOLUMEFLIP

% do a low-level check on the input data
if ~isfield(volume, 'transform'), ft_error('the input volume needs a transformation matrix'); end
if ~isfield(volume, 'dim'),       ft_error('the input volume needs a dim field');             end

% define some variable locally
T   = volume.transform;
dim = volume.dim;

% determine which fields can be permuted
fnames = fieldnames(volume);
sel    = false(1,numel(fnames));
for k = 1:numel(fnames)
  tmp = volume.(fnames{k});
  if isnumeric(tmp) && numel(tmp)==prod(dim) && ndims(tmp)==3
    sel(k) = true;
  end
end
fnames = fnames(sel);

if nargin<2
  permutevec = 'auto';
end

if isequal(permutevec, 'auto')
  % determine the order of permutation to make the transformatiom matrix approximately diagonal
  [dum, m1]  = max(abs(T(1,1:3)));
  [dum, m2]  = max(abs(T(2,1:3)));
  % [dum, m3]  = max(abs(T(3,1:3)));
  [dum, m3]  = setdiff(1:3, [m1 m2]); % whichever dimension remains
  permutevec = [m1 m2 m3];
end

if ~all(permutevec==[1 2 3])
  % do the permutation on the numeric data
  for k = 1:numel(fnames)
    volume = setsubfield(volume, fnames{k}, permute(getsubfield(volume, fnames{k}), permutevec));
  end

  % do the permutation on the transformation matrix and dim
  volume.transform = volume.transform(:,[permutevec 4]);
  volume.dim       = volume.dim(permutevec);
else
  % do nothing
end


text:
function r = ge(p,q)
% this function compares the greater than or equality condition of x, y, and z values of the structures that I use in
% the FVtool.
%
% SYNOPSIS:
%
%
% PARAMETERS:
%
%
% RETURNS:
%
%
% EXAMPLE:
%
% SEE ALSO:
%

% Copyright (c) 2012-2016 Ali Akbar Eftekhari
% See the license file

if (isa(p, 'FaceVariable')&&isa(q, 'FaceVariable'))
    r=p;
    r.xvalue = p.xvalue>=q.xvalue;
    r.yvalue = p.yvalue>=q.yvalue;
    r.zvalue = p.zvalue>=q.zvalue;
elseif isa(p, 'FaceVariable')
    r=p;
    r.xvalue = p.xvalue>=q;
    r.yvalue = p.yvalue>=q;
    r.zvalue = p.zvalue>=q;
else
    r=q;
    r.xvalue = p>=q.xvalue;
    r.yvalue = p>=q.yvalue;
    r.zvalue = p>=q.zvalue;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [tetIsoIndex]=tetIsotropyIndex(E,V)
%% function for computing  triangular faces isotropy index.
% This is an index of the  "regularity" of the triangle. Equilateral triangle (perfectly regular)
%
% INPUTS:
% * E: nElements-by-4 array representing the list of vertices for tetrahedral elements
% * V: nVertices-by-3 array representing the 3D positions of the vertices of E
%
% OUTPUTS:
% * tetIsoIndex: isotropy index for each element. value of 1 means
% perfectly isotropic (equilateral tetrahedron). value of 0 means vertices are coplanar.
%
%From the paper: Surface-Marker Cluster Design Criteria for 3-D Bone Movement Reconstruction (1997)
% Aurelio Cappozzo, Angelo Cappello, Ugo Della Croce, and Francesco Pensalfini
%%

tetIsoIndex=zeros(size(E,1),1);
for itet=1:size(E,1)

    x1=V(E(itet,1),:)';
    x2=V(E(itet,2),:)';
    x3=V(E(itet,3),:)';
    x4=V(E(itet,3),:)';

    if any(any(isnan([x1 x2 x3 x4]))) % if any nan, iso index=nan
        tetIsoIndex(itet)=NaN;
    else
        xa=(x1+x2+x3+x4)/3;
        X=[x1-xa,x2-xa,x3-xa,x4-xa]; %cluster position model
        K=X*X'/4;
        Keig=real(eig(K)); %eigenvalues of X
        tetIsoIndex(itet)=3*Keig(1)/(Keig(1)+Keig(2)+Keig(3)); % isotropy index
    end
end

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
classdef task_pf_legacy < mp.task_pf & mp.task_shared_legacy
%MP.TASK_PF_LEGACY  MATPOWER task for legacy power flow (PF).
%   MP.TASK_PF_LEGACY provides implementation for power flow problem.
%
%   Properties
%       ?
%
%   Methods
%       ?
%
%   See also MP.TASK_PF

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

    properties
    end

    methods
        %%-----  task methods  -----
        function [d, mpopt] = run_pre(obj, d, mpopt)
            [d, mpopt] = obj.run_pre_legacy(d, mpopt);
            [d, mpopt] = run_pre@mp.task_pf(obj, d, mpopt);
        end

        function obj = run_post(obj, mm, nm, dm, mpopt);
            if obj.nm.np ~= 0
                obj.dm.source = obj.dmc.export(obj.dm, obj.dm.source);
            end
        end

        %%-----  other methods  -----
        function [results, success] = legacy_post_run(obj, mpopt)
            success = obj.success;
            results = obj.dm.source;
            results.task = obj;
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


function [prices, alphas] = Heston1993KahlJaeckelLordRev3(PC, S,K,T,t,r,q,v0,theta,rho,kappa,sigma, alphas)
% Heston pricing function based on the implementation suggested by
% Roger Lord and Chrisitan Kahl in "Optimal Fourier inversion in
% semi-analytical option pricing"
%
%
% Input: (PC till q can be vectorized)
%       PC: 1 for Calls, 2 for Puts
%       S: Spot
%       K: Strike
%       T: Maturity
%       t: start date
%       r: interest rate
%       q: dividend
%       v0: initial variance
%       theta: long run mean variance
%       kappa: mean reversion speed of  volatility
%       sigma: volatility of volatility
%       rho: correlation between returns volatility
%       alpha: alpha can be a vector supplied by the user, otherwise the
%       function attempts to find a payoff-dependent optimal alpha
%
%   Output: Price for each option, optionally generated alphas
%
%   Usage: Heston1993KahlJaeckelLordRev3(1, 100, 100, 20,0, 0.05, 0.0,
%           0.00003, 0.00003,-0.3, 0.5, 0.0008)
%
%   Author: Jonathan Frei, 2015
%

    % force column vector
    PC=PC(:);
    S=S(:);
    K=K(:);
    T=T(:);
    t=t(:);
    r=r(:);
    q=q(:);

    nos = numel(S);
    prices=NaN(nos,1);
    tau=T-t;
    mu=(r-q);
    F = S.*exp(mu.*tau);

    if(~exist('alphas','var'))
        alphas = NaN(numel(S),1);
    elseif(numel(alphas)==1)
        alphas = repmat(alphas,numel(S),1);
    end

    alpha0=0.75;

    for(ind=1:nos)
        if(isnan(alphas(ind)))
            try
                % using fzero here instead of fminsearch
               alphas(ind) = fzero( @(a) psi(a,K(ind), F(ind), kappa, theta, rho, sigma, tau(ind), v0), alpha0);
            catch
               alphas(ind) = alpha0;
            end
        end
        prices(ind) =  Ralpha(F(ind), K(ind), alphas(ind))+1/pi*integral(@(x) phi(x, K(ind), alphas(ind), F(ind), kappa, theta, rho, sigma, tau(ind), v0) , 0, Inf);
        if (PC(ind)==2)
            prices(ind) = prices(ind) + K(ind)*exp(-r(ind)*tau(ind))-S(ind)*exp(-q(ind)*tau(ind));
        end
    end

end

function p = psi(alpha, K, F, kappa, theta, rho, sigma, tau, v0)
    k = log(K);
    p = -alpha*k+0.5*log(phi(-(alpha+1)*1i, K, alpha, F, kappa, theta, rho, sigma, tau, v0)^2);
end

function r = Ralpha(F, K, alpha)
    r = F*(alpha<=0)-K*(alpha<=-1)-0.5*(F*(alpha==0)-K*(alpha==-1));
end

function y = phi(v, K, alpha, F, kappa, theta, rho, sigma, tau, v0)
    k = log(K);
    y = real(exp(-1i*(v-1i*alpha)*k).*( cf(v-1i*(alpha+1), F, kappa, theta, rho, sigma, tau, v0)./(-(v-1i*(alpha+1)).*(v-1i*alpha))));
end

function c = cf(u, F, kappa, theta, rho, sigma, tau, v0)
    f = log(F);
    c = exp(1i*u*f+ A(u, kappa, theta, rho, sigma, tau)+Bv(u, rho, sigma, kappa, tau)*v0);
end

function b = Bv(u, rho, sigma, kappa, tau)
    b = ((beta(u,rho,sigma,kappa)-D(u, rho, sigma, kappa)).*(1-exp(-D(u, rho, sigma, kappa)*tau)))./(sigma.^2*(1-G(u, rho, sigma, kappa).*exp(-D(u, rho, sigma, kappa)*tau)));
end

function a = A(u, kappa, theta, rho, sigma, tau)
    a = (kappa*theta*((beta(u,rho,sigma,kappa)-D(u, rho, sigma, kappa))*tau-2*log(phi2(u, rho, sigma, kappa, tau))))/sigma.^2;
end

function p = phi2(u, rho, sigma, kappa, tau)
    p = (G(u, rho, sigma, kappa).*exp(-D(u, rho, sigma, kappa)*tau)-1)./(G(u, rho, sigma, kappa)-1);
end

function g = G(u, rho, sigma, kappa)
    g = (beta(u,rho,sigma,kappa)-D(u, rho, sigma, kappa))./(beta(u,rho,sigma,kappa)+D(u, rho, sigma, kappa));
end

function d = D(u, rho, sigma, kappa)
    d = sqrt(beta(u,rho,sigma,kappa).^2-4*alphahat(u)*gamma(sigma));
end

function a = alphahat(u)
    a = -0.5*u.*(1i+u);
end

function b = beta(u,rho,sigma,kappa)
    b = kappa-rho*sigma*u*1i;
end

function y = gamma(sigma)
    y = 0.5*sigma.^2;
end


text:
function [newW] = update(eta, Y, W)
%UPDATE gradient value W = W + delW

newW = zeros(size(W));
newW = W + gradient(eta, Y, W);

end

author: vsubhashini

repo: ica

save_path: github-repos/MATLAB/vsubhashini-ica


text:
function fighandle=vf(varargin)
% fighandle=vf(figs)      View one or more Figures.
% This simple function brings the group of figures specified by "figs" to
% the front. No tiling or rearrangment of the figures is done.
% If "figs" is empty then all of the figures are brought to the front.
%   "fighandle" has the figure handles of the figures brought forward.
%   Examples:
%       vf(a) where "a" is a vector of figure handles.
%       a=vf, vf a, vf [1,6,8], a=vf(1,6,8), vf 1 a 6
% Note that the command line format is supported.

%   Copyright 2006 Mirtech, Inc.
%   created 08/22/2006  by Mirko Hrovat on Matlab Ver. 7.2
%   Mirtech, Inc.       email: mhrovat@email.com

figs=[];
if nargin>0
    for n=1:length(varargin),
        arg=varargin{n};
        if ischar(arg),
            tmp=str2num(arg);   %#ok
            if isempty(tmp),    % since it is empty then it may be a variable
                tmp = evalin('base',arg);
                if isempty(tmp),
                    error ('  Syntax error in argument!'),
                end
            end  % isempty(tmp)
        else
            tmp=arg;
        end  % ischar(arg)
        figs=[figs,tmp];
    end  % for
end

% set this to "on" so that all figures with "HandleVisibility" set to off can be seen.
set(0,'ShowHiddenHandles','on')
% Below are two alternative methods to find all of the figures.
% Second method should be more foolproof.
% allhands=get(0,'Children');
allhands=findobj('Type','figure');
set(0,'ShowHiddenHandles','off')
if isempty(figs),
    figs=sort(allhands)';
else
    [tmp,indx]=intersect(figs,allhands);    %#ok tmp not used, view figures that exist
    figs=figs(sort(indx));                  % order the figures as specified
end
if isempty(figs),
    warning('MATLAB:EmptyValue','There are no figures to select.')
else
    for n=numel(figs):-1:1,
        figure(figs(n))                         % bring figure to front.
    end
end

if nargout==1,
    fighandle=figs;
end

author: Sable

repo: mcbench-benchmarks


text:
function [ r, info ] = r8pp_fa ( n, a )

%*****************************************************************************80
%
%% R8PP_FA factors a R8PP matrix.
%
%  Discussion:
%
%    The R8PP storage format is appropriate for a symmetric positive
%    definite matrix.  Only the upper triangle of the matrix is stored,
%    by successive partial columns, in an array of length (N*(N+1))/2,
%    which contains (A11,A12,A22,A13,A23,A33,A14,...,ANN)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2004
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A((N*(N+1))/2), a R8PP matrix.
%
%    Output, real R((N*(N+1))/2), an upper triangular matrix R, stored
%    in packed form, so that A = R'*R.
%
%    Output, integer INFO, error flag.
%    0, for normal return.
%    K, if the leading minor of order K is not positive definite.
%
  info = 0;

  r(1:(n*(n+1))/2) = a(1:(n*(n+1))/2);

  jj = 0;

  for j = 1 : n

    s = 0.0E+00;
    kj = jj;
    kk = 0;

    for k = 1 : j-1

      kj = kj + 1;
      t = r(kj);
      for i = 1 : k-1
        t = t - r(kk+i) * r(jj+i);
      end
      kk = kk + k;
      t = t / r(kk);
      r(kj) = t;
      s = s + t * t;

    end

    jj = jj + j;
    s = r(jj) - s;

    if ( s <= 0.0E+00 )
      info = j;
      return;
    end

    r(jj) = sqrt ( s );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%==========================================================================
% Walks through a tree structure data. Whenever it finds a structure, which
% have field named 'import' it assumes that in that field is a cell array
% and merges all structures found in that array. Parameter verb is used for
% debugging purposes.
%==========================================================================
function result = mergeimports(data, verb)
    if ~exist('verb','var')
        verb = 0;
    end;
    result = recurse(data, 0, [], verb);
end

%--------------------------------------------------------------------------
% Recursion dispatcher, calls appropriate method for cell/structure or
% displays data if the parameter data is not of mentioned type.
%   addit ... for possible future use, now unused
%   verb  ... for debugging
%
function result = recurse(data, level, addit, verb)
    indent = repmat(' | ',1,level); % for debugging
    if iscell(data)
        result = iter_cell(data, level, addit, verb);
    elseif isstruct(data)
        result = iter_struct(data, level, addit, verb);
    else
        if any(verb == 1) % for debugging
            fprintf([indent,'Some data: ']);
            disp(data);
        end;
        result = data;
    end;
end

%--------------------------------------------------------------------------
% Walks through a cell array and calls recurse on every field.
%   data        ... Assumed to be a cell. Data to be walked.
%   level       ... Level in the tree, root has zero.
%   addit, verb ... for debugging
%
function result = iter_cell(data, level, addit, verb)
    indent = repmat(' | ',1,level); % for debugging
    result = {};
    if any(verb == 1); fprintf([indent,'cell {\n']); end; % for debugging
    for i = 1:length(data)
        itemcontent = recurse(data{i}, level + 1, addit, verb);
        result{end + 1} = itemcontent;
    end;
    if any(verb == 1); fprintf([indent,'} cell\n']); end; % for debugging
end

%--------------------------------------------------------------------------
% Walks through a struct and calls recurse on every field. If there is a
% field called 'import' it calls process_import_field on its content. Then
% merges processed import with the rest of the structure. Meaning of all
% parameters is similar to those of iter_cell.
%
function result = iter_struct(data, level, addit, verb)
    indent = repmat(' | ',1,level); % for debugging
    result = struct();
    collected_imports = {};
    if any(verb == 1); fprintf([indent,'struct {\n']); end; % for debugging
    for i = fields(data)'
        fld = char(i);
        if any(verb == 1); fprintf([indent,' +-field ',fld,':\n']); end; % for debugging
        result.(fld) = recurse(data.(fld), level + 1, addit, verb);
        % Tree back-pass - all potential underlying imports were processed,
        % so process import here, if needed.
        if isequal(fld, 'import')
            processed_import = process_import_field(result.(fld));
            result = rmfield(result, 'import');
            if isstruct(processed_import)
                collected_imports{end+1} = processed_import;
                % It is maybe useless to collect imports to the array since
                % there can be only one field named 'import' per structure.
                % collected_imports is proposed to be changed to a simple
                % variable.
            else
                % One of imports was not a struct. In following versions it
                % probably won't be an error and will merge to a cell with
                % structs.
                disp(processed_import);
                error('Expected struct, otherwise it cannot be merged with the rest.');
            end;
        end;
    end;
    for i = 1:length(collected_imports)
        result = merge_struct(result, collected_imports{i}, {}, 'deep');
    end;
    if any(verb == 1); fprintf([indent,'} struct\n']); end; % for debugging
end

%--------------------------------------------------------------------------
% Walks through the data parameter, which is assumed to be a cell. Merges
% all structures in that cell and returns them as a struct or possibly as a
% cell of merged struct and unmeregeable data.
%
function result = process_import_field(data)
    if iscell(data)
        merged_structs = struct();
        collected_nonstruct = {};
        for i = 1:length(data)
            if isstruct(data{i})
                merged_structs = merge_struct(merged_structs, data{i}, {}, 'deep');
            else
                collected_nonstruct{end+1} = data{i};
            end;
        end;
        if isempty(collected_nonstruct)
            result = merged_structs;
        elseif isempty(merged_structs)
            result = collected_nonstruct;
        else
            result = {merged_structs; collected_nonstruct};
        end;
    else
        % For clarity and simplicity, the whole transformation is done so
        % that every import field in a struct is cell array even there is
        % only one object to be imported.
        error('BUG: import field should always contain a cell.');
    end;
end

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

author: ziyuw

repo: rembo

save_path: github-repos/MATLAB/ziyuw-rembo


text:
function Q = up_expert_findQ(b, data, del0, delta)
% This program finds Q for KV02
%
% function Q = up_expert_findQ(b, data, del0, delta, N)
%
% Q: A modification of the random walk wealth distribution P
%
% b: portfolio
% data: market sequence
% del0, delta: parameters
%
% Example: Q = up_expert_findQ(b, data, del0, delta, N);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[~, N] = size(data);  % Number of stocks,Time period.

P = prod(data(:, :)*b);

Q = P*min(1, exp((b(N)-(2*del0))/(N*delta)));

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function [z,w] = spm_DEM_z(M,N)
% creates hierarchical innovations for generating data
% FORMAT [z w] = spm_DEM_z(M,N)
% M    - model structure
% N    - length of data sequence
%
% z{i} - innovations for level i (N.B. z{end} corresponds to causes)
% w{i} - innovations for level i (state noise)
%
% If there is no fixed or hyper parameterized precision, then unit noise is
% created. It is assumed that this will be later modulated by state
% dependent terms, specified by M.ph and M.pg in spm_DEM_int
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_DEM_z.m 7540 2019-03-11 10:44:51Z karl $

% temporal convolution matrix (with unit variance)
%--------------------------------------------------------------------------
s  = M(1).E.s + exp(-16);
dt = M(1).E.dt;
t  = ((1:N) - 1)*dt;
K  = toeplitz(exp(-t.^2/(2*s^2)));
K  = K*diag(1./sqrt(diag(K*K')));

% create innovations z{i} and w{i}
%--------------------------------------------------------------------------
for i = 1:length(M)

    % precision of causes
    %======================================================================
    P     = M(i).V;

    % plus prior expectations
    %----------------------------------------------------------------------
    try
        for j = 1:length(M(i).Q)
            P = P + M(i).Q{j}*exp(M(i).hE(j));
        end
    end

    % create causes: assume i.i.d. if precision is zero
    %----------------------------------------------------------------------
    if norm(P,1) == 0;
        z{i}  = randn(M(i).l,N)*K;
    elseif norm(P,1) >= exp(16)
        z{i}  = sparse(M(i).l,N);
    else
        z{i}  = spm_sqrtm(inv(P))*randn(M(i).l,N)*K;
    end

    % precision of states
    %======================================================================
    P     = M(i).W;

    % plus prior expectations
    %----------------------------------------------------------------------
    try
        for j = 1:length(M(i).R)
            P = P + M(i).R{j}*exp(M(i).gE(j));
        end
    end

    % create states: assume i.i.d. if precision (P) is zero
    %----------------------------------------------------------------------
    if ~isempty(P)
        if norm(P,1) == 0;
            w{i} = randn(M(i).n,N)*K*dt;
        elseif norm(P,1) >= exp(16)
            w{i} = sparse(M(i).n,N);
        else
            w{i} = spm_sqrtm(inv(P))*randn(M(i).n,N)*K*dt;
        end
    else
        w{i} = sparse(0,0);
    end

end

author: spm


text:
function x = proxTVSquare(x,lambda,alpha,varargin)

[n,m] = size(x);
lambda = alpha * lambda; % only the product matters

% apply proxTV horizontally
mMax = 2*floor((m+1)/2)-1; % round down to an odd number
[x(:,2:2:mMax),x(:,3:2:mMax)] = proxTV(x(:,2:2:mMax),x(:,3:2:mMax),lambda,varargin{:});

mMax = 2*floor(m/2); % round down to an even number
[x(:,1:2:mMax),x(:,2:2:mMax)] = proxTV(x(:,1:2:mMax),x(:,2:2:mMax),lambda,varargin{:});

% apply proxTV vertically
nMax = 2*floor((n+1)/2)-1; % round down to an odd number
[x(2:2:nMax,:),x(3:2:nMax,:)] = proxTV(x(2:2:nMax,:),x(3:2:nMax,:),lambda,varargin{:});

nMax = 2*floor(n/2); % round down to an even number
[x(1:2:nMax,:),x(2:2:nMax,:)] = proxTV(x(1:2:nMax,:),x(2:2:nMax,:),lambda,varargin{:});

author: mtex-toolbox

repo: mtex

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


text:
function [saveLogic,savePath]=Qsave3DDICPPresults(structPath)
%% function for asking the user for saving options for the 3D-DIC results in STEP3.
% Including overwriting options
%
% INPUTS:
% * structPaths: 1-by-Npairs cell array, containing in each
%   cell the path to the DIC2DpairResults_C_#1_C_#2 structure file
%
% OUTPUTS:
% * saveLogic: logic true=save, false=not save
% * savePath: a path were to save the results. [] if saveLogic=false.

%%
saveButton = questdlg('Save 3D-DIC Post-Processing results?', 'Save?', 'Yes', 'No', 'Yes');

switch saveButton
    case 'Yes'
        saveLogic=true(1);
        continueLog=false(1);

        while ~continueLog
            % Path where to save the camera parameters
            savePath = uigetdir(fileparts(structPath),'Select a folder for saving the results');

            % check if DIC3DpairResults in this folder already exist
            fileExist=false;
            [~,name,~]=fileparts(structPath);
            filenames=[savePath '\DIC3DPPresults.mat'];
            if exist(filenames,'file')
                fileExist=true;
            end

            if fileExist
                promptMessage=cell(2,1);
                promptMessage{1}='The following output file already exists: ';
                promptMessage{5}='Do you want to overwrite?';
                [~,fileName,ext]=fileparts(filenames);
                promptMessage{3}=[fileName ext];

                titleBarCaption = 'Overwrite?';
                buttonText = questdlg(promptMessage, titleBarCaption, 'Yes', 'No', 'Yes');
                switch buttonText
                    case 'Yes'
                        % User wants to overwrite. % Set flag to continue
                        continueLog=true(1);
                    case 'No'
                        % User does not want to overwrite. % Set flag to not do the write.
                        waitfor(warndlg('Choose another folder to save the results','OK'));
                        % let user choose a different folder to save the file
                end
            else
                % if file doesn't esixt, cotinue.
                continueLog=true(1);
            end
        end

    case 'No'
        saveLogic=false(1);
        savePath=[];
end

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>


text:
function model = gpReversibleDynamicsExpandParam(model, params)

% GPREVERSIBLEDYNAMICSEXPANDPARAM Place the parameters vector into the model for GP dynamics.

% FGPLVM

model = gpDynamicsExpandParam(model, params);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat