text:
classdef cheboppref < chebpref
%CHEBOPPREF   Class for managing preferences for the Chebfun ODE suite.
%   CHEBOPPREF is a class for managing CHEBOP construction-time and solver
%   preferences, such as what solver is used for linear problem, the error or
%   residual tolerance for nonlinear problems, whether damped Newton iteration
%   should be performed for nonlinear problems, and how much information is to
%   be printed to the console during while the solver is active.
%
% Available Preferences:
%
%   bvpTol                      - Error tolerance for boundary value problems
%     [5e-13]
%
%     This is the relative tolerance used to test convergence during the
%     adaptive solution of linear boundary value problems. For nonlinear
%     boundary value problems the Newton convergence tolerance is set to
%     200*bvpTol.
%
%   domain                     - Construction domain.
%     [-1, 1]
%
%     This sets the default domain that will be used for CHEBOP construction if
%     no domain argument is explicitly passed to the constructor.
%
%   discretization             - Discretization of linear problems
%     ['values']
%     'coeffs'
%     @chebcolloc1
%     @chebcolloc2
%     @ultraS
%     @trigcolloc
%     @trigspec
%
%     This options determines whether linear operators are discretized using
%     rectangular collocation methods or the ultraspherical method. Please
%     observe that
%         * 'values' and 'coeffs' are convenient ways of specifying the
%           @chebcolloc2 and @ultraS options respectively (when the boundary
%           conditions are not periodic), and @trigcolloc and @trigspec
%           respectively (when the boundary conditions are periodic).
%         * The @trigcolloc/@trigspec options are only supported for problems
%           that are specified to have periodic boundary conditions.
%         * Specifying the @chebcolloc1 option causes the CHEBFUN solution
%           returned to be based on the @chebtech1 tech. The @chebtech2/@ultraS
%           option causes the CHEBFUN solution returned to be based on the
%           @chebtech2 tech. The @trigcolloc/@trigspec option causes the
%           CHEBFUN solution to be periodic, based on the @trigtech tech.
%
%   damping                     - Should Newton's method be damped?
%     [true]
%     false
%
%     If true, damped Newton iteration in function space is performed for
%     nonlinear problems. If false, undamped Newton iteration is performed, that
%     is, the solver will always take full Newton steps.
%
%   display                     - How much information is to be printed
%     'final'
%     ['iter']
%      'off'
%
%     If 'final', information is only printed after the solver of BVPs has
%     finished. If 'iter', information is printed at every Newton step. If
%     'off', no information is printed.
%
%   happinessCheck              - Routine for checking that solution converged
%     [@standardCheck]
%     @basicCheck
%     @plateauCheck
%     @classicCheck
%     @looseCheck
%     @strictCheck
%     @happinessCheck
%     @linopV4Check
%
%     This options determines which routine is used to determine that the
%     approximate solution has converged. Any of the above options may be
%     used, as well as any user defined function handle that conforms to
%     the happinessCheck standards.
%
%   ivpAbsTol                    - Absolute tolerance for the ivpSolver
%     [1e5*eps]
%
%     This options specifies the option for the absolute tolerance passed as an
%     option to the built-in MATLAB ODE solver when solving IVPs.
%
%   ivpRelTol                    - Relavtive tolerance for the ivpSolver
%     [100*eps]
%
%     This options specifies the option for the relative tolerance passed as an
%     option to the built-in MATLAB ODE solver when solving IVPs.
%
%   ivpRestartSolver             - Restart IVP solvers at breakpoints
%     false
%     [true]
%
%     This option specifies whether the MATLAB built in solvers should be
%     restarted at breakpoints. That is, whether each subinterval of a piecewise
%     problem will get integrated separately. This can be very useful for e.g.
%     short forcing pulses, which otherwise might get overlooked.
%
%   ivpSolver                  - Solver for IVPs
%     ['ode113']
%     'ode15s'
%     'ode45'
%     'values'
%     'coeffs'
%
%     This options determines which of the MATLAB built-in IVP solvers is used
%     for solving IVPs posed with the CHEBOP class. Any option of
%     CHEBOPPREF.discretization (see above) is allowed, which causes IVPs to be
%     solved globally via spectral methods, rather than reformulating them as
%     first-order problems and then solved via time-stepping method.
%
%   lambdaMin                   - Minimum allowed step-size for Newton's method
%     [1e-6]
%
%     The value of lambdaMin determines the minimum allowed step-size that the
%     damped Newton iteration is allowed to take.
%
%   maxDimension
%     [4096]
%
%     The maximum number of gridpoints/coefficients used as linear operators are
%     discretized at finer and finer grids to resolve the solution. The
%     intermediate values for the discretization between cheboppref.minDimension
%     and cheboppref.maxDimension depend on the discretization used for the
%     operator.
%
%   maxIter                     - Maximum number of Newton steps
%     [25]
%
%   The maximum number of steps that the (damped) Newton iteration is allowed to
%   take, before it is considered to be non-convergent.
%
%   minDimension
%     [32]
%
%     The minimum number of gridpoints/coefficients used as linear operators are
%     discretized at finer and finer grids to resolve the solution. The
%     intermediate values for the discretization between cheboppref.minDimension
%     and cheboppref.maxDimension depend on the discretization used for the
%     operator.
%
%   plotting                    - Plotting of intermediate Newton steps
%     DELAY
%     'on'
%     ['off']
%     'pause'
%
%   If plotting = 'on', the current iterate in the Newton solution is plotted at
%   every step, as well as the current Newton correction. If plotting = DELAY,
%   where DELAY has a numerical value, the iteration is paused and the plots are
%   shown for the time DELAY seconds. If plotting = 'pause', the iteration is
%   paused and the plots are shown until the user presses a button. If plotting
%   = 'off', no plots are shown during the Newton iteration.
%
%   vectorize                   - Automatic vectorization of anon. functions
%     [true]
%     false
%
%   Determines whether the CHEBOP class should try to automatically try to
%   vectorize anonymous functions used for describing the differential equation
%   and boundary condition(s).
%
%
% The default values for any of these preferences may be globally overridden
% using CHEBOPPREF.SETDEFAULTS(); see the documentation for that function for
% further details.
%
% Constructor inputs:
%   P = CHEBOPPREF() creates a CHEBOPPREF object with the default values of the
%   preferences.  For a list of all available preferences, see above.
%
%   P = CHEBOPPREF(Q), where Q is a MATLAB structure uses the field/value pairs
%   in Q to set the properties of P.  If a field of Q has a name which matches
%   a property of P, the value of that property of P is set to the value
%   associated to that field in Q.  If a field of Q has a name that does not
%   correspond to a known preference, then an error is thrown.
%
%   P = CHEBOPPREF(Q), where Q is a CHEBOPPREF, sets P to be a copy of Q.
%
% See also CHEBFUNPREF.

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

% TODO:  Further documentation of CHEBOPPREF preferences.

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS CONSTRUCTOR:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

        function outPref = cheboppref(inPref, varargin)
            if ( (nargin == 1) && isa(inPref, 'cheboppref') )
                outPref = inPref;
                return
            elseif ( nargin < 1 )
                inPref = struct();
            elseif ( ischar(inPref) )
                if ( nargin == 1 )
                    error('CHEBFUN:CHEBOPPREF:cheboppref:deprecated', ...
                        ['cheboppref() no longer supports queries of ', ...
                         'the form cheboppref(''prop'').\n', ...
                         'Please use cheboppref().prop.']);
                else
                    error('CHEBFUN:CHEBOPPREF:cheboppref:deprecated', ...
                        ['cheboppref() no longer supports assignment ', ...
                         'via cheboppref(''prop'', val).\n', ...
                         'Please use cheboppref.setDefaults(''prop'', val).']);
                end
            elseif ( nargin > 1 )
                error('CHEBFUN:CHEBOPPREF:cheboppref:inputs', ...
                    'Too many input arguments.')
            end

            % Initialize default preference values.
            outPref.prefList = cheboppref.manageDefaultPrefs('get');

            % Copy fields from q, merging incomplete substructures.
            for field = fieldnames(inPref).'
                field1 = field{1};
                if ( isfield(outPref.prefList, field1) )
                    if ( isstruct(outPref.prefList.(field1)) )
                        outPref.prefList.(field1) = ...
                            chebpref.mergePrefStructs(...
                                outPref.prefList.(field1), ...
                                inPref.(field1));
                    else
                        outPref.prefList.(field1) = inPref.(field1);
                    end
                else
                    error('CHEBFUN:CHEBOPPREF:cheboppref:badPref', ...
                        'Unrecognized preference name.');
                end
            end
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% CLASS METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = false )

       function display(pref)
       %DISPLAY   Display a CHEBOPPREF object.
       %   DISPLAY(PREF) prints out a list of the preferences stored in the
       %   CHEBOPPREF object PREF.

            % Compute the screen column in which pref values start.
            valueCol = 24; % length('    blowup:   ');

            % A subfunction to pad strings for formatting.
            function s = padString(s)
            %PADSTRING   Add whitespace to string for formatting.
                s = [s repmat(' ', 1, valueCol - length(s))];
            end

            % Print values of "known" preferences.
            prefList = pref.prefList;

            fprintf('cheboppref object with the following preferences:\n');
            fprintf([padString('    bvpTol:') '%g\n'], ...
                prefList.bvpTol);
            fprintf([padString('    domain:') '[%g, %g]\n'], ...
                prefList.domain(1), prefList.domain(end));
            if ( isa(prefList.discretization,'function_handle') )
                fprintf([padString('    discretization:') '%s\n'], ...
                    func2str(prefList.discretization));
            elseif ( isa(prefList.discretization,'char') )
                fprintf([padString('    discretization:') '%s\n'], ...
                    prefList.discretization);
            end

            fprintf([padString('    damping:') '%d\n'], ...
                prefList.damping);
            fprintf([padString('    display:') '%s\n'], ...
                prefList.display);
            fprintf([padString('    happinessCheck:') '%s\n'], ...
                func2str(prefList.happinessCheck));
            fprintf([padString('    ivpAbsTol:') '%g\n'], ...
                prefList.ivpAbsTol);
            fprintf([padString('    ivpRelTol:') '%g\n'], ...
                prefList.ivpRelTol);
            fprintf([padString('    ivpRestartSolver:') '%d\n'], ...
                prefList.ivpRestartSolver);
            fprintf([padString('    ivpSolver:') '%s\n'], ...
                func2str(prefList.ivpSolver));
            fprintf([padString('    lambdaMin:') '%g\n'], ...
                prefList.lambdaMin);
            fprintf([padString('    maxDimension:') '%d\n'], ...
                prefList.maxDimension);
            fprintf([padString('    maxIter:') '%d\n'], ...
                prefList.maxIter);
            fprintf([padString('    minDimension:') '%d\n'], ...
                prefList.minDimension);
            fprintf([padString('    plotting:') '%s\n'], ...
                prefList.plotting);
            fprintf([padString('    vectorize:') '%i\n'], ...
                prefList.vectorize);
       end

        function pref = subsasgn(pref, ind, val)
        %SUBSASGN   Subscripted assignment for CHEBOPPREF.
        %   P.PROP = VAL, where P is a CHEBOPPREF object, assigns the value
        %   VAL to the CHEBOPPREF property PROP stored in P.  If PROP is not a
        %   CHEBOPPREF property, an error will be thrown.
        %
        %   CHEBOPPREF does not support any other subscripted assignment types,
        %   including '()' and '{}'.

            % Support user-friendlier syntax for specifying discretization
            % choice:
            val = cheboppref.parseDiscretization(val);

            % Support user-friendlier syntax for specifying IVP solver choice:
            val = cheboppref.parseIVPsolver(val);

            % Call the superclass method.
            pref = subsasgn@chebpref(pref, ind, val);
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% STATIC METHODS:
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Access = public, Static = true )

        function pref = getFactoryDefaults()
        %GETFACTORYDEFAULTS   Get factory default preferences.
        %   PREF = CHEBOPPREF.GETFACTORYDEFAULTS() returns a CHEBOPPREF
        %   object with the preferences set to their factory defaults,
        %   irrespective of the currently defined values of the default
        %   preferences.  This function is useful if the user wishes to
        %   solve ODEs with CHEBOP using the factory defaults when other
        %   user-set defaults are currently in force.
        %
        % See also SETDEFAULTS.

            fd = cheboppref.factoryDefaultPrefs();
            pref = cheboppref(fd);
        end

        function setDefaults(varargin)
        %SETDEFAULTS   Set default preferences.
        %   CHEBOPPREF.SETDEFAULTS(PREF1, VAL1, PREF2, VAL2, ...) sets the
        %   default values for the preferences whose names are stored in the
        %   strings PREF1, PREF2, ..., etc. to VAL1, VAL2, ..., etc.  All
        %   subsequently constructed CHEBOPPREF objects will use these values
        %   as the defaults.
        %
        %   CHEBOPPREF.SETDEFAULTS(PREF) sets the default values to the
        %   preferences stored in the CHEBOPPREF object PREF.  PREF can also
        %   be a MATLAB structure, in which case it is converted to a
        %   CHEBOPPREF as described in the documentation for the CHEBOPPREF
        %   constructor first.
        %
        %   CHEBOPPREF.SETDEFAULTS('factory') resets the default preferences to
        %   their factory values.
        %
        % See also GETFACTORYDEFAULTS.

            % The reason we don't just use manageDefaults as the second
            % argument to chebpref.setDefaults and wrap it in an additional
            % anonymous function instead is to get around what seems to be a
            % bug in MATLAB.  See commit messages for more information.
            manageDefaults = @cheboppref.manageDefaultPrefs;
            chebpref.setDefaults(@(inPref) cheboppref(inPref), ...
                @(varargin) manageDefaults(varargin{:}), varargin{:});
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% PRIVATE STATIC METHODS
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods ( Static = true, Access = private )

        function varargout = manageDefaultPrefs(varargin)
        %MANAGEDEFAULTPREFS   Private method for handling default preferences.
        %   CHEBOPPREF.MANAGEDEFAULTPREFS('get') returns a structure suitable
        %   for storing in the prefList property of a CHEBOPPREF with all of
        %   the currently stored default preferences suitable for initializing
        %   a CHEBOPPREF object.
        %
        %   CHEBOPPREF.MANAGEDEFAULTPREFS('set-factory') restores the default
        %   preferences to their "factory" values.
        %
        %   CHEBOPPREF.MANAGEDEFAULTPREFS('set', PREFLIST) sets the default
        %   values to those stored in the structure PREFLIST.  PREFLIST should
        %   be a structure suitable for use as a CHEBOPPREF prefList.
        %
        %   CHEBOPPREF.MANAGEDEFAULTPREFS('set', PREF1, VAL1, PREF2, VAL2, ...)
        %   sets the default values for PREF1, PREF2, ..., etc. to VAL1, VAL2,
        %   ..., etc.

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % DEVELOPER NOTE:
        %  - MATLAB has no equivalent to what might be called a "static" class
        %    variable in other languages, so a persistent variable is the best
        %    we can do for providing this feature.  Persistent variables are
        %    local to a specific function, so we can have only a single
        %    function for managing it.  As a result, this function has a mildly
        %    awkward syntax and so is not user-facing.
        %  - More importantly, this function is also not user-facing because
        %    its inputs and outputs depend on the internal representation of a
        %    CHEBOPPREF as a MATLAB structure, and that's not something with
        %    which anyone outside of CHEBOPPREF should be concerned.
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            persistent defaultPrefs;

            if ( isempty(defaultPrefs) )
                defaultPrefs = cheboppref.factoryDefaultPrefs();
            end

            if ( strcmp(varargin{1}, 'get') )
                varargout{1} = defaultPrefs;
            elseif ( strcmp(varargin{1}, 'set-factory') )
                defaultPrefs = cheboppref.factoryDefaultPrefs();
            elseif ( strcmp(varargin{1}, 'set') )
                    varargin(1) = [];
                if ( isstruct(varargin{1}) )
                    defaultPrefs = varargin{1};
                else
                    while ( ~isempty(varargin) )
                        prefName = varargin{1};
                        prefValue = varargin{2};

                        % Support user-friendlier syntax for specifying
                        % discretization choice:
                        prefValue = cheboppref.parseDiscretization(prefValue);
                        prefValue = cheboppref.parseHappinessCheck(prefValue);
                        prefValue = cheboppref.parseIVPsolver(prefValue);
                        if ( isfield(defaultPrefs, prefName) )
                            defaultPrefs.(prefName) = prefValue;
                        else
                            error('CHEBFUN:CHEBOPPREF:cheboppref:badPref', ...
                                'Unrecognized preference name.');
                        end
                        varargin(1:2) = [];
                    end
                end
            end
        end

        function factoryPrefs = factoryDefaultPrefs()
        %FACTORYDEFAULTPREFS   Get structure of factory default preferences.
        %   S = CHEBOPPREF.FACTORYDEFAULTPREFS() returns a structure suitable
        %   for storing in the prefList property of a CHEBOPPREF object that
        %   contains all of the "factory default" values of the CHEBOP
        %   preferences.

            factoryPrefs.bvpTol = 5e-13;
            factoryPrefs.domain = [-1 1];
            factoryPrefs.discretization = 'values';
            factoryPrefs.scale = NaN;
            factoryPrefs.damping = 1;
            factoryPrefs.display = 'off';
            factoryPrefs.happinessCheck = @standardCheck;
            factoryPrefs.ivpAbsTol = 1e5*eps;
            factoryPrefs.ivpRelTol = 100*eps;
            factoryPrefs.ivpRestartSolver = true;
            factoryPrefs.ivpSolver = @chebfun.ode113;
            factoryPrefs.lambdaMin = 1e-6;
            factoryPrefs.maxDimension = 4096;
            factoryPrefs.maxIter = 25;
            factoryPrefs.minDimension = 32;
            factoryPrefs.plotting = 'off';
            factoryPrefs.vectorize = true;
        end

        function val = parseDiscretization(val)
        %PARSEDISCRETIZATION    Allow different syntax for specifying
        %                       discretization.

            % We want to allow user-friendly syntax for specifying the
            % discretization (#433). So check whether we have some of the
            % strings we want to allow, and convert them to the correct function
            % handle:
            if ( any(strcmpi(val, {'ultraspherical', 'ultraS'})) )
                warning('CHEBOPPREF:PARSEDISCRETIZATION', ...
                    ['''ULTRAS''/''ULTRASPHERICAL'' is deprecated. \n' ...
                    'Please use ''COEFFS''/@ultraS.']);
                val = @ultraS;

            elseif ( any(strcmpi(val, {'chebcolloc2', 'collocation', 'colloc2'})) )
                warning('CHEBOPPREF:PARSEDISCRETIZATION', ...
                    ['''COLLOCATION''/''COLLOC2''/''CHEBCOLLOC2'' is deprecated. \n' ...
                    'Please use ''VALUES''/@chebcolloc2.']);
                val = @chebcolloc2;

            elseif ( any(strcmpi(val, {'chebcolloc1', 'colloc1'})) )
                warning('CHEBOPPREF:PARSEDISCRETIZATION', ...
                    ['''COLLOC1''/''CHEBCOLLOC1'' is deprecated. \n' ...
                    'Please use ''VALUES''/@chebcolloc2.']);
                val = @chebcolloc1;

            elseif ( any(strcmpi(val, {'trigcolloc', 'periodic'})) )
                warning('CHEBOPPREF:PARSEDISCRETIZATION', ...
                    ['''TRIGCOLLOC''/''PERIODIC'' is deprecated. \n' ...
                    'Please use ''VALUES''/@trigcolloc.']);
                val = @trigcolloc;
            end

        end

        function val = parseHappinessCheck(val)
        %PARSEHAPPINESSCHECK    Allow different syntax for specifying
        %                       happinessCheck.

            % handle:
            if ( any(strcmpi(val, {'classic', 'classicCheck'})) )
                val = @classicCheck;

            elseif ( any(strcmpi(val, {'plateau', 'plateauCheck'})) )
                val = @plateauCheck;

            elseif ( any(strcmpi(val, {'strict', 'strictCheck'})) )
                val = @strictCheck;

            elseif ( any(strcmpi(val, {'loose', 'looseCheck'})) )
                val = @looseCheck;

            elseif ( any(strcmpi(val, {'happiness', 'happinessCheck'})) )
                val = @happinessCheck;

            elseif ( any(strcmpi(val, {'linopV4', 'linopV4Check'})) )
                val = @linopV4Check;

            end

        end

        function val = parseIVPsolver(val)
        %PARSEIVPSOLVER   Allow different syntax for specifying the IVPsolver.

            % Check whether we got pref.ivpSolver = @ode113/@ode45/@ode15s, that
            % is, a function handle, but not the CHEBFUN overload of it.
            if ( isa(val, 'function_handle') && ...
                    any(strcmpi(func2str(val), {'ode113', 'ode15s', 'ode45'})) )
                val = eval(['@chebfun.', func2str(val)]);

            % Check whether we got a string argument, e.g.
            % pref.ivpSolver = 'ode113'.
            elseif ( any(strcmpi(val, {'ode113', 'ode15s', 'ode45'})) )
                val = eval(['@chebfun.', val]);
            end
        end

    end

end

author: chebfun

repo: chebfun


text:
% Test file for ADCHEBFUN CUMSUM, DIFF, MEAN and SUM

function pass = test_cumsumDiffSumMean

% List of trigonometric functions to test.
funcList = {@diff, @(u)diff(u, 2), @(u)diff(u, 4), ...
                 @sum, @(u) sum(u, -.25, .8), ...
                 @cumsum, @(u)cumsum(u,2), ...
                 @mean, ...
                 @(u) deriv(u, 0:0.01:0.1, 2)};

% Tolerance for Taylor testing
tolOrder = 1e-2;
tolDiff = 1e-12;
% Initialise vector with pass information
pass = zeros(2, numel(funcList));

% Do the tests.
for k = 1:numel(funcList)
    % Call the valueTesting method, which also returns linearity information
    [err, lin] = adchebfun.valueTesting(funcList{k});

    % First, check that the computed function values match what we expect
    pass(1, k) = ( err == 0 );

    % Call the taylorTesting method
    [order1, order2, nDiff2] = adchebfun.taylorTesting(funcList{k});

    % We expect all elements of ORDER1 to be close to 1. Since the methods being
    % tested in this case are all linear, ORDER2 will be noise. However, since
    % the methods are indeed linear, we should expect nDiff2 to have values all
    % close to machine epsilon, which we can use to check for the correctness of
    % the derivative computed.
    pass(2,k) = ( (max(abs(order1 - 1)) < tolOrder) && ...
        (max(abs(nDiff2)) < tolDiff) );

    % Check that we received the correct linearity information
    pass(3, k) = ( lin == 1 );
end

end

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function out = cfg_run_runjobs(job)

% Initialise, fill, save and run a job with repeated inputs.
% To make use of possible parallel execution of independent jobs, all
% repeated jobs are filled first and (if successfully filled) run as one
% large job.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_run_runjobs.m 6514 2015-08-06 10:07:53Z volkmar $

rev = '$Rev: 6514 $'; %#ok

if isfield(job.save, 'savejobs')
    [p, n, e] = fileparts(job.save.savejobs.outstub);
    outfmt = strrep(fullfile(job.save.savejobs.outdir{1}, sprintf('%s_%%0%dd.m', n, ceil(log10(numel(job.inputs))+1))), '\', '\\');
end;
hjobs = cell(size(job.inputs));
sts = false(size(job.inputs));
out.jobfiles = {};
for cr = 1:numel(job.inputs)
    cjob = cfg_util('initjob', job.jobs);
    inp = cell(size(job.inputs{cr}));
    for ci = 1:numel(job.inputs{cr})
        fn = fieldnames(job.inputs{cr}{ci});
        inp{ci} = job.inputs{cr}{ci}.(fn{1});
    end;
    sts(cr) = cfg_util('filljob', cjob, inp{:});
    if sts(cr)
        [un, hjobs{cr}] = cfg_util('harvest', cjob);
    end;
    if isfield(job.save, 'savejobs')
        out.jobfiles{cr} = sprintf(outfmt, cr);
        cfg_util('savejob', cjob, out.jobfiles{cr});
    end;
    cfg_util('deljob', cjob);
end;
if all(sts)
    % keep all hjobs;
elseif any(sts) && strcmp(job.missing,'skip')
    hjobs = hjobs(sts);
else
    hjobs = {};
end
if ~isempty(hjobs)
    cjob = cfg_util('initjob', hjobs);
    cfg_util('run', cjob);
    out.jout = cfg_util('getalloutputs', cjob);
%    if isfield(job.save, 'savejobs')
%        [p n e] = fileparts(job.save.savejobs.outstub);
%        out.jobrun{1} = fullfile(p, [n '_run.m']);
%        cfg_util('saverun', cjob, out.jobrun{1});
%    end;
    cfg_util('deljob', cjob);
end;

author: spm


text:
% Implementation of the multiple extended object tracking algorithm based on the article
%
% "Linear-Time Joint Probabilistic Data Association for Multiple Extended Object Tracking (to appear)"
% S. Yang, K. Thormann, and M. Baum
% 2018 IEEE Sensor Array and Multichannel Signal Processing Workshop (SAM 2018), Sheffield, United Kingdom, 2018.
%
%
%
% Further information:
% http://www.fusion.informatik.uni-goettingen.de
% https://github.com/Fusion-Goettingen
%
% Source code written by Shishan Yang
% =============================

clc
close all
clear
dbstop error

scenario = 1; % turn and closely-spaced
% scenario = 2; % cross
clambda = 40; % clutter rate

% motion and measurement parameters used for multiplicative noise model
Cv = diag([10 10]);
Crw = diag([10 10 10 10]);
Cpw(:,:,1) = diag([0.02 1 1]);
Cpw(:,:,2) = diag([0.02 1 1]);
Ar = [1 0 1 0; 0 1 0 1; 0 0 1 0; 0 0 0 1];
Ap = eye(3);
H = [eye(2);zeros(2,2)]';

figure
hold on
box on
axis equal

[gt, meas,mlambda, xbound, ybound,cp,nr_timesteps] = getMeasGt(scenario,clambda,Cv);

% first guess
r = [-20 -250 10 10;-20 250 10 -10];
p = [0 30 30;0 15 15 ];

N = size(r,1);

Cr(:,:,1) = diag([900 900 10 10]);
Cp(:,:,1) = diag([.2 400 400]);

Cr(:,:,2) = diag([900 900 10 10]);
Cp(:,:,2) = diag([.02 100 100]);

%% plot first guess
for n = 1:N
    plot_extent([r(n,1:2) p(n,:)],'--','r',1);
end

for t = 1:nr_timesteps

    [r,p,Cr,Cp] = MEOT_JPDA(meas{t},r,p,Cr,Cp,cp,H,Cv,mlambda,clambda);

    %% Visulize
    if mod(t,3)==1

        pMeas = plot(meas{t}(1,:),meas{t}(2,:),'k.');
        for n = 1:N
            plotGT = plot_extent(gt(n,1:5,t),'-','k',1);
            plotEst = plot_extent([r(n,1:2),p(n,:)],'-','g',1);
        end
        pause(0.001)
    end

    %% prediction
    for n = 1:N
        r(n,:) = Ar*r(n,:)';
        p(n,:) = Ap*p(n,:)';
        Cr(:,:,n) = Ar*Cr(:,:,n)*Ar'+Crw;
        Cp(:,:,n) = Ap*Cp(:,:,n)*Ap'+Cpw(:,:,n);
    end
end
    legend([plotGT plotEst],{'Ground Truth','Estimates'})

author: Fusion-Goettingen

repo: ExtendedObjectTracking


text:
function normals = points2normals(points)
    % estimating a normal vector based on nearby 100 points
    % points is 3 * n matrix for n points

    if size(points,2)==3 && size(points,1)~=3
        points = points';
    end

    normals = lsqnormest(points, 100);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% functions from http://www.mathworks.com/matlabcentral/fileexchange/27804-iterative-closest-point

% Least squares normal estimation from point clouds using PCA
%
% H. Hoppe, T. DeRose, T. Duchamp, J. McDonald, and W. Stuetzle.
% Surface reconstruction from unorganized points.
% In Proceedings of ACM Siggraph, pages 71:78, 1992.
%
% p should be a matrix containing the horizontally concatenated column
% vectors with points. k is a scalar indicating how many neighbors the
% normal estimation is based upon.
%
% Note that for large point sets, the function performs significantly
% faster if Statistics Toolbox >= v. 7.3 is installed.
%
% Jakob Wilm 2010

function n = lsqnormest(p, k)
m = size(p,2);
n = zeros(3,m);

v = ver('stats');
if str2double(v.Version) >= 7.5
    neighbors = transpose(knnsearch(transpose(p), transpose(p), 'k', k+1));
else
    neighbors = k_nearest_neighbors(p, p, k+1);
end

for i = 1:m
    x = p(:,neighbors(2:end, i));
    p_bar = 1/k * sum(x,2);

    P = (x - repmat(p_bar,1,k)) * transpose(x - repmat(p_bar,1,k)); %spd matrix P
    %P = 2*cov(x);

    [V,D] = eig(P);

    [~, idx] = min(diag(D)); % choses the smallest eigenvalue

    n(:,i) = V(:,idx);   % returns the corresponding eigenvector
end

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

% Program to find the k - nearest neighbors (kNN) within a set of points.
% Distance metric used: Euclidean distance
%
% Note that this function makes repetitive use of min(), which seems to be
% more efficient than sort() for k < 30.

function [neighborIds,neighborDistances] = k_nearest_neighbors(dataMatrix, queryMatrix, k)

numDataPoints = size(dataMatrix,2);
numQueryPoints = size(queryMatrix,2);

neighborIds = zeros(k,numQueryPoints);
neighborDistances = zeros(k,numQueryPoints);

D = size(dataMatrix, 1); %dimensionality of points

for i=1:numQueryPoints
    d=zeros(1,numDataPoints);
    for t=1:D % this is to avoid slow repmat()
        d=d+(dataMatrix(t,:)-queryMatrix(t,i)).^2;
    end
    for j=1:k
        [s,t] = min(d);
        neighborIds(j,i)=t;
        neighborDistances(j,i)=sqrt(s);
        d(t) = NaN; % remove found number from d
    end
end

author: jianxiongxiao

repo: ProfXkit


text:
function v = value( x, data )
global cvx___
if nargin == 1,
    data = cvx___.x;
end
nx = size( data, 1 );
nb = size( x.basis_, 1 );
if nx < nb,
    data( end + 1 : nb, : ) = NaN;
elseif nx > nb,
    data( nb + 1 : end, : ) = [];
end
v = cvx_reshape( data.' * x.basis_, x.size_ );
if any( x.size_ == 1 ), v = full( v ); end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


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

clear all;
%close all;
clf

num=1;
den=conv([1 0],[1 2 1]);
rlocus(num,den);
title('Fig. 6.14(b)');
grid;

author: Sable

repo: mcbench-benchmarks


text:
classdef testInverseSymmetricFourthOrderTensor < testInverseFourthOrderTensor

    methods (Access = protected)
        function createRandomFourthOrderTensor(obj)
            obj.tensor = Stiffness3DTensor;
            obj.tensor.createRandomTensor();
        end
    end

    methods (Static, Access = protected)

        function Id = computeIdentityTensor(I,i,j,k,l)
            Id = 0.5*(I(i,k)*I(j,l) + I(i,l)*I(j,k));
        end
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function spy(disc, varargin)
%TODO: Either remove this method, or uncomment lines below and ensure this
%method works. The following is broken:
%   D = chebop(@(x,u) diff(u));
%   spy(chebcolloc2(linop(D)))

error('CHEBFUN:OPDISCRETIZATION:spy:noSupported', ...
    'OPDISCRETIZATION currently does not support spy().');
end
% %SPY    Visualize a LINOP.
% %   SPY(DISC) creates a picture of the nonzero pattern of the given LiNOP
% %   discretization DISC. Block boundaries are indicated by gray lines, and
% %   side condition rows are marked off by dashed lines (boundary and
% %   continuity conditions).
% %
% %   SPY(DISC,DIM) uses the dimension vector DIM to override the dimension
% %   property of DISC.
% %
% %   SPY(DISC,S) or SPY(DISC,DIM,S) allows modification of the plot
% %   attributes, as with the built in method.
% %
% % See also LINOP.SPY.
%
% % Copyright 2017 by The University of Oxford and The Chebfun Developers.
% % See http://www.chebfun.org for Chebfun information.
%
% % Obtain domain information.
% dom = disc.domain;
%
% % Parse the inputs:
% if ( ( nargin > 1 ) && isnumeric(varargin{1}) )
%     dim = varargin{1};
%     if ( numel(dim) == 1 )
%         % Scalar expand the dimension here.
%         dim = repmat(dim, 1, length(dom) - 1);
%     end
%
%     disc.dimension = dim;
%     varargin = varargin(2:end);
% end
% dim = disc.dimension;
%
% plotOpt = varargin;
%
% % Override hold state.
% holdState = ishold;
%
% % Check whether we need to derive continuity conditions.
% L = disc.source;
% if ( isempty(L.continuity) )
%     L = deriveContinuity(L, dom);
% end
% disc.source = L;
%
% % Spy the matrix, with a useful label.
% data = matrix(disc);
% spy(data, plotOpt{:}), hold on
% s =  sprintf('%i,', disc.dimension);    % list of sizes
% s = [ 'discretization = [', s(1:end-1), ']' ];
% xlabel(s)
%
% % Find the number of constraints and continuity conditions
% nbc = length(L.constraint);
% ncon = length(L.continuity);
% numInts = numel(dom) - 1;
%
% % Find all block sizes, substituting in the discretization size for Inf.
% [m, n] = blockSizes(L);
% m(isinf(m)) = sum(dim);
% n(isinf(n)) = sum(dim) + numInts*disc.dimAdjust(isinf(n));
%
% % Draw vertical block boundaries.
% cscol = cumsum(n(1,:));
% coldiv = cscol(1:end-1) + 1/2;
% rowmax = size(data, 1);
% plot([coldiv ; coldiv], [ 0; rowmax+1]*ones(size(coldiv)), 'color', [.6 .6 .6])
%
% % Draw horizontal block boundaries. Account for the down-sampling of each row.
% csrow = cumsum( m(:, 1)');                          % remove down-sampling
% rowdiv = csrow(1:end - 1) + 1/2;                    % boundary after each block
% rowdiv = nbc + ncon + rowdiv;                       % offset from top rows
% colmax = size(data, 2);
% % plot([0; colmax+1]*ones(size(rowdiv)), [rowdiv; rowdiv], 'color', [.6 .6 .6])
%
% % Draw horizontal BC and continuity boundaries.
% y = nbc + 1/2;
% plot([0; colmax+1], [y; y], '--', 'color', [.6 .6 .6])
% if ( ncon > 0 )
%     plot([0; colmax + 1], [y + ncon; y + ncon], '--', 'color', [.6 .6 .6])
% end
%
% % Reset hold state.
% if ( ~holdState )
%     hold off
% end
%
% end
%


text:
function savesurfpoly(v,f,holelist,regionlist,p0,p1,fname,forcebox)
%
% savesurfpoly(v,f,holelist,regionlist,p0,p1,fname)
%
% save a set of surfaces into poly format (for tetgen)
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
% date: 2007/11/21
%
% input:
%      v: input, surface node list, dimension (nn,3)
%         if v has 4 columns, the last column specifies mesh density near each node
%      f: input, surface face element list, dimension (be,3)
%      holelist: list of holes, each hole is represented by an internal point
%      regionlist: list of regions, similar to holelist
%      p0: coordinate of one of the end of the bounding box
%      p1: coordinate for the other end of the bounding box
%      fname: output file name
%      forcebox: non-empty: add bounding box, []: automatic
%                if forcebox is a 8x1 vector, it will be used to
%                specify max-edge size near the bounding box corners
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%
dobbx=0;
if(nargin>=8)
	dobbx=~isempty(forcebox) & all(forcebox);
end

if(~iscell(f) & size(f,2)==4)
    faceid=f(:,4);
    f=f(:,1:3);
end

if(~iscell(f))
    edges=surfedge(f);
else
    edges=[];
end
bbxnum=0;

nodesize=[];
if(size(v,2)==4)
   nodesize=v(:,4);
   v=v(:,1:3);
end
node=v;
loopid=[];
if(~isempty(edges))
    loops=extractloops(edges);
    if(length(loops)<3)
        error('degenerated loops detected');
    end
    seg=[0,find(isnan(loops))];
    segnum=length(seg)-1;
    newloops=[];
    for i=1:segnum
       if(seg(i+1)-(seg(i)+1)==0) continue; end
       newloops=[newloops nan bbxflatsegment(node,loops(seg(i)+1:seg(i+1)-1))];
    end
    loops=[newloops nan];

    seg=[0,find(isnan(loops))];
    segnum=length(seg)-1;
    bbxnum=6;
    loopcount=zeros(bbxnum,1);
    loopid=zeros(segnum,1);
    for i=1:segnum     % walk through the edge loops
        subloop=loops(seg(i)+1:seg(i+1)-1);
        if(isempty(subloop)) continue; end
        boxfacet=find(sum(abs(diff(v(subloop,:))))<1e-8); % find a flat loop
        if(length(boxfacet)==1)   % if the loop is flat along x/y/z dir
            bf=boxfacet(1);    % no degeneracy allowed
            if(sum(abs(v(subloop(1),bf)-p0(bf)))<1e-2)
                loopcount(bf)=loopcount(bf)+1;
                v(subloop,bf)=p0(bf);
                loopid(i)=bf;
            elseif(sum(abs(v(subloop(1),bf)-p1(bf)))<1e-2)
                loopcount(bf+3)=loopcount(bf+3)+1;
                v(subloop,bf)=p1(bf);
                loopid(i)=bf+3;
            end
        end
    end
end

if(dobbx & isempty(edges))
    bbxnum=6;
    loopcount=zeros(bbxnum,1);
end

if(dobbx|~isempty(edges))
    nn=size(v,1);
    boxnode=[p0;p1(1),p0(2:3);p1(1:2),p0(3);p0(1),p1(2),p0(3);
              p0(1:2),p1(3);p1(1),p0(2),p1(3);p1;p0(1),p1(2:3)];
    boxelem=[
        4 nn nn+3 nn+7 nn+4;   % x=xmin
        4 nn nn+1 nn+5 nn+4;   % y=ymin
        4 nn nn+1 nn+2 nn+3;   % z=zmin
        4 nn+1 nn+2 nn+6 nn+5; % x=xmax
        4 nn+2 nn+3 nn+7 nn+6; % y=ymax
        4 nn+4 nn+5 nn+6 nn+7];% z=zmax

    node=[v;boxnode];
end

node=[(0:size(node,1)-1)',node];

fp=fopen(fname,'wt');
fprintf(fp,'#node list\n%d 3 0 0\n',length(node));
fprintf(fp,'%d %f %f %f\n',node');

if(~iscell(f))
    fprintf(fp,'#facet list\n%d 1\n',length(f)+bbxnum);
    elem=[3*ones(length(f),1),f-1,ones(length(f),1)];
    if(~isempty(elem))
	fprintf(fp,'1 0\n%d %d %d %d %d\n',elem');
    end
else % if the surface is recorded as a cell array
    totalplc=0;
    for i=1:length(f)
        if(~iscell(f{i}))
            totalplc=totalplc+size(f{i},1);
        else
            totalplc=totalplc+size(f{i}{1},1);
        end
    end
    fprintf(fp,'#facet list\n%d 1\n',totalplc+bbxnum);
    for i=1:length(f)
        plcs=f{i};
        faceid=-1;
        if(iscell(plcs)) % if each face is a cell, use plc{2} for face id
            if(length(plcs)>1)
                faceid=plcs{2};
            end
            plcs=plcs{1};
        end
        for row=1:size(plcs,1);
         plc=plcs(row,:);
         if(any(isnan(plc))) % we use nan to separate outter contours and holes
            holeid=find(isnan(plc));
            if(faceid>0)
                fprintf(fp,'%d %d %d\n%d',length(holeid)+1,length(holeid),faceid,holeid(1)-1);
            else
                fprintf(fp,'%d %d\n%d',length(holeid)+1,length(holeid),holeid(1)-1);
            end
            fprintf(fp,'\t%d',plc(1:holeid(1)-1)-1);
            fprintf(fp,'\t1\n');
            for j=1:length(holeid)
                if(j==length(holeid))
                    fprintf(fp,'%d',length(plc(holeid(j)+1:end)));
                fprintf(fp,'\t%d',plc(holeid(j)+1:end)-1);
                else
                    fprintf(fp,'%d',length(plc(holeid(j)+1:holeid(j+1)-1)));
                    fprintf(fp,'\t%d',plc(holeid(j)+1:holeid(j+1)-1)-1);
                end
                fprintf(fp,'\t1\n');
            end
            for j=1:length(holeid)
                if(j==length(holeid))
                    fprintf(fp,'%d %f %f %f\n',j,mean(node(plc(holeid(j)+1:end),2:4)));
                else
                    fprintf(fp,'%d %f %f %f\n',j,mean(node(plc(holeid(j)+1:holeid(j+1)-1),2:4)));
                end
            end
         else
	    	if(faceid>0)
                fprintf(fp,'1 0 %d\n%d',faceid,length(plc));
            else
                fprintf(fp,'1 0\n%d',length(plc));
            end
    		fprintf(fp,'\t%d',plc-1);
            fprintf(fp,'\t1\n');
         end
        end
    end
end

if(dobbx|~isempty(edges))
    for i=1:bbxnum
        fprintf(fp,'%d %d 1\n',1+loopcount(i),loopcount(i));
        fprintf(fp,'%d %d %d %d %d\n',boxelem(i,:));
        if(~isempty(edges) & loopcount(i) &~isempty(find(loopid==i)))
            endid=find(loopid==i);
            for k=1:length(endid)
                j=endid(k);
                subloop=loops(seg(j)+1:seg(j+1)-1);
                fprintf(fp,'%d ',length(subloop));
                fprintf(fp,'%d ',subloop-1);
                fprintf(fp,'\n');
            end
            for k=1:length(endid)
                j=endid(k);
                subloop=loops(seg(j)+1:seg(j+1)-1);
                fprintf(fp,'%d %f %f %f\n',k,internalpoint(v,subloop)); %mean(v(subloop,:)));
            end
        end
    end
end

if(size(holelist,1))
        fprintf(fp,'#hole list\n%d\n',size(holelist,1));
        for i=1:size(holelist,1)
                fprintf(fp,'%d %f %f %f\n', i, holelist(i,:));
        end
else
	fprintf(fp,'#hole list\n0\n');
end

if(size(regionlist,1))
	fprintf(fp,'#region list\n%d\n',size(regionlist,1));
	for i=1:size(regionlist,1)
		fprintf(fp,'%d %f %f %f %d\n', i, regionlist(i,:),i);
	end
end
fclose(fp);

if(~isempty(nodesize))
	if(size(nodesize,1)+size(forcebox(:),1)==size(node,1))
		nodesize=[nodesize;forcebox(:)];
	end
	fid=fopen(regexprep(fname,'\.poly$','.mtr'),'wt');
	fprintf(fid,'%d 1\n',size(nodesize,1));
	fprintf(fid,'%f\n',nodesize);
	fclose(fid);
end


function dataout=modeWPE_NLMS_BSS(nummics, numrefs, datain)
%
% Perform wpe, nlms, and bss.
% nummics:              no. of mic channels
% numrefs:              no. of reference channels
% datain:               input data
% dataout:              output data
%

fs=16000;

%% perform wpe
addpath('wpe_v1.33');

mic=datain(:, 1:2);
ref=datain(:, 3);

cfgs='wpe_v1.33/settings/local.m';
early=wpe(mic, cfgs);

%% perform nlms
addpath('Speex-AEC-matlab-master');

%    Usage:
%
%       speex_mdf_out = speex_mdf(Fs, u, d, filter_length, frame_size, dbg_var_name);
%
%       Fs                  sample rate
%       u                   speaker signal, column vector in range [-1; 1]
%       d                   microphone signal, column vector in range [-1; 1]
%       filter_length       typically 250ms, i.e. 4096 @ 16k FS
%                           must be a power of 2
%       frame_size          typically 8ms, i.e. 128 @ 16k Fs
%                           must be a power of 2
%       dbg_var_name        internal state variable name to trace.
%                           Default: 'st.leak_estimate'.
%
%    Jonathan Rouach <jonr@waves.com>
%
filter_length=512;
frame_size=256;
nearend1=speex_mdf(fs, ref, early(:, 1), filter_length, frame_size);
nearend2=speex_mdf(fs, ref, early(:, 2), filter_length, frame_size);
nearend=[nearend1.e, nearend2.e];

%% perform stft
addpath('stft2');
% fft size
fftsize=512;
stftshift=fftsize/2;

M=nummics;
N=M;

Xtf=cell(M, 1);
for m=1:M
    Xtf{m}=stft(nearend(:, m), stftshift, fftsize, false);
end
[K, T]=size(Xtf{1});

Ytf=cell(M, 1);
for m=1:M
    Ytf{m}=zeros(K, T);
end

%% params
% forgetting factor for bss
BF_FORGET=0.999;
% the shape parameter of the source prior
GAMMA=0.2;
%
% used to keep stable
%
VAR_BIAS=0.01;
STABLE_EPS=1e-3;
BF_DIAGLOAD=1e-6;

%% space for bss
% the weighted correlation matrices
C1=cell(K, 1);
C2=cell(K, 1);
for k=1:K
    C1{k}=STABLE_EPS*eye(M, M);
    C2{k}=STABLE_EPS*eye(M, M);
end

% demixing matrices
Demix=cell(K, 1);
for k=1:K
    Demix{k}=eye(N, M);
end

%% perform bss
for tau=1:T
    %% perform bss
    Bssout=zeros(K, M);

    %
    % calculate nonlinearity
    %
    phi1=0;
    phi2=0;

    for k=1:K
        x=zeros(M, 1);
        for m=1:M
            x(m)=Xtf{m}(k, tau);
        end

        y=Demix{k}*x;
        % output data
        Bssout(k, :)=y.';

        phi1=phi1+abs(y(1))^2;
        phi2=phi2+abs(y(2))^2;
    end

    phi1=(1-BF_FORGET)*(phi1+VAR_BIAS)^((GAMMA-2)/2);
    phi2=(1-BF_FORGET)*(phi2+VAR_BIAS)^((GAMMA-2)/2);

    % update the demixing matrices
    for k=1:K
        %
        % accumulate the weighted correlation
        %
        x=zeros(M, 1);
        for m=1:M
            x(m)=Xtf{m}(k, tau);
        end

        C1{k}=BF_FORGET*C1{k}+phi1*(x*x');
        C2{k}=BF_FORGET*C2{k}+phi2*(x*x');

        %
        % solve gev problem
        %
        [Ev, Ed]=eig(C2{k}+BF_DIAGLOAD*eye(M, M), C1{k}+BF_DIAGLOAD*eye(M, M));
        if Ed(1, 1)>=Ed(2, 2)
            e1=Ev(:, 1);
            e2=Ev(:, 2);
        else
            e1=Ev(:, 2);
            e2=Ev(:, 1);
        end

        D=[e1'; e2'];

        %
        % solve the scaling ambiguity
        %
        A=inv(D);

        if abs(A(1, 1))>=abs(A(2, 1))
            a1=A(1, 1);
        else
            a1=A(2, 1);
        end

        if abs(A(2, 2))>=abs(A(1, 2))
            a2=A(2, 2);
        else
            a2=A(1, 2);
        end

        D=diag([a1; a2])*D;
        Demix{k}=D;
    end

    for m=1:M
        Ytf{m}(:, tau)=Bssout(:, m);
    end
end

%% perform istft and output signal
dataout=zeros(dataLength(T, stftshift, fftsize ), N);
for n=1:N
    dataout(:, n)=istft(Ytf{n}, stftshift, false);
end

end


text:
function nlp = Time(nlp, varargin)

    plant = nlp.Plant;
    T = SymVariable('t',[2,1]);

    tf = SymFunction(['time_' plant.Name],T(2),{T});
    addNodeCost(nlp,tf,{'T'},'first');
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function [dv, dvArr, flyBy1DVVect, flyBy2DVVect, flyBy1Rp, flyBy2Rp, ...
          xferOrbitP1, xferOrbitP2, xferOrbitP3, flyBy1OrbitIn, flyBy1OrbitOut, ...
          flyBy2OrbitIn, flyBy2OrbitOut] = ...
          findOptimalTwoFlyByTraj(x, departBodyInfo, flyby1BodyInfo, flyby2BodyInfo, arrivalBodyInfo, gmuXfr)
%findOptimalTwoFlyByTraj Summary of this function goes here
%   Detailed explanation goes here
    if(length(x)~=4)
        error('Length of x in findOptimalTwoFlyByTraj must be 4');
    end

    departDate = x(1);
    phase1TOF = x(2);
    phase2TOF = x(3);
    phase3TOF = x(4);

    xDates = [departDate,
              departDate + phase1TOF,
              departDate + phase1TOF + phase2TOF,
              departDate + phase1TOF + phase2TOF + phase3TOF];

    preFunc1 = @(arrivalUT, departUT) findOptimalDepartureArrivalObjFunc(arrivalUT, departUT, departBodyInfo, flyby1BodyInfo, gmuXfr, 'departDVRadioBtn');
    objFuncDepart = @(x) preFunc1(x(2), x(1));

    objFuncFlyby1DV = @(x) flybyDVObjFunc([xDates(1) xDates(2) xDates(3)], departBodyInfo, flyby1BodyInfo, flyby2BodyInfo, gmuXfr);
    objFuncFlyby2DV = @(x) flybyDVObjFunc([xDates(2) xDates(3) xDates(4)], flyby1BodyInfo, flyby2BodyInfo, arrivalBodyInfo, gmuXfr);

    preFunc2 = @(arrivalUT, departUT) findOptimalDepartureArrivalObjFunc(arrivalUT, departUT, flyby2BodyInfo, arrivalBodyInfo, gmuXfr, 'arrivalDVRadioBtn');
    objFuncArrival = @(x) preFunc2(x(4), x(3));

    dv1 = objFuncDepart(xDates);
    [dv2, flyBy1DVVect, flyBy1Rp, xferOrbitP1, xferOrbitP2, flyBy1OrbitIn, flyBy1OrbitOut] = objFuncFlyby1DV(xDates);
    [dv3, flyBy2DVVect, flyBy2Rp, xferOrbitP2, xferOrbitP3, flyBy2OrbitIn, flyBy2OrbitOut] = objFuncFlyby2DV(xDates);
    dv4 = objFuncArrival(xDates);
    dvArr = [dv1, dv2, dv3, dv4];

    dv = dv1 + dv2 + dv3 + dv4;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [xVals,tVals,dxdtVals,exitCode,numRejections]=RosenbrockAdaptiveOverRange(xStart,tSpan,f,order,JacobianFun,dfdtFun,initStepSize,RelTol,AbsTol,maxSteps)
%%ROSENBROCKADAPTIVEOVERRANGE Integrate an ordinary differential equation
%                   using a modified Rosenbrock method with an adaptive
%                   step size. That is, integrating dx/dt=f(x,t) given
%                   initial conditions (xStart,tSpan(1)). Rosenbrock
%                   methods are better then Runge-Kutta methods when the
%                   differential equations are stiff. However, Rosenbrock
%                   methods require a derivative. More information
%                   on the algorithms used with different orders is given
%                   in the comments to the function RosenbrockStep.
%
%INPUTS: xStart The NX1 state vector at time tSpan(1).
%         tSpan  A 2X1 or 1X2 vector where tSpan(1) is the starting time
%                and tSpan2 is the desired stopping time for the
%                integration.
%             f  The function handle for f(x,t)=dxdt over which integration
%                is to be performed.
%          order The order of the main Rosenbrock routine to use. See the
%                function RosenbrockStep for possible values and algorithms
%                used as this function calls RosenbrockStep. If omitted or
%                an empty matrix is passed, the default value of 2 is used.
%   initStepSize An optional initial step size (in t) to use for the
%                integration. If omitted or an empty matrix is passed, an
%                ad-hoc method described below is used to find an initial
%                step size.
%    JacobianFun This is either a function handle having the form
%                JacobianFun(x,t) that provides the Jacobian of the
%                function f --the derivatives with respect to the parameter
%                x-- or this is a matrix that is the constant Jacobian
%                matrix. If this parameter is omitted or an empty matrix is
%                passed, then the NXN Jacobian matrix is computed using the
%                numjac function with the tolerance set to AbsTol.
%        dfdtFun This is either a function handle having the form
%                dfdtFun(x,t) that provides the derivative of f with
%                respect to the scalar parameter t, or this is a constant
%                NX1 vector for the derivative of f with respect to t. In
%                autonomous problems, a zero matrix should be passed. If
%                omitted, a central difference formula is used to
%                numerically approximate the derivative vector.
%   initStepSize An optional initial step size (in t) to use for the
%                integration. If omitted or an empty matrix is passed, an
%                ad-hoc method described below is used to find an initial
%                step size.
%         RelTol The maximum relative error tolerance allowed, a
%                positive scalar (its use is explained in more detail
%                below). If omitted or an empty matrix is passed, the
%                default value of 1e-3 is used.
%         AbsTol The absolute error tolerance allowed, a positive scalar,
%                of the same for all components of x, or a positive NX1
%                vector. If omitted or an empty matrix is passed, the
%                default value of 1e-6 is used.
%       maxSteps The maximum allowable number of steps to perform the
%                integration. If omitted, the default of 4096 is used.
%
%OUTPUTS:xVals The NXnumSteps set of values of x along the path.
%                  xVals(:,1) is xStart and xVals(:,end) is the value at
%                  final time tSpan(2). If the integration failed, e.g. due
%                  to encountering a NaN or being unable to get a sufficient
%                  step size, an empty matrix is returned.
%            tVals A numStepsX1 vector of values of t corresponding to
%                  the values of x in xVals. tVals(1) is equal to tSpan(1)
%                  and tVals(end) is equal to tSpan(2). If integration
%                  fails, an empty matrix is returned.
%         dxdtVals A numStepsX1 array containing deivatives of x evaluated
%                  at the times in tVals. The derivatives are just the
%                  result of evaluating f(x,t) at each point in
%                  (xVals,tVals). This combined with xVals and tVals could
%                  be used in functions to perform Hermite interpolation.
%                  If integration fails, an empty matrix is returned.
%         exitCode A code indicating how the algorithm terminated. Possible
%                  values are
%                  0: Integration was successful.
%                  1: Unable to get a small enough step size.
%                  2: Maximum number of steps reaced without completion.
%                  3: Non-finite number encountered.
%    numRejections The number of times a stepsize hypothesis was rejected.
%
%The same type of adaptive stepsize control is used with the Rosenbrock
%algorithm (a diagonal implicit Runge-Kutta routine) as with standard
%Runge-Kutta methods. Thus, the comments to the stepsize adaptation routine
%in RKAdaptiveOverRange explain how the stepsize is adjusted here. The onyl
%notable change is that a certain matrix has to be inverted for each of the
%Rosenbrock methods. The matrix depends on the stepsize. Thus, if the
%matrix is singular, this function has an additional rejection of the step
%and the stepsize is halved.
%
%The central difference formula used for the numeric derivative of f with
%respect to t if dfdtFun is not provided is taken from Chapter 4.1 of [1],
%whereby the three-point midpoint formula is used.
%
%REFERENCES:
%[1] R. L. Burden and J. D. Faires, Numerical Analysis, 9th ed. Boston, MA:
%    Brooks/ Cole, 2011.
%
%February 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<10)
    maxSteps=4096;
end

if(nargin<9||isempty(AbsTol))
    AbsTol=1e-6;
end

if(nargin<8||isempty(RelTol))
   RelTol=1e-3;
end

xDim=size(xStart,1);
tStart=tSpan(1);
tEnd=tSpan(2);
tDiff=diff(tSpan);
tDiffMag=abs(tDiff);
deltaTSign=sign(tDiff);
%The sign is separated out so that one can use the algorithm running
%backwards as well as forwards in time.

%The maximum step size is arbitraily set to 1/3 the total distance.
maxStepSize=tDiffMag/3;

if(nargin<7||isempty(initStepSize))
    %If no initial step size is given, then just use the smallest uniform
    %step size for the given value of maxSteps.
    initStepSize=tDiffMag/maxSteps;
elseif(initStepSize<eps(tStart))
   error('The initial step size provided is too small')
end

if(nargin<6)
    %Tell it to perform numerical differentiation for dfdt.
    dfdtFun=[];
end

if(nargin<4||isempty(order))
    order=2;
end

%The lowest-order formula is used for the convergence rate. The isFSAL flag
%indicates whether k(:,end) is the value of f evaluated at the next step.
[orders,isFSAL]=RungeKStep(order);
RKOrder=min(orders);

%Allocate the maximum possible space needed for the return parameters.
xVals=zeros(xDim,maxSteps);
dxdtVals=zeros(xDim,maxSteps);
tVals=zeros(maxSteps,1);

xVals(:,1)=xStart;
tVals(:,1)=tStart;
dxdtVals(:,1)=f(xStart,tStart);

if(nargin<5||isempty(JacobianFun))
    %Indicate that the jacobian should be determined numerically
    JacobianType=0;
    FAC=[];%A pointer to store results for the jacobian function between
           %calls to make it faster.

    %The numjac function requires a minimum tolerance that is a vector.
    if(isscalar(AbsTol))
       AbsTol=repmat(AbsTol,[xDim,1]);
    end
elseif(isa(JacobianFun,'function_handle'))
    %A function handle is provided for the Jacobian.
    JacobianType=1;
else
    JacobianType=2;
    %A constant Jacobian matrix is provided.
	J=JacobianFun;
end

deltaTMag=initStepSize;
numRejections=0;
for curStep=2:maxSteps
    xCur=xVals(:,curStep-1);
    tCur=tVals(curStep-1);
    dxdtCur=dxdtVals(:,curStep-1);

    %The absolute value of the minimum allowable step size. It is set so
    %that the step must makes something of a difference compared to the
    %numerical precision.
    deltaTMinMag=2^4*eps(tCur);
    %Since th minimum step size changes every loop, this makes sure that
    %deltaT does not go beneath it just because the loop changed.
    deltaTMag=max(deltaTMag,deltaTMinMag);
    deltaT=deltaTMag*deltaTSign;

    %If we would have overstepped the end, change deltaT to be the end.
    %Then, we can terminate at the end.
    if(deltaTSign*(tCur+deltaT)>tEnd*deltaTSign)
        deltaTMag=abs(tEnd-tCur);
        deltaT=deltaTMag*deltaTSign;
        %Allow the last step to be very small.
        deltaTMinMag=min(deltaTMinMag,deltaTMag);
    end

    if(JacobianType==0)
        %Numerically find the Jacobian. The numjac function is used.
        fRev=@(t,x)f(x,t);
        [J,FAC] = numjac(fRev,tCur,xCur,dxdtCur,AbsTol,FAC,false);
    elseif(JacobianType==1)
        %User-provided function
        J=JacobianFun(xCur,tCur);
    end

    %The first time the choice in step size fails, it is adjusted the
    %"optimal" way in the Runge-Kutta-Fehlberg method. Additional times,
    %the step size is just halved in the hope that it will reach an
    %accepted value more quickly.
    failedReducingStepSize=false;
    moveOnToNextStep=false;
    while(moveOnToNextStep==false)
        %If an explicit function (or matrix) for the derivative with
        %respect to time is provided.
        if(isempty(dfdtFun))
            %Just use a central finite-difference formula with a set being
            %a quarter of the current step size.
            dfdt=(f(xCur,tCur+deltaT/8)-f(xCur,tCur-deltaT/8))/(2*8);
        elseif(isa(dfdtFun,'function_handle'))
            dfdt=dfdtFun(xCur,tCur);
        else%It is a constant matrix.
            dfdt=dfdtFun;
        end

        [xPredMain,xPredSubsid,k]=RosenbrockStep(xCur,tCur,f,deltaT,dxdtCur,J,dfdt,order);
        %If the step failed because the W matrix was singular, reduce the
        %step size and the matrix should not stay singular.
        if(isempty(xPredMain))
            numRejections=numRejections+1;

            if(deltaTMag<deltaTMinMag)
                %If the step size got too small, then return.
                xVals=[];
                dxdtVals=[];
                tVals=[];
                exitCode=1;
                return;
            end

            %Just halve the step size to make W non-singular.
            deltaTMag=deltaTMag/2;
            deltaT=deltaTSign*deltaTMag;
            continue;
        end

        %Integration can only be over finite functions.
        if(any(~isfinite(xPredMain))||any(~isfinite(xPredSubsid)))
            xVals=[];
            tVals=[];
            dxdtVals=[];
            exitCode=3;
            return;
        end

        %The step-size adaptation is the same as in RKAdaptiveOverRange
        %The local error estimate. This must be transformed into a
        %combination relative/ absolute error term to determine whether
        %the step should be rejected.
        normFactor=max(max(abs(xPredMain),abs(xCur)),AbsTol/RelTol);
        theError=max(abs((xPredMain-xPredSubsid)./normFactor));
        if(theError>RelTol)
            %The step should be rejected.
            numRejections=numRejections+1;

            if(deltaTMag<deltaTMinMag)
                %If the step size got too small, then return.
                xVals=[];
                dxdtVals=[];
                tVals=[];
                exitCode=1;
                return;
            end

            if(failedReducingStepSize==false)
                failedReducingStepSize=true;

                %The Fehlberg step reduction (using the relative error).
                deltaTMag=max(deltaTMinMag, deltaTMag * max(0.1, 0.8*(RelTol/theError)^(1/RKOrder)));
            else
                %Just halve the step size.
                deltaTMag=deltaTMag/2;
            end
            deltaT=deltaTSign*deltaTMag;
        else
            %If a step is successful, then increase the step size for the
            %next step in the standard manner used with Runge-Kutta-
            %Fehlberg methods, but limit the maximum size of the increase
            %to a scale factor of 4. This avoid huge step sizes when the
            %predicted error is very small.
            deltaTMag=min(maxStepSize,deltaTMag*min(4,0.8*(RelTol/theError)^(1/RKOrder)));

            %Save the results from the step.
            xVals(:,curStep)=xPredMain;
            tVals(curStep)=tCur+deltaT;
            if(isFSAL)
                dxdtVals(:,curStep)=k(:,end);
            else
                dxdtVals(:,curStep)=f(xVals(:,curStep),tVals(curStep));
            end

            moveOnToNextStep=true;
        end
    end

    if(tVals(curStep)==tEnd)
        %If integration completed, reshape the outputs to match the actual
        %number of steps.
        xVals=xVals(:,1:curStep);
        tVals=tVals(1:curStep);
        dxdtVals=dxdtVals(:,1:curStep);
        exitCode=0;
        return
    end
end

%If we get here, then the maximum number of iterations was reached without
%reaching the end.
xVals=[];
dxdtVals=[];
tVals=[];
exitCode=2;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
filename='Bridge_hexahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP'; incrementFactor = 1;
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.05;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;


text:
function [MVMResult,TWAResult] = eval_mvm_twa(ecg,ann,fs)
%OVERVIEW This function evaluates MVM and TWA for ecg data sampled at 1000
% Hz and the following list of annotations, qrs onset, q, r, s, qrs offsett,
% t offset.
%
%   INPUTS          MANDATORY       DESCRIPTION
%                   ecg             N by M array of ecg data. Contains M
%                                   channels of ecg with N data points.
%
%                   fs              Sampling frequency for ecg, needs to 1000 Hz.
%
%                   OPTIONAL
%                   ann             Structure containing the following
%                                   annotations, qrs onset, q, r, s, qrs
%                                   offset, t offset. If annotations are
%                                   not available pass an empty variable
%                                   ([]) as input for ann. Annotations will
%                                   be generated.
%
%   OUTPUTS
%                   MVMResult                       Structure containing the results of MVM
%                                                   analysis. The fields in the structure
%                                                   are listed below.
%
%                   Structure.Field
%                   MVMResult.energyinband_array    Array containing the energy in the every 2-7
%                                                   beat interval computed in the beatquency domain, for each analysis window.
%
%                   MVMResult.sqi_array             signal quality index
%                                                   for each analysis window.
%
%                   MVMResult.heart_rate_est_arr    average heart rate
%                                                   estimate for each analysis window.
%
%                   TWAResult               Structure containing the
%                                           results of TWA analysis. The fields in the structure
%                                           are listed below.
%
%                   TWAResult.HR            Estimate of average heart rate for each analysis window.
%                   TWAResult.VAlt          Estimate of TWA amplitude for
%                                           each analysis window.
%                   TWAResult.VAlt_Sig      Estimate of TWA amplitude for
%                                           each analysis window which are statistically
%                                           significant compared to the noise threshold.
%                   TWAResult.Noise_Median  The median of the gamma
%                                           distribution used to model the noise.
%                   TWAResult.Noise_95      The 95th percentile of the
%                                           gamma distribution used to model the noise.
%                   TWAResult.VAltPt        The location on the S-T offset
%                                           segment corresponding to the maximum difference between
%                                           the average even and average odd beat.
%                   TWAResult.successful    Flag if set to 1 means TWA
%                                           analysis is successful.
%
%   REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Ismail Sadiq
%	COPYRIGHT (C) 2019
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information. The license may be found in
%       the Documents folder of the Physionet-Cardiovascular-Signal-Toolbox.
%

% compute minimum amplitude resolution
amplitude_res = diff(ecg);
amplitude_res = min(amplitude_res(amplitude_res > eps));
if (amplitude_res > 0.006) % if amplitude resolution > 6 uV give warning on accuracy of twa results.
    disp(['Warning: Amplitude resolution for input ecg is ' num2str(amplitude_res) '. This greater than the suggested minimum resolution of 6 uV. Results for MV analysis may not be reliable.']);
end

if (fs == 1000)

    % Add dependencies
    addpath(genpath('../../../../PhysioNet-Cardiovascular-Signal-Toolbox-master')); % Cardio vascular toolbox
    addpath(genpath('./Tools/')); % Helper functions for mvtoolbox

    % check if ann provided as input
    annfields = {'QRSon', 'Q', 'R', 'S', 'QRSoff', 'Toff'};
    if (~isempty(ann))
        ann_check = zeros(1,length(annfields));
        for ii = 1:length(ann_check)
            if (isfield(ann,annfields{1,ii}) && ~isempty(getfield(ann,annfields{1,ii})))
                ann_check(ii) = 1;
            end
        end

        % check if any field missing
        missingi = find(ann_check == 0);
        if (~isempty(missingi))
            anntxt = [];
            for ii = 1:length(missingi)
                anntxt = [anntxt annfields{1,missingi(ii)} ', '];
            end
            disp(['Warning: Missing annotations, ' anntxt(1:end-2)])
            % generate annotations if missing
            disp('Incomplete annotations provided as input. Generating missing annotations.')
            HRVparams = InitializeHRVparams('MVanalysis'); %ecg = double(ECG_signal(:,7))/max(double(ECG_signal(:,7)));%std(double(ECG_signal(:,7)));
            [qrs_pos,sign,en_thres] = jqrs(ecg,HRVparams);
            ECG_header.nsig = 1; ECG_header.freq = fs; ECG_header.nsamp = length(ecg);
            wavedet_config.setup.wavedet.QRS_detection_only = 0;
            % plot(ecg); hold on; scatter(qrs_pos, ecg(qrs_pos)); hold off;
            [ann,~,~] = wavedet_3D_detector(ecg, qrs_pos', ECG_header, wavedet_config );

            [detection_flg,ann] = improvfiducials(ann, fs, ecg);

            if(~detection_flg)
                disp('Warning: Unable to improve fiducial point detection.')
            end

        end

    else
        disp('Annotations not provided as input. Generating missing annotations.')
        HRVparams = InitializeHRVparams('MVanalysis'); %ecg = double(ECG_signal(:,7))/max(double(ECG_signal(:,7)));%std(double(ECG_signal(:,7)));
        [qrs_pos,sign,en_thres] = jqrs(ecg,HRVparams);
        ECG_header.nsig = 1; ECG_header.freq = fs; ECG_header.nsamp = length(ecg);
        wavedet_config.setup.wavedet.QRS_detection_only = 0;
        % plot(ecg); hold on; scatter(qrs_pos, ecg(qrs_pos)); hold off;
        [ann,~,~] = wavedet_3D_detector(ecg', qrs_pos', ECG_header, wavedet_config );

        [detection_flg,ann] = improvfiducials(ann, fs, ecg);

        if(~detection_flg)
            disp('Warning: Unable to improve fiducial point detection.')
        end

    end

    % figure(1); plot(ecg); hold on; % Check annotation if needed
    % stem(ann.QRSon, ecg(ann.QRSon)); stem(ann.Q, ecg(ann.Q));
    % stem(ann.R, ecg(ann.R)); stem(ann.S, ecg(ann.S));
    % stem(ann.QRSoff, ecg(ann.QRSoff)); stem(ann.Toff, ecg(ann.Toff)); hold off
    % legend('ECG','QRSon','Q','R','S','QRSoff','Toff')

    % Variable for storing mvm for qrs
    numofleads = size(ecg,2);
    energyinband_array = cell(1,numofleads);
    sqi_array = cell(1,numofleads);
    heart_rate_est_arr = cell(1,numofleads);
    % Compute Morph Var QRS
    disp(['Evaluating QRS MVM ...']);
    normalize = 1;
    segment_size = 5;
    for lead = 1:size(ecg,2)
        [energyinband,sqi,heart_rate_est] = ComputeMVM(ecg(:,lead),ann,fs,segment_size,normalize);
        energyinband_array(1,lead) = mat2cell(energyinband, size(energyinband,1), size(energyinband,2));
        sqi_array(1,lead) = mat2cell(sqi, size(sqi,1), size(sqi,2));
        heart_rate_est_arr(1,lead) = mat2cell(heart_rate_est, size(heart_rate_est,1), size(heart_rate_est,2));
    end
    MVMResult.energyinband_array = energyinband_array;
    MVMResult.sqi_array = sqi_array;
    MVMResult.heart_rate_est_arr = heart_rate_est_arr;

    % Twa analysis
    pause(2)
    disp('Evaluating TWA ...');
    TWAResult = ComputeTWA(ecg,ann,fs);

else
    disp(['Input ecg is sampled at ' num2str(fs) ' Hz, needs to be sampled at 1000 Hz. Unable to perform MV analysis.'])
    MVMResult = NaN; TWAResult = NaN;
end

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function ndx = subv2ind(siz, subv)
% SUBV2IND Like the built-in sub2ind, but the subscripts are given as row vectors.
% ind = subv2ind(siz,subv)
%
% siz can be a row or column vector of size d.
% subv should be a collection of N row vectors of size d.
% ind will be of size N * 1.
%
% Example:
% subv = [1 1 1;
%         2 1 1;
%         ...
%         2 2 2];
% subv2ind([2 2 2], subv) returns [1 2 ... 8]'
% i.e., the leftmost digit toggles fastest.
%
% See also IND2SUBV.

if isempty(subv)
  ndx = [];
  return;
end

if isempty(siz)
  ndx = 1;
  return;
end

[ncases ndims] = size(subv);

%if length(siz) ~= ndims
%  error('length of subscript vector and sizes must be equal');
%end

if all(siz==2)
  %rbits = subv(:,end:-1:1)-1; % read from right to left, convert to 0s/1s
  %ndx = bitv2dec(rbits)+1;
  twos = pow2(0:ndims-1);
  ndx = ((subv-1) * twos(:)) + 1;
  %ndx = sum((subv-1) .* twos(ones(ncases,1), :), 2) + 1; % equivalent to matrix * vector
  %ndx = sum((subv-1) .* repmat(twos, ncases, 1), 2) + 1; % much slower than ones
  %ndx = ndx(:)';
else
  %siz = siz(:)';
  cp = [1 cumprod(siz(1:end-1))]';
  %ndx = ones(ncases, 1);
  %for i = 1:ndims
  %  ndx = ndx + (subv(:,i)-1)*cp(i);
  %end
  ndx = (subv-1)*cp + 1;
end

%%%%%%%%%%%

function d = bitv2dec(bits)
% BITV2DEC Convert a bit vector to a decimal integer
% d = butv2dec(bits)
%
% This is just like the built-in bin2dec, except the argument is a vector, not a string.
% If bits is an array, each row will be converted.

[m n] = size(bits);
twos = pow2(n-1:-1:0);
d = sum(bits .* twos(ones(m,1),:),2);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function result=hl(x)

%HL computes the Hodges-Lehmann location estimate on the columns of x.
% The Hodges-Lehmann estimator is defined as
%           hl(x)=med {(x_i + x_j)/2} with 1<=i<j<=n
% It can resist about 29% outliers.
% If x is a matrix, the location estimate is computed on the columns of x. The
% result is then a row vector. If x is a row or a column vector,
% the output is a scalar.
%
% The behavior of the HL estimator in small samples is discussed in:
%   Rousseeuw, P.J. and Verboven, S. (2002),
%   "Robust estimation in very small samples",
%   Computational Statistics and Data Analysis, 40, 741-758.
%
% Required input argument:
%    x: either a data matrix with n observations in rows, p variables in columns
%       or a vector of length n.
%
% I/O:  result=hl(x);
%
% This function is part of LIBRA: the Matlab Library for Robust Analysis,
% available at:
%              http://wis.kuleuven.be/stat/robust.html
%
% Written by S. Verboven
% Last revision: 28/08/03 by N. Smets

[n,p]=size(x);

if n==1 & p==1
    result=x;     %when X is a one by one matrix, all location estimators must be equal to that matrix
    return
elseif n==1
    x=x';      %we only want to work with column vectors
    n=p;
    p=1;
end

if n==2
    result=mean(x,1);  % all location estimators must equal to the average for n=2
    return
end

for k=1:p
    m=0;
    y=zeros(1,n*(n-1)/2); %initializing help vector
    X=x(:,k);
    for i=1:n              %calculating all possible pairwise means
        for j=(i+1):n
            m=m+1;
            y(m)=(X(i)+X(j))/2;
        end
    end
    result(k)=median(y);
    y=[];                 %clear help vector
end

author: marianux

repo: ecg-kit

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


text:
function out = checkOption(optStruct,option,value)
% check for option in option list
%
% Syntax
%   out = checkOption(option_list,option_name)
%   out = checkOption(option_list,option_name,value)
%
% Input
%  option_list - Cell Array
%  option_name - String
%  option      - String
%  type        - class
%
% Output
%  out         - true / false
%
%% See also
% get_option set_option find_option

if ~isfield(optStruct,option)
  out = false;
elseif nargin == 2

  if islogical(optStruct.(option))
    out = optStruct.(option);
  end

else

  switch class(optStruct.(option))

    case 'char'

      out = strcmpi(optStruct.(option),value);

    otherwise %{'logical','double'}

      out = all(optStruct.(option) == value);

  end
end

author: mtex-toolbox

repo: mtex

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


text:
%SurfPointFeature  SURF point corner feature object
%
% A subclass of OrientedScalePointFeature for SURF features.
%
% Methods::
% plot         Plot feature position
% plot_scale   Plot feature scale
% distance     Descriptor distance
% ncc          Descriptor similarity
% match        Match features
% uv           Return feature coordinate
% display      Display value
% char         Convert value to string
%
% Properties::
% u             horizontal coordinate
% v             vertical coordinate
% strength      feature strength
% scale         feature scale
% theta         feature orientation [rad]
% descriptor    feature descriptor (vector)
% image_id      index of image containing feature
%
% Properties of a vector of SurfCornerFeature objects are returned as a vector.
% If F is a vector (Nx1) of SurfCornerFeature objects then F.u is a 2xN matrix
% with each column the corresponding u coordinate.
%
% Notes::
% - SurfCornerFeature is a reference object.
% - SurfCornerFeature objects can be used in vectors and arrays
%
% Reference::
% "SURF: Speeded Up Robust Features",
% Herbert Bay, Andreas Ess, Tinne Tuytelaars, Luc Van Gool,
% Computer Vision and Image Understanding (CVIU),
% Vol. 110, No. 3, pp. 346--359, 2008
%
% See also ISURF, PointFeature, ScalePointFeature, OrientedScalePointFeature, SiftPointFeature.

classdef SurfPointFeature < OrientedScalePointFeature

    properties
        %image_id_
    end % properties

    methods
        function f = SurfPointFeature(varargin)
        %SurfPointFeature.SurfPointFeature Create a SURF point feature object
        %
        % F = SurfPointFeature() is a point feature object with null parameters.
        %
        % F = SurfPointFeature(U, V) is a point feature object with specified
        % coordinates.
        %
        % F = SurfPointFeature(U, V, STRENGTH) as above but with specified strength.
        %
        % F = SurfScalePointFeature(U, V, STRENGTH, SCALE) as above but with specified
        % feature scale.
        %
        % F = SurfPointFeature(U, V, STRENGTH, SCALE, THETA) as above but with specified
        % feature orientation.
        %
        % See also isurf, OrientedScalePointFeature.

            f = f@OrientedScalePointFeature(varargin{:});  % invoke the superclass constructor
        end

        function val = image_id(features)
            val = [features.image_id_];
        end

        function [m,corresp] = match(f1, f2, varargin)
        %SurfPointFeature.match Match SURF point features
        %
        % M = F.match(F2, OPTIONS) is a vector of FeatureMatch objects that
        % describe candidate matches between the two vectors of SURF
        % features F and F2.  Correspondence is based on descriptor
        % similarity.
        %
        % [M,C] = F.match(F2, OPTIONS) as above but returns a correspodence
        % matrix where each row contains the indices of corresponding features
        % in F and F2  respectively.
        %
        % Options::
        % 'thresh',T    Match threshold
        % 'top',N       Take strongest N matches
        %
        % Notes::
        % - to obtain all matches use 'top', Inf
        %
        % See also FeatureMatch.

            if isempty(f2)
                m = [];
                corresp = [];
                return;
            end

            % HACK, pg 463 requires median (default)
            opt.thresh = [];
            opt.top = [];
            opt.all = false;
            opt = tb_optparse(opt, varargin);

            % Put the landmark descriptors in a matrix
            D1 = f1.descriptor;
            D2 = f2.descriptor;

            % Find the best matches

%             err=zeros(1,length(f1));
%             cor1=1:length(f1);
%             cor2=zeros(1,length(f1));
%             for i=1:length(f1)
%                 distance = sum((D2-repmat(D1(:,i),[1 length(f2)])).^2,1);
%                 [err(i),cor2(i)] = min(distance);
%             end

            % vectorized code (much faster)
            cor1 = 1:length(f1);
            [cor2,err] = closest(D1, D2);
            err = err.^2;  % closest returns distance, old code used distance squared

            % Sort matches on vector distance
            [err, ind] = sort(err);
            cor1=cor1(ind);
            cor2=cor2(ind);

            % Build a list of FeatureMatch objects
            m = [];
            cor = [];
            for i=1:length(f1)
                k1 = cor1(i);
                k2 = cor2(i);
                mm = FeatureMatch(f1(k1), f2(k2), err(i));
                m = [m mm];
                cor(:,i) = [k1 k2]';
            end

            % find the strongest matches
            if ~opt.all
                if ~isempty(opt.top)
                    % take the N strongest

                    k = min(opt.top, numcols(cor));
                    cor(:,k+1:end) = [];
                    m(k+1:end) = [];
                else
                    % use a threshold

                    if isempty(opt.thresh)
                        % use median of errors if no threshold given
                        thresh = median(err);
                    else
                        thresh = opt.thresh;
                    end

                    k = err > thresh;
                    cor(:,k) = [];
                    m(k) = [];
                end
            end

            if nargout > 1
                corresp = cor;
            end
        end

    end % methods

    methods(Static)

        % the MEX functions live in a private subdirectory, so these static methods
        % provide convenient access to them

        function Ipts = surf(im, opt)

            if exist('detectSURFFeatures') && 0
                % Use the CVST version
                fprintf('Using CVST\n');

                % do some option translation
                points = detectSURFFeatures(im);

                if ~isinf(opt.nfeat)
                    % choose strongest
                    points = points.selectStrongest(opt.nfeat);
                end

                [features,valid_points] = extractFeatures(im, points);

                p = valid_points.Location;
                Ipts = struct(...
                    'x',           num2cell(p(:,1)), ...
                    'y',           num2cell(p(:,2)), ...
                    'scale',       num2cell(valid_points.Scale), ...
                    'strength',    num2cell(valid_points.Metric), ...
                    'orientation', num2cell(valid_points.Orientation), ...
                    'descriptor',  num2cell(features, 2)   )';

            elseif exist('OpenSurf')
                % Use OpenSurf (MATLAB)
                params.octaves = opt.octaves;   % for OpenSurf
                if ~isempty(opt.thresh)
                    params.tresh = opt.thresh;      % for OpenSurf, (sic)
                end

                Ipts = OpenSurf(im, params);

            if false
                % Use surfmex (mex OpenCV wrapper)

                %if exist('surfpoints') == 3
                fprintf('MEX\n');
                % do the OpenCV/MEX version
                % put the results into the same return format as OpenSurf
                params.extended = 0;
                params.nOctaves = opt.octaves;
                if ~isempty(opt.thresh)
                    params.hessianThreshold = opt.thresh;
                end

                try
                    [p,d,l,info] = surfpoints(iint(im), params);
                catch me
                    if strcmp(me.identifier, 'MATLAB:UndefinedFunction')
                        error('MVTB:SurfPointFeature:notinstalled', 'Contributed software for SURF features is not installed')
                    end
                end

                % returns
                % p    point coordinates, one per column
                % d    SURF descriptor, one per column
                % l    sign of the Laplacian (light or dark feature)
                % info other parameters, per row: scale, strength,
                %      orientation

                % put the data into a vector of structs format to
                % match OpenSurf
                Ipts = struct('x', num2cell(p(1,:)), ...
                              'y', num2cell(p(2,:)), ...
                              'scale', num2cell(info(1,:)), ...
                              'strength', num2cell(info(2,:)), ...
                              'orientation', num2cell(info(3,:)), ...
                              'descriptor', num2cell(d,1)   );

            end
        end
    end
    end
end % classdef

author: petercorke

repo: machinevision-toolbox-matlab


text:
function pde = LinearFun2HM(am,ap,bm,bp,kappa,r,x0,y0,z0)
%% USAGE: polynomial solution for Poisson equation
%  Last Modified: 02/21/2020 by Xu Zhang

%% PDE Structure
pde = struct('intf',@intf,...
    'exactu1',@exactu1,'exactu2',@exactu2,'exactu3',@exactu3,...
    'um1',@um1,'um2',@um2,'um3',@um3,'up1',@up1,'up2',@up2,'up3',@up3,...
    'Dxu',@Dxu,'Dxum',@Dxum,'Dxup',@Dxup,'Dyu',@Dyu,...
    'Dyum',@Dyum,'Dyup',@Dyup,'Dzu',@Dzu,'Dzum',@Dzum,'Dzup',@Dzup,...
    'f1',@f1,'f2',@f2,'f3',@f3,...
    'fm1',@fm1,'fm2',@fm2,'fm3',@fm3,...
    'fp1',@fp1,'fp2',@fp2,'fp3',@fp3,...
    'A',@A,'Am',@Am,'Ap',@Ap,'one',@one,...
    'B',@B,'Bm',@Bm,'Bp',@Bp);

pde.bm = bm;
pde.bp = bp;
pde.am = am;
pde.ap = ap;
pde.kappa = kappa;
%% interface function
    function u = intf(x,y,z)
        u = ((x-x0).^2 + (y-y0).^2 + (z-z0).^2).^(1/2)/r-1;
    end

%% exact solution
    function u = exactu1(x,y,z)
        u = um1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up1(x(id),y(id),z(id));
    end
    function u = exactu2(x,y,z)
        u = um2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up2(x(id),y(id),z(id));
    end
    function u = exactu3(x,y,z)
        u = um3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = up3(x(id),y(id),z(id));
    end
    function u = um1(x,y,z)
        u = 2*(x-x0)/bm;
    end
    function u = um2(x,y,z)
        u = 2*(y-y0)/bm;
    end
    function u = um3(x,y,z)
        u = 2*(z-z0)/bm;
    end
    function u = up1(x,y,z)
        u = 2*(x-x0)/bp;
    end
    function u = up2(x,y,z)
        u = 2*(y-y0)/bp;
    end
    function u = up3(x,y,z)
        u = 2*(z-z0)/bp;
    end

%% Boundary Function
    function u = gD1(x,y,z)
        u = exactu1(x,y,z);
    end
    function u = gD2(x,y,z)
        u = exactu2(x,y,z);
    end
    function u = gD3(x,y,z)
        u = exactu3(x,y,z);
    end
%% Derivative of the exact solution
    function u = Dxu(x,y,z)
        u = Dxum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dxup(x(id),y(id),z(id));
    end
    function u = Dyu(x,y,z)
        u = Dyum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dyup(x(id),y(id),z(id));
    end
    function u = Dzu(x,y,z)
        u = Dzum(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Dzup(x(id),y(id),z(id));
    end
    function u = Dxum(x,y,z)
        u = zeros(size(x));
    end
    function u = Dyum(x,y,z)
        u = zeros(size(x));
    end
    function u = Dzum(x,y,z)
        u = zeros(size(x));
    end
    function u = Dxup(x,y,z)
        u = zeros(size(x));
    end
    function u = Dyup(x,y,z)
        u = zeros(size(x));
    end
    function u = Dzup(x,y,z)
        u = zeros(size(x));
    end

    function u = Duker(x,y,z)
        u = zeros(size(x));
    end

%% right hand side function
    function u = f1(x,y,z)
        u = fm1(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp1(x(id),y(id),z(id));
    end
    function u = f2(x,y,z)
        u = fm2(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp2(x(id),y(id),z(id));
    end
    function u = f3(x,y,z)
        u = fm3(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = fp3(x(id),y(id),z(id));
    end

    function u = fm1(x,y,z)
        u = -2*(x-x0)*kappa;
    end
    function u = fm2(x,y,z)
        u = -2*(y-y0)*kappa;
    end
    function u = fm3(x,y,z)
        u = -2*(z-z0)*kappa;
    end
    function u = fp1(x,y,z)
        u = -2*(x-x0)*kappa;
    end
    function u = fp2(x,y,z)
        u = -2*(y-y0)*kappa;
    end
    function u = fp3(x,y,z)
        u = -2*(z-z0)*kappa;
    end

%% Diffusion coefficient function
    function u = A(x,y,z)
        u = Am(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Ap(x(id),y(id),z(id));
    end
    function u = Am(x,y,z)
        u = am*ones(size(x));
    end
    function u = Ap(x,y,z)
        u = ap*ones(size(x));
    end
%% Mass coefficient function
    function u = B(x,y,z)
        u = Bm(x,y,z);
        id = intf(x,y,z) > 0;
        u(id) = Bp(x(id),y(id),z(id));
    end
    function u = Bm(x,y,z)
        u = bm*ones(size(x));
    end
    function u = Bp(x,y,z)
        u = bp*ones(size(x));
    end

%% Other function
    function u = one(x,y,z)
        u = ones(size(x));
    end
end

repo: ifem


text:
% TIMEDEP This code solves the nonlinear parabolic pde
%
% u_t = u_xx + exp(u); u(0,t) = u(1,t) = 0; u(x,0) = 0; 0 < t < 1.
%
% with the backward Euler discretization. Newton's method is used
% for the nonlinear solver. The Jacobian is tridiagonal, so we
% use the banded differencing function.
%
% The value of u at the current time and the time step are passed
% to the nonlinear residual as MATLAB global variables.
%
% This problem is 1-D, so we can store the time history of the
% integration and draw a surface plot.
%
global uold dt
dt=.1;
nx=63; nt=1+1/dt;
dx=1/(nx+1);
tval=0:dt:1;
xval=0:dx:1;
%
% Use tight tolerances, Newton's method, and a tridiagonal Jacobian.
%
tol=[1.d-6,1.d-6];
parms=[40, 1, 0, 1, 1, 1];
uhist=zeros(nx+2,nt);
uold=zeros(nx,1);
for it=1:nt-1
    [unew,it_hist,ierr]=nsold(uold,'ftime',tol,parms);
    uhist(2:nx+1,it+1)=unew;
    uold=unew;
end
%
% Plot the results.
%
mesh(tval,xval,uhist)

repo: Aero-matlab