text:
function [alpha, beta, stability] = evaluateAlphaBetaParam(process, noisy, dt)
% evaluateAlphaBetaParam - evaluates alpha and beta parameters for alpha-beta filter
% With this parameters, alpha-beta filter becomes a steady-state Kalman filter
%
% Syntax:          [xkp,vkp] = alphaBetaFilter(xm, dt, xk, vk, alpha, beta)
%   [alpha, beta, stability] = evaluateAlphaBetaParam(process, noisy, dt)
%     [alpha,beta,stability] = evalABGParam([alpha,beta])
%
% Inputs:
%   process - real system state
%     noisy - measured system state
%
% Outputs:
%   alpha - alpha parameter
%    beta - beta parameter
%
% Other m-files required: none
% Subfunctions: none
% MAT-files required: none
%
% See also: alphaBetaFilter;

% Author: Marco Borges, Ph.D. Student, Computer/Biomedical Engineer
% UFMG, PPGEE, Neurodinamica Lab, Brazil
% email address: marcoafborges@gmail.com
% Website: http://www.cpdee.ufmg.br/
% June 2013; Version: v2; Last revision: 2013-09-18
% Changelog:
%  v2 - add Stability test
%
%------------------------------- BEGIN CODE -------------------------------

if nargin == 3
    varProcess = var(process);
    varNoise = var(noisy);
    l = varProcess * dt / varNoise; % lambda
    r = (4+l-sqrt(8*l+l^2))/4;
    alpha = 1 - r^2;
    beta = 2*(2-alpha)-4*sqrt(1-alpha);
elseif nargin == 1 && length(process) == 2
    alpha = process(1);
    beta = process(2);
else
    error('evaluateAlphaBetaParam : Incorrect Parameters!');
end

if ( alpha > 0 && alpha < 2 && beta > 0 && beta < (4-2*alpha) )
    stability = 'Stable';
else
    stability = 'Unstable';
end

end
%-------------------------------- END CODE --------------------------------

author: Sable

repo: mcbench-benchmarks


function types = guitypes()
% GUITYPES returns a cell array of MATLAB/Octave GUI object types
%
% Syntax
%   types = guitypes()
%
% These types are ignored by matlab2tikz and figure2dot.
%
% See also: matlab2tikz, figure2dot

types = {'uitoolbar', 'uimenu', 'uicontextmenu', 'uitoggletool',...
         'uitogglesplittool', 'uipushtool', 'hgjavacomponent',...
         'matlab.graphics.shape.internal.Button', ...
         'matlab.ui.controls.ToolbarPushButton', ...
         'matlab.ui.controls.ToolbarStateButton', ...
         'matlab.ui.controls.ToolbarDropdown', ...
         'matlab.ui.controls.AxesToolbar', ...
         };
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a script that runs the motion compensation algorithm and
% computes the MSE and the MI error for all the dataSets in the workDir
% directory (can be manually changed).
%gb 05/15/05
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% temp file storage is on biac2 - change this below too!!!
% if isunix
%     networkPath = '/biac2/wandell';
% else
%     networkPath = '\\White\biac2-wandell\';
% end
% workDir = fullfile(networkPath,'data','reading_longitudinal_study','To_correct');
%%%% temp file storage is on teal - change this below too!!!
if isunix
    networkPath = '/teal/scr1';
else
    networkPath = '\\teal\scr1\';
end
workDir = fullfile(networkPath,'fmri','mcToCorrect');

cd(workDir);
[fileNum dirName] = countDirs(pwd);
for count = 3:fileNum
    currentDir = dirName{count};
    cd(fullfile(workDir,currentDir));
    transformAll;
	close all
	clear all
	mrVista
	vw = getSelectedInplane;
	computeError

    % Redefining variables because they are cleared in transform
    if isunix
        networkPath = '/teal/scr1';
    else
        networkPath = '\\teal\scr1\';
    end
    workDir = fullfile(networkPath,'fmri','mcToCorrect');
    cd(workDir);
    [fileNum dirName] = countDirs(pwd);
    pause(0.001);
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% Plot P

a = getElement(logsout, 'Kdata');
t = squeeze(a.Values.time);
K = squeeze(a.Values.data);

plot(t, K(1,:), '-', ...
	 t, K(2,:), '-.',...
     t, K(3,:), ':', ...
     'LineWidth', 2)

legend('K_1', 'K_2', 'K_3')

xlabel('Time (sec)')

%%
print('Ch2_ex1_fig3_K', '-depsc')

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
classdef MimEditDelete < MimGuiPlugin
    % MimEditDelete. Gui Plugin for setting paint colour to transparent
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimEditDelete is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Delete'
        SelectedText = 'Delete'
        ToolTip = 'Paint over edits'
        Category = 'Segmentation label'
        Visibility = 'Dataset'
        Mode = 'Edit'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 5
        ButtonHeight = 1
        Icon = 'paint.png'
        IconColour = GemMarkerPoint.DefaultColours{8}
        Location = 38
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.PaintBrushColour = 0;
        end

        function enabled = IsEnabled(gui_app)
             enabled = gui_app.IsDatasetLoaded && gui_app.ImagePanel.OverlayImage.ImageExists && ...
                isequal(gui_app.ImagePanel.SelectedControl, 'Paint');
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.PaintBrushColour == 0;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function y = cvx_isaffine( x, full )
error( nargchk( 1, 2, nargin ) );
if nargin == 1,
    y = true;
else
    y = true( size( x ) );
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:
% Initialize a test suite using test functions in the body of a function
%
% This function is not intended to be used directly; instead, it should be
% used at the top of a '.m'-file where the subfunctions may contain test
%
% To include subfunctions as tests and return a test-suite based on them,
% use the following layout:
%
%     function test_suite=my_test_of_abs
%         try % assignment of 'localfunctions' is necessary in Matlab >= 2016
%             test_functions=localfunctions();
%         catch % no problem; early Matlab versions can use initTestSuite fine
%         end
%         initTestSuite;
%
%     function test_abs_scalar
%         assertTrue(abs(-1)==1)
%         assertEqual(abs(-NaN),NaN);
%         assertEqual(abs(-Inf),Inf);
%         assertEqual(abs(0),0)
%         assertElementsAlmostEqual(abs(-1e-13),0)
%
%     function test_abs_vector
%         assertEqual(abs([-1 1 -3]),[1 1 3]);
%
%     function test_abs_exceptions
%         if moxunit_util_platform_is_octave()
%             assertExceptionThrown(@()abs(struct),'');
%         else
%             assertExceptionThrown(@()abs(struct),...
%                                    'MATLAB:UndefinedFunction');
%         end
%
% After suite=my_test_of_abs(), suite will be a test suite testing
% the test_abs_scalar and test_abs_vector subfunctions
%
% Each test function should either start with 'test_' or end with '_test',
% and have no output arguments
%
% It is absolutely necessary that the main function has one output, and
% that the output variable is named 'test_suite'
%
% NNO Jan 2014

    get_example_syntax=@()sprintf('%s\n',...
                            '',...
                            '    function test_suite=test_foo',...
                            '    % tests for function ''foo''',...
                            ['        try % assignment of '...
                                    '''localfunctions'' is necessary '...
                                    'in Matlab >= 2016'],...
                            ['            test_functions='...
                                        'localfunctions();'],...
                            ['        catch % no problem; early Matlab '...
                                                'versions can use '...
                                                'initTestSuite fine'],...
                            '        end',...
                            '        initTestSuite;',...
                            '',...
                            '    function test_foo1',...
                            '        % your test code here',...
                            '',...
                            '    function test_foo2',...
                            '        % your test code here'...
                            );

    e=lasterror();
    [call_stack, idx] = dbstack('-completenames');
    if numel(call_stack)==1
        example_syntax=get_example_syntax();

        error(['The script ''%s'' must be called from within '...
                    'another function, typically using the '...
                    'following syntax:\n\n%s'...
                    ],...
                    call_stack(1).name,...
                    example_syntax);
    end

    caller_fn=call_stack(idx+1).file;
    sub_func_struct=moxunit_util_mfile_subfunctions(caller_fn);

    has_error=false;

    if moxunit_util_platform_supports('localfunctions_in_script')
        n_sub_func=numel(sub_func_struct);

        handle_candidates=cell(n_sub_func,1);
        for k=1:n_sub_func
            handle_candidates{k}=str2func(sub_func_struct(k).name);
        end

    else
        if ~exist('test_functions','var')
            func_name=call_stack(2).name;
            example_syntax=get_example_syntax();

            func_error=@()error(...
                    ['In %s:\nThe variable ''test_functions'' is '...
                    'not assigned and this version\nof Matlab does '...
                    'not support the assignment of function handles\n'...
                    'through calls by functions. \nAdapt the function '...
                    '%s so that it resembles the following '...
                    'syntax:\n%s'],...
                    func_name, func_name, example_syntax);
            handle_candidates={func_error};
            has_error=true;
        else
            handle_candidates=test_functions;
        end
    end

    test_re=moxunit_util_get_test_name_regexp();

    n_sub_func=numel(sub_func_struct);
    sub_func_names=arrayfun(@(i)sub_func_struct(i).name,1:n_sub_func,...
                            'UniformOutput',false);

    n_handles=numel(handle_candidates);
    test_suite=MOxUnitTestSuite();

    for k=1:n_handles
        handle=handle_candidates{k};

        add_test=false;

        if has_error
            name=call_stack(1).name;
            add_test=true;
        else
            name=func2str(handle);
            if moxunit_util_regexp_matches(name,test_re)
                idx=find(strcmp(name,sub_func_names));
                if ~isempty(idx)
                    assert(numel(idx)==1, 'name match should be unique');
                    add_test=sub_func_struct(idx).nargout==0;
                end
            end
        end

        if add_test
            test_case=MOxUnitFunctionHandleTestCase(...
                                                name,...
                                                caller_fn, handle);
            test_suite=addTest(test_suite, test_case);
        end

    end

    % If not called from another function, execute the test directly
    % and remove the test_suite variable
    if numel(call_stack)==2
        disp(run(test_suite));

        % Avoid showing "ans = MOxUnitTestSuite object: 1-by-1"
        % when run without explicitly assigning output to a variable
        clear test_suite;
    end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function seed = hammersley_seed_get ( )

%*****************************************************************************80
%
%% HAMMERSLEY_SEED_GET gets the seed vector of the leaped Hammersley subsequence.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer SEED(1:DIM_NUM), the Hammersley seed vector.
%
  global hammersley_BASE
  global hammersley_LEAP
  global hammersley_DIM_NUM
  global hammersley_SEED
  global hammersley_STEP

  dim_num = hammersley_DIM_NUM;

  if ( ~halham_dim_num_check ( dim_num ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HAMMERSLEY_SEED_GET - Fatal error!\n' );
    error ( 'HAMMERSLEY_SEED_GET - Fatal error!' );
  end

  if ( length ( hammersley_SEED ) ~= dim_num )

    hammersley_SEED = [];

    for i = 1 : dim_num
      hammersley_SEED(i) = 0;
    end

  end

  seed(1:dim_num) = hammersley_SEED(1:dim_num);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function g = p14_g ( n, x )

%*****************************************************************************80
%
%% P14_G evaluates the gradient for problem 14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the values of the variables.
%
%    Output, real G(N), the gradient of the objective function.
%
  g = zeros ( n, 1 );

  for j = 1 : n

    if ( mod ( j, 2 ) == 1 )
      g(j) = - 2.0 * ( 1.0 - x(j) );
    else
      g(j) = 200.0 * ( x(j) - x(j-1) * x(j-1) );
      g(j-1) = g(j-1) - 400.0 * x(j-1) * ( x(j) - x(j-1) * x(j-1) );
    end

  end

  return
end


text:
% Fig. 5.28  Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
% script for LEFT side of Figure 5.28.  Use fig5_28b to see RIGHT side.

clf
n=[1 2];
d=conv([1 1 0],[1 13]);
nc=91*[1 .05];
dc=[1 .01];
nol=conv([0 0 n],nc);
dol=conv(d,dc);
 rlocus(nol,dol);
 hold on
 title('Fig.5.28a Root locus for lead plus lag')
 axis([-20 4 -9 9])
 z=0:.1:.9;
 wn=2:2:19;
 sgrid(z, wn)
 dcl=nol+dol;
 r=roots(dcl);
 plot(r,'*')
 hold off


text:
function mask = pm_mask(angvar,mthres,ndil)
% Creating a mask that will determine how far
% we should proceed with phase unwrapping.
% FORMAT: mask = pm_mask(angvar,mthrea,ndil)
%
% Input:
% angvar     : Map of variance of angle estimate.
% mthres     : Threshold for variance beyond which
%              phase unwrapping is considered too
%              uncertain. Default value (pi^2)/6
%              is half the variance of a U[-pi,pi]
%              distribution.
% ndil       : We can optionally specify a no. of
%              erodes-dilates to apply to the mask
%              in order to exclude areas connected
%              only by thin bridges to the rest of
%              the brain.
%
% Output:
% mask       : Well...
%
% Jesper Andersson 26/9-03.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Jesper Andersson
% $Id: pm_mask.m 1317 2008-04-08 16:16:38Z chloe $

if nargin < 2
   mthres = (pi^2)/6;
end
if nargin < 3
   ndil = 0;
end

%
% Threshold map of angular variance and pick
% largest connected component.
%
mask = double(angvar < mthres);
[lmask,num] = spm_bwlabel(mask,6);
n = histc(lmask(:),[0:num]+0.5);
[mv,mi] = max(n);
indx = lmask(:)==mi;
mask(~indx) = 0;

if ndil
   dmask = mask;
   for i=1:ndil
      dmask = spm_erode(dmask);
   end
   [lmask,num] = spm_bwlabel(dmask,6);
   n = histc(lmask(:),[0:num]+0.5);
   [mv,mi] = max(n);
   indx = lmask(:)==mi;
   dmask(~indx) = 0;
   for i=1:ndil
      dmask = spm_dilate(dmask);
   end
   mask = mask.*dmask;
end

return

author: spm


text:
function [simple_rules,skipped] = simplifyGrRules(grRules,expanded)
%simplifyGrRules  Simplify and condense the logic of model grRules.
%
% simplifyGrRules uses Matlab's symbolic toolbox to mathematically simplify
% gene rules. Each rule is converted into a symbolic boolean expression,
% and these expressions are simplified yielding rules that are as reduced/
% simplified as possible, while still maintaining functional equivalency to
% the original grRule logic.
%
% The function will also remove unnecessary parentheses, trailing ANDs/ORs,
% and duplicate genes in model grRules. The function requires that the gene
% names/IDs do not contain spaces, parentheses, "&", "|", or "#".
% Furthermore, grRules must separate gene names/IDs from ANDs and ORs by at
% least one space. For example,
%
%       OK: "GENE1 and GENE2", or "GENE1 & GENE2"
%   NOT OK: "GENE1andGENE2", or "GENE1&GENE2"
%
%
% *!WARNING!* Simplifying grRules with this function may eliminate some
%             gene-reaction associations.
%             For example: 'G1 or (G1 and G2)' simplifies to 'G1'
%
%
% USAGE:
%
%   [simple_rules,skipped] = simplifyGrRules(grRules,expanded);
%
%
% INPUT:
%
%   grRules        The grRules field from a genome-scale metabolic model.
%                  Note that grRules can use written ("and" "or") or
%                  symbolic ("&" "|") boolean operators.
%
%   expanded       (Optional, default FALSE). If TRUE, the grRule will be
%                  returned in its expanded form, where groups of ANDs
%                  (enzyme complexes) are separated by ORs, rather than
%                  using a shorter nested format.
%
%                  For example:
%
%                   input: G1 and G2 and (G3 or G4) and (G4 or G3 or G3) and G1
%
%              simplified: G1 and G2 and (G3 or G4)
%
%       simplify + expand: (G1 and G2 and G3) or (G1 and G2 and G4)
%
%                  NOTE: Expanding rules can make them very long, and in
%                        some cases will exceed Matlab's limit on the
%                        length of a symbolic expression, requiring them to
%                        be skipped.
%
%
% OUTPUT:
%
%   simple_rules   Updated/simplified grRules.
%
%   skipped        Logical vector indicating which grRules were skipped
%                  because they were too long to be converted into a
%                  symbolic equation. For these cases, the grRule will be
%                  copied into simple_rules without any simplification.
%

if nargin < 2
    expanded = false;
end

% perform a preliminary "clean" of the gene rules
grRules = cleanGrRules(grRules);

% check if the grRules use written or symbolic boolean operators
if any(contains(grRules,{'&','|'}))
    Btype = 'symbol';  % record rule type, to revert back at the end
else
    % convert "and" to "&" and "or" to "|"
    grRules = regexprep(grRules, ' or ', ' | ');
    grRules = regexprep(grRules, ' and ', ' & ');
    Btype = 'text';  % record rule type, to revert back at the end
end

% get list of all gene IDs present in grRules
rxnGenes = cellfun(@(r) unique(regexp(r,'[^&|\(\) ]+','match')),grRules,'UniformOutput',false);
nonEmpty = ~cellfun(@isempty,rxnGenes);
allGenes = unique([rxnGenes{nonEmpty}]');

% Replace all gene IDs with G1, G2, G3, etc. to avoid potential problems
% with certain gene ID formats (e.g., Entrez IDs are numbers, which will
% cause problems).
tempIDs = strcat('G',arrayfun(@num2str,(1:length(allGenes))','UniformOutput',false));  % construct new ID vector
convertGene = @(g) tempIDs{ismember(allGenes,g)};  % define conversion function
grRules = regexprep(grRules, '[^&|\(\) ]+', '${convertGene($0)}');  % find and convert all gene IDs to new temporary type

% identify all rules containing both ANDs and ORs
AndOrInd = find(contains(grRules,'&') & contains(grRules,'|'));

% initialize output specifying which rules were skipped due to complexity
skipped = false(size(grRules));

% iterate through complex rules (those containing both ANDs and ORs)
for i = 1:length(AndOrInd)

    r = grRules{AndOrInd(i)};

    % attempt to convert string to symbolic equation
    try
        reqn = str2sym(r);
    catch
        fprintf('grRule #%u was too long (%u characters), and therefore skipped.\n',AndOrInd(i),numel(r));
        skipped(AndOrInd(i)) = true;
        continue
    end

    % simplify eqn, allowing for a max of 10 steps
    reqn = simplify(reqn,10);

    % attempt to expand equation if requested
    reqn_noexpand = reqn;
    if ( expanded )
        try
            reqn = expand(reqn);
        catch
            fprintf('Failed to expand grRule #%u due to excess length. Rule will remain in simplified form.\n',AndOrInd(i));
        end
    end

    % Add parentheses around groups of genes. This may add some
    % unnecessary parentheses, but they will be cleaned up later.
    if is_or(reqn) || is_and(reqn)
        % the input and output of this function are of the type "char"
        try
            % sometimes an expanded equation will fail here
            r = add_parentheses(char(reqn));
        catch
            fprintf('Failed to expand grRule #%u due to excess length. Rule will remain in simplified form.\n',AndOrInd(i));
            reqn = reqn_noexpand;
            r = add_parentheses(char(reqn));
        end
        r = regexprep(r,'^\(|\)$','');  % remove the extra set of parentheses that always enclose the expression
    end

    % update rule in grRule vector
    grRules{AndOrInd(i)} = r;

end

% restore original gene IDs
restoreGene = @(g) allGenes{ismember(tempIDs,g)};  % define conversion function
grRules = regexprep(grRules, '[^&|\(\) ]+', '${restoreGene($0)}');  % find and restore all gene IDs to original ID

% change boolean operators back to original type
if strcmpi(Btype,'text')
    grRules = regexprep(grRules, ' \| ', ' or ');
    grRules = regexprep(grRules, ' & ', ' and ');
end

% assign output
simple_rules = grRules;

end

%% Add parentheses around groups of equation terms
function eqn_str_new = add_parentheses(eqn_str)
% This is a recursive function, which groups terms in an equation with
% parentheses in order to clearly indicate the order of operations. This is
% necessary to avoid ambiguity when the equations are converted back into
% gene rules.
%
% NOTE: both eqn_str_new and eqn_str are strings (char).

% obtain symbolic form of equation
eqn = str2sym(eqn_str);

% separate equation into terms/pieces ("children")
eqn_pieces = arrayfun(@char,children(eqn),'UniformOutput',false);

% find pieces that are further separable, and recursively separate until
% all pieces are single genes
ind = find(contains(eqn_pieces,{'|','&'}));
for i = 1:length(ind)
    eqn_pieces{ind(i)} = add_parentheses(eqn_pieces{ind(i)});
end

% re-assemble equation pieces, adding | or &, depending on the input eqn
if is_or(eqn)
    eqn_str_new = char(strcat('(',join(eqn_pieces,' | '),')'));
elseif is_and(eqn)
    eqn_str_new = char(strcat('(',join(eqn_pieces,' & '),')'));
else
    % if the input eqn was simply one gene, just return the input
    eqn_str_new = eqn_str;
end

end

%% Functions to test whether an equation is a collection of ORs or ANDs
function res = is_or(eqn)
% check if outermost operations contain ORs
res = length(regexp(char(eqn),'\|')) > length(regexp(char(children(eqn)),'\|'));
end

function res = is_and(eqn)
% check if outermost operations contain ANDs
res = length(regexp(char(eqn),'&')) > length(regexp(char(children(eqn)),'&'));
end

author: SysBioChalmers

repo: Human-GEM


text:
function [err,time,solver,eqn] = femStokesHdiv(node,elem,pde,bdFlag,option,varargin)
%% FEMPOISSON solve Poisson equation by various finite element methods
% Created by Ming Wang, at Nov., 2012.
%

%% Check input arguments
if ~exist('node','var') || ~exist('elem','var')
    [node,elem] = squaremesh([0,1,0,1],0.25);  % default mesh is a square
end
if ~exist('option','var'), option = []; end
if ~exist('pde','var')
   pde = StokesZulehnerdata;                   % default data
end
if ~exist('bdFlag','var')
    bdFlag = setboundary(node,elem,'Dirichlet');
end
% default elemType
if ~isfield(option,'elemType')
    option.elemType = 'RT0-P0';
end

%% Parameters
option = femoption(option);
elemType = option.elemType;
maxIt = option.maxIt;
maxN = option.maxN;
L0 = option.L0;
refType = option.refType;

%% Generate an initial mesh
for k = 1:L0
    if strcmp(refType,'red')
        [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
    elseif strcmp(refType,'bisect')
        [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
    end
end

%% Initialize err
erruL2 = zeros(maxIt,1); erruIuhH1 = zeros(maxIt,1); erruInf = zeros(maxIt,1);
errpL2 = zeros(maxIt,1); errpIphL2 = zeros(maxIt,1); errpInf = zeros(maxIt,1);
errpL2re = zeros(maxIt,1);
errwL2 = zeros(maxIt,1); errwIwh = zeros(maxIt,1);
errTime = zeros(maxIt,1); solverTime = zeros(maxIt,1);
assembleTime = zeros(maxIt,1); meshTime = zeros(maxIt,1);
itStep = zeros(maxIt,1);  stopErr = zeros(maxIt,1); flag = zeros(maxIt,1);
N = zeros(maxIt,1);

%% Finite Element Method
for k = 1:maxIt
    % solve the equation
    switch elemType
        case 'RT0-P0'   % RT0-P0 mixed FEM
            [u,p,w,edge,eqn,info] = StokesRT0(node,elem,bdFlag,pde,option);
        case 'BDM1B-P0' % (BDM1+bubble)-P0 mixed FEM
            [u,p,w,edge,eqn,info] = StokesBDM1B(node,elem,bdFlag,pde,option);
    end
    % compute error
    tic;
    % ================== error for velocity ==================
    if isfield(pde,'exactu')
        switch elemType
            case 'RT0-P0'
                erruL2(k) = getL2errorRT0(node,elem,pde.exactu,u);
                % interpolation
                uI = faceinterpolate(pde.exactu,node,edge,'RT0');
                ufreeDof = eqn.ufreeDof;
                u0 = u(ufreeDof);
                uI0 = uI(ufreeDof);
                erruIuhH1(k) = sqrt((u0-uI0)'*eqn.A*(u0-uI0));
                erruInf(k) = max(abs(u0-uI0));
            case 'BDM1B-P0'
                erruL2(k) = getL2errorBDM1(node,elem,pde.exactu,u);
                % interpolation
                uI = u;
                uI(1:2*size(edge,1)) = faceinterpolate(pde.exactu,node,edge,'BDM1');
                ufreeDof = eqn.ufreeDof;
                u0 = u(ufreeDof);
                uI0 = uI(ufreeDof);
                erruIuhH1(k) = sqrt((u0-uI0)'*eqn.A*(u0-uI0));
                erruInf(k) = max(abs(u0-uI0));
        end
    end
    % ================== error for pressure ==================
    if isfield(pde,'exactp')
        area = simplexvolume(node,elem);
        errpL2(k) = getL2error(node,elem,pde.exactp,p);
        pI = Lagrangeinterpolate(pde.exactp,node,elem,'P0');
        rp = recoverP02P1(node,elem,p,'LA');
        % --------------------------
        errpL2(k) = getL2error(node,elem,pde.exactp,p);
        errpL2re(k) = getL2error(node,elem,pde.exactp,rp);
        errpIphL2(k) = sqrt(dot((pI-p).^2,area));
        errpInf(k) = max(abs(p-pI));
    end
    % ================== error for vorticity ==================
    if isfield(pde,'exactw')
        switch elemType
            case 'RT0-P0'
                errwL2(k) = getL2error(node,elem,pde.exactw,w);
                wI = Lagrangeinterpolate(pde.exactw,node,elem,'P1');
                errwIwh(k) = sqrt((w-wI)'*eqn.Mv*(w-wI));
            case 'BDM1B-P0'
                errwL2(k) = getL2error(node,elem,pde.exactw,w);
                wI = Lagrangeinterpolate(pde.exactw,node,elem,'P2',edge);
                errwIwh(k) = sqrt((w-wI)'*eqn.Mv*(w-wI));
        end
    end
    % ================== record ==================
    errTime(k) = toc;
    % record time
    solverTime(k) = info.solverTime;
    assembleTime(k) = info.assembleTime;
    % record solver information
    itStep(k) = info.itStep;
    stopErr(k) = info.stopErr;
    flag(k) = info.flag;
    % plot
    N(k) = size(node,1);
    if option.plotflag && N(k) < 3e3 % show mesh and solution for small size
       figure(1); showresult(node,elem,p); pause(0.1);
       figure(2);
       showsolutionRT(node,elem,u); pause(0.1);
    end
    if N(k) > maxN
        break;
    end
    % refine mesh
    tic;
    if strcmp(refType,'red')
        [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
    elseif strcmp(refType,'bisect')
        [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
    end
    meshTime(k) = toc;
end

%% Plot convergence rates
if option.rateflag
    figure;
    set(gcf,'Units','normal');
    set(gcf,'Position',[0.25,0.25,0.80,0.40]);
    subplot(1,3,1)
    h = 1./sqrt(N(1:k));
    showrateh3(h,erruL2(1:k),2,'k-+','||u-u_h||',...
               h,erruIuhH1(1:k),2,'r-*','||u_I-u_h||_1',...
               h,erruInf(1:k),2,'b-*','||u_I-u_h||_{\infty}');
    subplot(1,3,2)
    showrateh4(h,errpL2(1:k),2,'k-+', '||p - p_h||',...
               h,errpL2re(1:k),2,'g-+','||p- p_h^r||',...
               h,errpIphL2(1:k),2,'r-+','||p_I - p_h||',...
               h,errpInf(1:k),2,'b-+','||p_I - p_h||_{\infty}');
    subplot(1,3,3)
    showrateh2(h,errwL2(1:k),2,'k-+','||w - w_h||',...
               h,errwIwh(1:k),2,'r-+','||w_I - w_h||');
end

% Output
err = struct('N',N,'uL2',erruL2(1:k),'uInf',erruInf,'uIuhH1',erruIuhH1(1:k),...
             'pL2',errpL2(1:k),'pIphL2',errpIphL2(1:k),'pInf',errpInf(1:k),'pL2re',errpL2re(1:k),...
             'wL2',errwL2(1:k),'wIwhL2',errwIwh(1:k));
time = struct('N',N,'err',errTime(1:k),'solver',solverTime(1:k), ...
              'assmble',assembleTime(1:k),'mesh',meshTime(1:k));
solver = struct('N',N(1:k),'itStep',itStep(1:k),'time',solverTime(1:k),...
                'stopErr',stopErr(1:k),'flag',flag(1:k));

%% Display error on screen
ts = zeros(k,3); ts = char(ts);
% error of u
fprintf('===========================================================\n');
display('   N   ||u_I-u_h||_1  ||u-u_h||   ||u_I-u_h||_{max}');
display([num2str(err.N) ts num2str(err.uIuhH1,'%0.5e') ts num2str(err.uL2,'%0.5e') ts num2str(err.uInf,'%0.5e')]);
% error of p
fprintf('===========================================================\n');
display('   N   ||p_I-p_h||   ||p-p_h||    ||p_I-p_h||_{max}      ||p_I - p^r_h||');
display([num2str(err.N) ts num2str(err.pIphL2,'%0.5e') ts num2str(err.pL2,'%0.5e') ts num2str(err.pInf,'%0.5e') ts num2str(err.pL2re,'%0.5e')]);
% error of w
fprintf('===========================================================\n');
display('   N   ||w_I-w_h||   ||w-w_h||');
display([num2str(err.N) ts num2str(err.wIwhL2,'%0.5e') ts num2str(err.wL2,'%0.5e')]);
fprintf('===========================================================\n');

repo: Aero-matlab


text:
% m-file which generates data from a simple system with a localized
% nonlinearity and processes it using RLM's zeroing algorithm.
%
% Matt Allen, June 27, 2005
% MSA: Major Update - Nov, 2005 - changed to two arbitrary subsystems.
%
clear all; close all;

% Set up system parameters
% m1 = 1; m2 = 1; m3 = 1; m4 = 1; m5 = 1;
    M1 = eye(5);
    M2 = 0.1*eye(2); % Subsystem is about 10% the weight of normal system.

    kall = 1e6;
    k1 = kall; k2 = kall; k3 = kall; k4 = kall; k5 = kall; k6 = kall; kat = kall;
    K1 = [k1, -k1, 0, 0, 0; % System 1
        -k1, k1+k2, -k2, 0, 0
        0, -k2, k2+k3, -k3, 0;
        0, 0, -k3, k3+k4, -k4;
        0, 0, 0, -k4, k4];
    Ns1 = 5; Ns2 = 2; Ntot = 7;
    K2 = [k5, -k5; % System 2 - attached subsystem - dof's
        -k5, k5];
    Mtot = [M1, zeros(size(M1,1),size(M2,2));
        zeros(size(M2,1),size(M1,2)), M2];
    Ktot = [K1, zeros(size(K1,1),size(K2,2));
        zeros(size(K2,1),size(K1,2)), K2];
    % Specify how systems are connected
    at_ns = [2,6]; % attachment happens between this pair of nodes.
    fext_ns = 5; % Node at which external force is applied
    fnl_vec = zeros(Ntot,1); fnl_vec(at_ns(1)) = -1; fnl_vec(at_ns(2)) = 1;
    fext_vec = zeros(Ntot,1); fext_vec(fext_ns(1)) = 1;
    dnodes = [1:Ntot];
    vnodes = [(Ntot+1):(Ntot*2)];

    % Damping Matrix
    cfactk = 0.00003; % multiplied by K to get damping.
    cfactm = 8; % Multiplied by M to get damping
    % for k = 1:5; eval(['c',num2str(k),' = cfact*k',num2str(k),';']); end
    C1 = cfactk*K1 + cfactm*M1;
    C2 = cfactk*K2 + cfactm*M2; % proportional damping
    Ctot = [C1, zeros(size(C1,1),size(C2,2));
        zeros(size(C2,1),size(C1,2)), C2];

    % Linearized System Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % add a linear spring between attachment nodes
        Mlin = Mtot;
        Klin = Ktot;
        Klin(at_ns,at_ns) = Klin(at_ns,at_ns) + kat*[1, -1; -1, 1];
        % Linear Damping matrix:
        Clin = Ctot;
        Clin(at_ns,at_ns) = Clin(at_ns,at_ns) + cfactk*kat*[1, -1; -1, 1];

    % State Space Eigenanalysis
    Slin = [Clin, Mlin; Mlin, zeros(size(Mlin))];
    Rlin = [-Klin, zeros(size(Mlin));
        zeros(size(Mlin)), Mlin];
    Alin = (Slin\Rlin);
    [Philin,lamlin] = eig(Alin);
    lamlin = diag(lamlin);
    [junk,sind] = sort(abs(lamlin) - 0.001*min(abs(lamlin))*(imag(lamlin) > 0));
    lamlin = lamlin(sind);
    Philin = Philin(:,sind);

%         % Plot Mode Shapes
%         figure(1)
%         hls = plot([1:Ntot], imag(Philin(1:Ntot,1:2:end ))); grid on;
%         legend(hls, num2str([1:Ntot].'));
%         xlabel('X-coordinate');
%         ylabel('Im\{Mode Shape\}');

    wns = abs(lamlin);
    fns = wns/2/pi;
    zts = -real(lamlin)./abs(lamlin);
    disp('Natural Frequencies:, Damping Ratios:');
    [fns, zts]
%     DispFnZt(lamlin) - replace 2 lines abouve with this if you have the EMA Functions toolbox

    % Nonlinear Parameters
    NLType = 'bang'
    if strcmp(NLType,'bang');
        % Bang (Contact) Nonlinearity
        delcont = 1e-3;
        k4mult = 20; % Factor by which k4 increases: k4_contact = k4*k4mult
        c4mult = 1; % Factor by which c4 increases: c4_contact = c4*c4mult
    elseif strcmp(NLType,'cubic');
        % Cubic Spring
        katnl = 1e8;
    else
        error('NLType not recognized');
    end
    % Force paramters
        % length of half-sine force pulse
        % This is normalized in the EOM to unit area and multiplied by Afnl
        tfp = 1e-4;
        Afnl = 4e9;

    % Which Response to Use in evaluating Nonlinearity (i.e. x1, x2, x3..?)
    respind = 6;

    % Number of numerical derivatives to evaluate.  This M-file simulates
    % the displacement response of the 5-DOF system.  To simulate the
    % measurement of the velocity or acceleration response, the
    % displacement response is differentiated using finite differences
    % (i.e. Matlab's 'diff' command.)  This parameter sets the number of
    % derivatives to perform:
    nders = 2;
        % nders = 0; => use displacement
        % nders = 1; => use velocity
        % nders = 2; => use acceleration

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Put all of this in a structure for the EOM
global eom
S = whos;
for k = 1:length(S)
    eval(['eom.',S(k).name,' = ',S(k).name,';']);
end

    % Compute TF of Linearized System
    F = zeros(Ntot,1); F(fext_ns) = 1;
    wlin = [1:1:500]*2*pi;
    Halin = zeros(Ntot,length(wlin));
    for k = 1:length(wlin)
        Halin(:,k) = wlin(k).^2*([Mlin*-wlin(k).^2 + Clin*i*wlin(k) + Klin]\F);
    end
    Halin = Halin.';

%     figure(2)
%     semilogy(wlin/2/pi,abs(Halin)); grid on;
%     xlabel('Frequency (Hz)')
%     ylabel('H_{ACCEL}(\omega)');
%     for k = 1:Ntot; leg_txt{k} = ['H_{',num2str(k),',5}']; end
%     legend(leg_txt);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nonlinear Simulation
% Compute sample rate and time
    elas_modes = find(imag(lamlin) > 0);
    dt = (5*2*max(fns(elas_modes)))^-1;
    Tmin = log(0.001)/(-min(abs(real(lamlin(elas_modes)))));
    N = 2^ceil(log2(Tmin/dt));

    ts = [0:dt:(N-1)*dt].';

% Specify initial time step so the solver doesn't miss the pulse:
    odeopts = odeset('InitialStep',tfp/10,'AbsTol',1e-9);

    % [ts_ode,hnl] = ode45(@nldetect_anex2_v1_eom,ts,zeros(10,1),odeopts,eom);

% Alternate Calling form:
sol = ode45(@nldt_bng_cub_v2_eom,[ts(1),ts(end)],zeros(Ntot*2,1),odeopts);
    [hnl, hnldot] = deval(sol,ts);
    hnl = hnl.';  hnldot = hnldot.';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Linear Simulation
ltisys = ss(Alin,Slin\[F; zeros(Ntot,1)],[eye(Ntot), zeros(Ntot)], zeros(Ntot,1));

    % Use fine sample increment over duration of pulse
    if tfp < dt
        tlini = [0:100].'*tfp/100;
        ulini = (Afnl*2*tfp/pi)*sin((pi/tfp)*tlini);
        % Over Pulse duation to first sample
        [ylini,junki,xlini] = lsim(ltisys, [ulini], [tlini], zeros(Ntot*2,1));
        [ylin, junk2, xlin2] = lsim(ltisys, zeros(2,1), [junki(end); ts(2)], xlini(end,:).');
        % Free Response
        [ylin, junk3, xlin3] = lsim(ltisys, zeros(size(ts(2:end))), ts(2:end), xlin2(end,:).');
        hlin = [xlini(1,:); xlin2(end,:); xlin3(2:end,:)];
    else
        error('tfp > dt not yet supported!');
    end

% Plot Results
figure(10); set(gcf,'Position',[532   189   560   567]);
subplot(2,1,1)
plot(ts, hlin(:,respind), ts, hnl(:,respind), 'o-',sol.x, sol.y(respind,:), '.:'); grid on;
legend('h_{(r,5)} LIN','h_{(r,5)} NL','h_{(r,5)} ODE45');
xlabel('time (s)');
title('Comparison of Linear and Nonlinear Time Responses')

% Approximate derivatives of the response:
disp(['Number of Derivatives Requested: ',num2str(nders)]);
ndhlin = [zeros(nders,1); dt^-nders*diff(hlin(:,respind ),nders,1);];
ndhnl = [zeros(nders,1); dt^-nders*diff(hnl(:,respind ),nders,1); ];
%   Derivatives computed by "diff" don't appear to be ideal.  I don't know
%   that they are equivalent to a forward difference scheme.
% Select ODE45 derivatives - USE ODE45 Derivatives in the analysis!
if nders == 0;
    dhnlode = hnl(:,respind);
elseif nders == 1;
    dhnlode = hnl(:,respind + Ntot);
elseif nders == 2;
    dhnlode = hnldot(:,respind + Ntot);
end
tsdh = ts;%ts((1+nders):end);
subplot(2,1,2)
plot(tsdh, ndhlin, tsdh, ndhnl, ts, dhnlode, '.'); grid on;
legend(['h_D_{(',num2str(respind),',5)} LIN'],['h_D_{(',num2str(respind),',5)} NL'],...
    ['h_D_{(',num2str(respind),',5)} NL-ODE45']);
xlabel('time (s)');
title('Comparison of Derivatives of Linear and Nonlinear Responses')

% Find Actual force curve:
% This must be modified manually
% Fnl2 = Mtot(at_ns(2),at_ns(2))*hnldot(:,Ntot+at_ns(2)) - k5*(hnl(:,at_ns(2)) - hnl(:,at_ns(2)+1));
Fnl2 = Mtot(at_ns(2),:)*(hnldot(:,vnodes).') + Ktot(at_ns(2),:)*(hnl(:,dnodes).');
    if abs(fext_vec(at_ns(1))) > 0; warning('Equation for Force may not be correct'); end
Fnl3 = - Mtot(at_ns(1),:)*(hnldot(:,vnodes).') - Ktot(at_ns(1),:)*(hnl(:,dnodes).');
    if abs(fext_vec(at_ns(1))) > 0; warning('Equation for Force may not be correct'); end
delt_nl = hnl(:,at_ns(1)) - hnl(:,at_ns(2));

% Track the magnitude of the displacement of spring 4 - for NL
figure(11); set(gcf, 'Units','Normalized'); set(gcf,'Position',[0.069444,  0.48444, 0.43819, 0.39111]);
    set(gcf,'Name','Force vs. Delta Curve','Toolbar','figure');
subplot(1,3,1)
deltas = [-100:100]*(max(delt_nl)/100);
if strcmp(NLType,'bang');
    for k = 1:length(deltas)
        if deltas(k) < eom.delcont % No Contact Occurs
            Fnl(k) = eom.kat*deltas(k);
        else                % Contact Occurs
            Fnl(k) = eom.kat*eom.k4mult*deltas(k);
        end
    end
elseif strcmp(NLType,'cubic');
    Fnl = kat*deltas.*(1+katnl.*deltas.^2);
end
plot(Fnl, deltas*1e3, ...
    k4*deltas, deltas*1e3, '--',...
    Fnl2,delt_nl*1e3,'.',...
    Fnl3,delt_nl*1e3,'k.'); grid on;
ax1h = gca;
ylabel('\bf\Delta - NL Spring (mm)'); xlabel('\bfSpring Force')
title('\bfSpring Curve');
hx = subplot(1,3,2:3)
plot(ts*1e3, delt_nl*1e3); grid on;
% Set Axes to Equal
    ylim2 = get(gca,'Ylim');
    set(ax1h,'Ylim',ylim2);
xlabel('\bftime (ms)'); ylabel('\bf\Delta - NL Spring (mm)');
title('\bfDisplacement of NL Spring');

% Add button to make y-axes equal
set(hx,'tag','Main');
    PointsPerPixel = 72/get(0,'ScreenPixelsPerInch');
    hf = gcf;
    h1 = uicontrol('Parent',hf, ...
        'Units','points', ...
        'Position',[168, 51, 45, 25]*PointsPerPixel, ...
        'Callback','set_equal_ylims(hf);', ...
        'String','EqYlim', ...
        'Tag','YAxEqual');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Analysis

if strcmp(NLType,'bang');
    % ZNLDetect(hlin(:,respind),ts,[0,0.2],400)
    [Hmat,fs,ts_zc] = ZNLDetect(dhnlode,tsdh,[0,0.1],500);
    % ZNLDetectC(hnldot(:,6:10),tsdh,[0,0.3],30,400)
elseif strcmp(NLType,'cubic');
    [Hmat,fs,ts_zc] = ZNLDetect(dhnlode,tsdh,[0,0.15],500,4);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Alternate analysis - use all sensors with ZNLDetectC
% This gives 20 ZEFFTs equally distributed between 0 and 0.1 s, 0 to 500 Hz
%
% [Hmat,fs,ts_zc] = ZNLDetectC(hnldot(:,8:end),tsdh,[0,0.1],500,20);


text:
function a = moler3 ( m, n )

%*****************************************************************************80
%
%% MOLER3 returns the MOLER3 matrix.
%
%  Formula:
%
%    if ( I == J )
%      A(I,J) = I
%    else
%      A(I,J) = min(I,J) - 2
%
%  Example:
%
%    N = 5
%
%     1 -1 -1 -1 -1
%    -1  2  0  0  0
%    -1  0  3  1  1
%    -1  0  1  4  2
%    -1  0  1  2  5
%
%  Properties:
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A is positive definite.
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A has a simple Cholesky factorization.
%
%    A has one small eigenvalue.
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of A.
%
%    Output, real A(M,N), the matrix.
%
  for i = 1 : m
    for j = 1 : n
      if ( i == j )
        a(i,j) = i;
      else
        a(i,j) = min ( i, j ) - 2;
      end
    end
  end

  return
end


text:
function closewindow(name)
%CLOSEWINDOW Close a window.
%   CLOSEWINDOW(NAME) closese the window with a certain NAME.

% Matthew J. Simoneau, June 2003
% Copyright 2003 The MathWorks, Inc.

if ~ispc
    error('This only works on Windows.');
end

if ~libisloaded(mfilename);
    loadlibrary('user32.dll',@userproto,'alias',mfilename);
end

h = calllib(mfilename,'FindWindowA',[],name);
calllib(mfilename,'CloseWindow',h);

%===============================================================================
function [fcns,structs,enuminfo] = userproto

fcns=[]; structs=[]; enuminfo=[]; fcns.alias={};

%  HWND _stdcall FindWindowA(LPCSTR,LPCSTR);
fcns.name{1} = 'FindWindowA';
fcns.calltype{1} = 'stdcall';
fcns.LHS{1} = 'voidRef';
fcns.RHS{1} = {'int8Ref', 'string'};

fcns.name{2} = 'CloseWindow';
fcns.calltype{2} = 'stdcall';
fcns.LHS{2} = 'int32';
fcns.RHS{2} = {'voidRef'};


text:
function [MSLPoints,coeffData]=ellips2MSLHelmert(points,useNGAApprox,modelType,coeffData)
%%ELLIPS2MSLHELMERT  Given points in WGS-84 ellipsoidal coordinates,
%                    convert the ellipsoidal height components of the
%                    points into heights above mean-sea level (MSL)
%                    using Helmert's projection method with the EGM2008
%                    gravitational (and terrain correction) models. This
%                    differs from the more precise but seldom-used
%                    Pizzetti's projection method in that the curvature of
%                    the plumb line is completely ignored.
%
%INPUTS: points One or more points given in WGS-84 geodetic latitude and
%               longitude, in radians, and height, in meters for which the
%               corresponding MSL heights are desired. To convert N points,
%               points is a 3XN matrix with each column having the format
%               [latitude;longitude; height].
%  useNGAApprox If one wishes the intermediate tide-free geoid computation
%               to match the results of the National Geospatial
%               Intelligence Agency's code to three digits after the
%               decimal point, then this should be true. Setting this to
%               false might produce results that are marginally more
%               accurate. The default is false if this parameter is
%               omitted.
%     modelType An optional parameter specifying coefficient model to
%               load if coeffData is not provided. Possible values are
%               0 (The default if omitted) Load the EGM2008 model.
%               1 Load the EGM96 model.
%     coeffData A set of pre-loaded coefficients that can speed up the
%               computation by eliminating the need to compute them on the
%               fly. coeffData is as defined for the coeffData return value
%               of getEGMGeoidHeight.
%
%OUTPUTS: MSLPoints A 3XN array of the converted points, where each vector
%                   contains [latitude;longitude;MSLHeight]; The latitude
%                   and longitude are unchanged from the input data, but
%                   the heights have been converted from WGS-84 ellipsoidal
%                   heights to MSL heights.
%         coeffData The coeffData coefficients that can be passed to
%                   another call of ellips2MSLHelmert or MSL2EllipseHelmert
%                   to make it faster.
%
%Note that this function will be very slow if one hasn't called
%CompileCLibraries to compile the spherical harmonic synthesis functions.
%
%As described in Chapter 5.5 of [1], the MSL height using Helmert's
%projection is just the ellipsoidal height minus the geoid height. This
%function calls getEGMGeoidHeight to get the geoidal height and
%subtracts it off.
%
%Helmert's projection calculated the height using a line from the point to
%the reference ellipsoid, subtracting off the distance from the point to
%the geoid along the line. Pizzetti's projection follows the curved plumb
%line down to the geoid. Thus, the latitude and longitude of the projected
%point on the geoid are not the same as that of the point being projected.
%The MSL height obtained is also slightly different. However, the
%difference is small and Pizzetti's projection is difficult to use, so
%Helmert's projection is almost always used.
%
%The inverse of this function is MSL2EllipseHelmert.
%
%REFERENCES:
%[1] B. Hofmann-Wellenhof and H. Moritz, Physical Geodesy, 2nd ed.
%    SpringerWienNewYork, 2006.
%
%January 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<2||isempty(useNGAApprox))
    useNGAApprox=false;
end

if(nargin<3||isempty(modelType))
    modelType=0;
end

if(nargin>3)
    [geoidHeight,coeffData]=getEGMGeoidHeight(points(1:2,:),1,useNGAApprox,modelType,coeffData);
else
    [geoidHeight,coeffData]=getEGMGeoidHeight(points(1:2,:),1,useNGAApprox,modelType);
end

MSLHeight=points(3,:)-geoidHeight';
numPoints=size(points,2);

MSLPoints=zeros(3,numPoints);
MSLPoints(1:2,:)=points(1:2,:);
MSLPoints(3,:)=MSLHeight;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function pass = test_toFileIVP(pref)
%TEST_TOFILEIVP   Test exporting all IVP demos to an m-file.
%
% This test only checks whether nothing breaks, it does not try to solve the
% problems.

% Find the folders which demos are stored in. The chebguiDemos folder lives in
% the trunk folder, find the path of the Chebfun trunk.
trunkPath = chebfunroot();

% Append directory information
bvppath = fullfile(trunkPath, 'chebguiDemos', 'ivpdemos');

% Set up file to exporting to
tempPath = fullfile(trunkPath, 'tests', 'chebgui');
tempFileName = 'tempExporterTest.m';

% Obtain directory information:
D = dir(bvppath);

% Create a BVP exporter object:
exporter = chebguiExporter.constructor('ivp');

% Wrap the test in a try-catch statement to ensure we don't files hanging around
% in the test folder:
try

    % Loop through the demos and export
    for demoCounter = 3:length(D)

        % The path to the .guifile
        demoPath = fullfile(bvppath, D(demoCounter).name);

        % Load demo
        cg = chebgui.demo2chebgui(demoPath);

        % Export the demo!
        toFile(exporter, cg, tempFileName, tempPath)
    end

    % Got here without crashing == success!
    pass = 1;

catch
    % Something went wrong == fail :(
    pass = 0;
end

% Delete the temporary file we wrote to:
delete(fullfile(tempPath, tempFileName))

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [ seed, xd, td ] = triangulation_order3_sample ( node_num, node_xy, ...
  triangle_num, triangle_node, num_ran, seed )

%*****************************************************************************80
%
%% TRIANGULATION_ORDER3_SAMPLE returns random points in a triangulation.
%
%  Discussion:
%
%    It is assumed that the triangulation consists of a set of non-overlapping
%    triangles.
%
%    The point is chosen uniformly in the area covered by the triangulation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the node coordinates.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, integer TRIANGLE_NODE(3,TRIANGLE_NUM), the nodes that make up the triangles.
%
%    Input, integer NUM_RAN, the number of points to sample.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, integer SEED, a seed for the random number generator.
%
%    Output, real XD(2,NUM_RAN), the sample points.
%
%    Output, integer TD(NUM_RAN), the triangle to which each sample point
%    belongs.
%
  dim_num = 2;
%
%  Compute the areas of the triangles.
%  Build a cumulative area vector.
%  Convert it to a relative cumulative area vector.
%
  area_cum(0+1) = 0.0;

  for i = 1 : triangle_num

    i1 = triangle_node(1,i);
    i2 = triangle_node(2,i);
    i3 = triangle_node(3,i);

    t(1:dim_num,1) = node_xy(1:dim_num,i1);
    t(1:dim_num,2) = node_xy(1:dim_num,i2);
    t(1:dim_num,3) = node_xy(1:dim_num,i3);

    area = triangle_area_2d ( t );

    area_cum(i+1) = area_cum(i-1+1) + area;

  end

  area_total = area_cum(triangle_num+1);

  area_cum(0+1:triangle_num+1) = area_cum(0+1:triangle_num+1) / area_total;
%
%  Pick random values.  A random value R indicates the corresponding triangle
%  whose cumulative relative area contains R.
%
%  Bracket the random value in the cumulative relative areas,
%  indicating a triangle.
%
%  Pick a random point in the triangle.
%
  for i = 1 : num_ran

    [ r, seed ] = r8_uniform_01 ( seed );

    [ left, right ] = r8vec_bracket ( triangle_num+1, area_cum, r );

    td(i) = right - 1;

    i1 = triangle_node(1,td(i));
    i2 = triangle_node(2,td(i));
    i3 = triangle_node(3,td(i));

    t(1:dim_num,1) = node_xy(1:dim_num,i1);
    t(1:dim_num,2) = node_xy(1:dim_num,i2);
    t(1:dim_num,3) = node_xy(1:dim_num,i3);

    [ xd(1:dim_num,i), seed ] = triangle_sample ( t, 1, seed );

  end

  return
end


text:
function varargout = BallDetectionGUI(varargin)
% BALLDETECTIONGUI MATLAB code for BallDetectionGUI.fig
%      BALLDETECTIONGUI, by itself, creates a new BALLDETECTIONGUI or raises the existing
%      singleton*.
%
%      H = BALLDETECTIONGUI returns the handle to a new BALLDETECTIONGUI or the handle to
%      the existing singleton*.
%
%      BALLDETECTIONGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in BALLDETECTIONGUI.M with the given input arguments.
%
%      BALLDETECTIONGUI('Property','Value',...) creates a new BALLDETECTIONGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before BallDetectionGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to BallDetectionGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help BallDetectionGUI

% Last Modified by GUIDE v2.5 30-Mar-2011 17:45:48

% Copyright 2011 The MathWorks, Inc.

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
   'gui_Singleton',  gui_Singleton, ...
   'gui_OpeningFcn', @BallDetectionGUI_OpeningFcn, ...
   'gui_OutputFcn',  @BallDetectionGUI_OutputFcn, ...
   'gui_LayoutFcn',  [] , ...
   'gui_Callback',   []);
if nargin && ischar(varargin{1})
   gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
   [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
   gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before BallDetectionGUI is made visible.
function BallDetectionGUI_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to BallDetectionGUI (see VARARGIN)

% Choose default command line output for BallDetectionGUI
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes BallDetectionGUI wait for user response (see UIRESUME)
% uiwait(handles.figure1);

movegui(handles.figure1, 'center');

warning('off', 'imaq:peekdata:tooManyFramesRequested');
imaqreset;

% --- Outputs from this function are returned to the command line.
function varargout = BallDetectionGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in StartStopBtn.
function StartStopBtn_Callback(hObject, eventdata, handles)
% hObject    handle to StartStopBtn (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

if strcmp(get(hObject, 'String'), 'Start')
   % create video object
   handles.vid = webcam;

   % set timer function
   handles.vid.TimerPeriod = 0.1;
   handles.vid.TimerFcn = {@videoTimerFunction, handles};

   % start acquisition
   start(handles.vid);

   % change button label
   set(hObject, 'String', 'Stop');

   % Update handles structure
   guidata(hObject, handles);
else
   % stop image acquisition
   stop(handles.vid);

   % delete video object
   delete(handles.vid);

   % clear image
   cla(handles.axes1);

   % change button label
   set(hObject, 'String', 'Start');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Timer Function
function videoTimerFunction(vid, eventdata, handles)

% get a single frame
pic = peekdata(vid,1);

% find ball
thresh = str2double(get(handles.threshEdit, 'String'));
h = get(handles.imageTypeSelection, 'SelectedObject');
imageType = get(h, 'UserData');
findBallFcn(pic, thresh, imageType, handles.axes1);

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

function threshEdit_Callback(hObject, eventdata, handles)
% hObject    handle to threshEdit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of threshEdit as text
%        str2double(get(hObject,'String')) returns contents of threshEdit as a double

% --- Executes during object creation, after setting all properties.
function threshEdit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to threshEdit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor','white');
end

% --- Executes when selected object is changed in imageTypeSelection.
function imageTypeSelection_SelectionChangeFcn(hObject, eventdata, handles)
% hObject    handle to the selected object in imageTypeSelection
% eventdata  structure with the following fields (see UIBUTTONGROUP)
%	EventName: string 'SelectionChanged' (read only)
%	OldValue: handle of the previously selected object or empty if none was selected
%	NewValue: handle of the currently selected object
% handles    structure with handles and user data (see GUIDATA)

% --- Executes during object deletion, before destroying properties.
function figure1_DeleteFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

warning('on', 'imaq:peekdata:tooManyFramesRequested');
imaqreset;


text:
function FEBioPath=getFEBioPath

% function FEBioPath=getFEBioPath
% ------------------------------------------------------------------------
% This function reads a user specified FEBio path name from a config text
% file. The text file may contain multiple paths (e.g. for multiple
% operational systems), however the first valid patch is always used.
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2015/10/05 Updated for GIBBON
% 2015/10/05 Enabled searching valid paths in multi-line config. file.
% 2017/06/10 changed empty output to an empty character array.
%------------------------------------------------------------------------

%%

filePath=mfilename('fullpath');
toolboxPath=fileparts(fileparts(filePath));
configPath=fullfile(toolboxPath,'config');

fileName=fullfile(configPath,'FEBioPath.txt');

%Import text file containing paths
try
    [T]=txtfile2cell(fileName);
catch
    T={};
end

%Get first valid path
FEBioPath='';
if ~isempty(T)
    for q=1:1:numel(T)
        pathName=T{q};
        if exist(pathName,'file')==2
            FEBioPath=pathName;
            return
        end
    end
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
%% Swan - Wiki
% The student's guide to clean code development
% Task 7: UML of FEM code in Swan repository

% Instructions: run the following code, selecting previously the 'Swan'
% main folder as your current matlab path

file = 'test2d_triangle';
a.fileName = file;
s = FemDataContainer(a);
fem = FEM.create(s);
fem.solve();

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
% This file contains a list of metabolites required by the AGORA microbial
% reconstructions to grow.

AGORAessential ={'EX_12dgr180[u]'
'EX_26dap_M[u]'
'EX_2obut[u]'
'EX_3mop[u]'
'EX_4hbz[u]'
'EX_ac[u]'
'EX_acgam[u]'
'EX_acmana[u]'
'EX_ade[u]'
'EX_adn[u]'
'EX_adocbl[u]'
'EX_ala_D[u]'
'EX_ala_L[u]'
'EX_arab_D[u]'
'EX_arg_L[u]'
'EX_asn_L[u]'
'EX_btn[u]'
'EX_ca2[u]'
'EX_cbl1[u]'
'EX_cgly[u]'
'EX_chor[u]'
'EX_chsterol[u]'
'EX_cit[u]'
'EX_cl[u]'
'EX_cobalt2[u]'
'EX_csn[u]'
'EX_cu2[u]'
'EX_cys_L[u]'
'EX_dad_2[u]'
'EX_dcyt[u]'
'EX_ddca[u]'
'EX_dgsn[u]'
'EX_fe2[u]'
'EX_fe3[u]'
'EX_fol[u]'
'EX_for[u]'
'EX_gal[u]'
'EX_glc_D[u]'
'EX_glu_L[u]'
'EX_gly[u]'
'EX_glyc3p[u]'
'EX_glyc[u]'
'EX_gthox[u]'
'EX_gua[u]'
'EX_h2s[u]'
'EX_h[u]'
'EX_his_L[u]'
'EX_hxan[u]'
'EX_ile_L[u]'
'EX_k[u]'
'EX_lanost[u]'
'EX_leu_L[u]'
'EX_lys_L[u]'
'EX_met_L[u]'
'EX_mg2[u]'
'EX_mn2[u]'
'EX_mqn7[u]'
'EX_mqn8[u]'
'EX_nac[u]'
'EX_ncam[u]'
'EX_nmn[u]'
'EX_ocdca[u]'
'EX_ocdcea[u]'
'EX_orn[u]'
'EX_phe_L[u]'
'EX_pheme[u]'
'EX_pi[u]'
'EX_pnto_R[u]'
'EX_pro_L[u]'
'EX_ptrc[u]'
'EX_pydx[u]'
'EX_pydxn[u]'
'EX_q8[u]'
'EX_rib_D[u]'
'EX_ribflv[u]'
'EX_ser_L[u]'
'EX_sheme[u]'
'EX_so4[u]'
'EX_spmd[u]'
'EX_thm[u]'
'EX_thr_L[u]'
'EX_thymd[u]'
'EX_trp_L[u]'
'EX_ttdca[u]'
'EX_tyr_L[u]'
'EX_ura[u]'
'EX_val_L[u]'
'EX_xan[u]'
'EX_xyl_D[u]'
'EX_zn2[u]'
% new
'EX_arg_L[u]'
'EX_ca2[u]'
'EX_cgly[u]'
'EX_cl[u]'
'EX_cobalt2[u]'
'EX_cu2[u]'
'EX_fe2[u]'
'EX_fol[u]'
'EX_k[u]'
'EX_mg2[u]'
'EX_mn2[u]'
'EX_mqn7[u]'
'EX_ocdca[u]'
'EX_pheme[u]'
'EX_ribflv[u]'
'EX_so4[u]'
'EX_spmd[u]'
'EX_thm[u]'
'EX_trp_L[u]'
'EX_zn2[u]'
'EX_pnto_R[u]'
'EX_sheme[u]'
'EX_thymd[u]'
'EX_pydx[u]'
'EX_26dap_M[u]'
'EX_ala_L[u]'
'EX_fe3[u]'
'EX_pi[u]'
'EX_ttdca[u]'
'EX_mqn8[u]'
'EX_q8[u]'
'EX_2dmmq8[u]'
'EX_btn[u]'
'EX_nmn[u]'
'EX_asn_L[u]'
'EX_adpcbl[u]'
'EX_hxan[u]'
'EX_gthrd[u]'
'EX_4hbz[u]'
'EX_12dgr180[u]'
'EX_nac[u]'
'EX_gthox[u]'
'EX_ptrc[u]'
'EX_ura[u]'
'EX_ddca[u]'
'EX_lys_L[u]'
'EX_ala_D[u]'
'EX_arab_D[u]'
'EX_4abz[u]'
'EX_ile_L[u]'
'EX_met_L[u]'
'EX_ocdcea[u]'
'EX_tyr_L[u]'
'EX_val_L[u]'
'EX_his_L[u]'
'EX_leu_L[u]'
'EX_pro_L[u]'
'EX_cit[u]'
'EX_phe_L[u]'
'EX_thr_L[u]'
'EX_adn[u]'
'EX_3mop[u]'
'EX_cys_L[u]'
'EX_gal[u]'
'EX_ser_L[u]'
'EX_orn[u]'
'EX_h2s[u]'
'EX_xyl_D[u]'
'EX_2obut[u]'
'EX_lanost[u]'
'EX_rib_D[u]'
'EX_dcyt[u]'
'EX_glu_L[u]'
'EX_ac[u]'
'EX_ncam[u]'
'EX_cytd[u]'
'EX_amet[u]'
'EX_no2[u]'
'EX_cbl1[u]'
'EX_gly[u]'
'EX_gln_L[u]'
'EX_acgam[u]'
'EX_dad_2[u]'
'EX_ade[u]'
'EX_gua[u]'
'EX_for[u]'
'EX_h[u]'
'EX_pydxn[u]'
'EX_dgsn[u]'
'EX_malt[u]'
'EX_fald[u]'
'EX_amp[u]'
'EX_chor[u]'
'EX_glyc[u]'
'EX_glyc3p[u]'
'EX_h2o[u]'
'EX_chsterol[u]'
'EX_acmana[u]'
'EX_xan[u]'
'EX_gsn[u]'};

AGORAessential = unique(AGORAessential);

author: opencobra

repo: cobratoolbox


text:
function [dt6,b0] = dtiDeformer(dt6,b0,deformField)
% function [outDt6,outb0] = dtiDeformer(inputImgFN,deformFieldFN,dfType,outputImgFN)
% DTI Deformer: Warps dt6 image according to given deformation field
%
% ARGUMENTS:
% dt6: Input dt6 image
% b0: input b0 image
% deformField: Deformation field (XxYxZx3)
%
% HISTORY:
% 2004.01.19 GSM (gmulye@stanford.edu) wrote it and cleaned up the interface
%

dimSrc = size(dt6);
H = zeros(dimSrc(1),dimSrc(2),dimSrc(3),3);
for x = 1:dimSrc(1)
    for y = 1:dimSrc(2)
        for z = 1:dimSrc(3)
            H(x,y,z,:) = [y x z];
        end
    end
end
newH = H + deformField;
%Giant matrix of coordinates

% CREATE NEW DT6 IMAGE
coordsList = reshape(newH,dimSrc(1)*dimSrc(2)*dimSrc(3),3);
newDt6List = zeros(dimSrc(1)*dimSrc(2)*dimSrc(3),6);
for i = 1:6
    newDt6List(:,i) = myCinterp3(dt6(:,:,:,i),[dimSrc(1),dimSrc(2)],dimSrc(3),coordsList)';
end
newDt6 = zeros(dimSrc(1),dimSrc(2),dimSrc(3),6);
for i = 1:6
    newDt6(:,:,:,i) = reshape(newDt6List(:,i),dimSrc(1),dimSrc(2),dimSrc(3));
end
dt6 = newDt6; clear newDt6 newDt6List;

%CREATE NEW B0 IMAGE, SAVE OUT
newB0List(:,:,:) = myCinterp3(double(b0),[dimSrc(1),dimSrc(2)],dimSrc(3),coordsList)';
newB0 = reshape(newB0List,dimSrc(1),dimSrc(2),dimSrc(3));
b0 = int16(newB0); clear newB0List newB0; clear coordsList;

% %OUTPUT FILE NAMING CONVENTION
% if (~exist('outputImgFN', 'var'))
%     [a imgFN b c] = fileparts(inputImgFN);
%     us=findstr('_',inputImgFN);
%     subjectCode = inputImgFN(1:us(1)-1);
%     outputImgFN = ['registered_',subjectCode];
% end
% deformationNotes = [deformFieldFN];
% save(outputImgFN,'b0','xformToAcPc','xformToAnat','notes','deformationNotes',...
%     'mmPerVox','anat','dt6');

return


text:
function [ a, seed ] = r8but_random ( n, mu, seed )

%*****************************************************************************80
%
%% R8BUT_RANDOM randomizes a R8BUT matrix.
%
%  Discussion:
%
%    The R8BUT storage format is used for a banded upper triangular matrix.
%    The matrix is assumed to be zero above the MU-th superdiagonal.
%    The matrix is stored in an MU+1 by N array.
%    Columns are preserved.
%
%    The diagonal is stored in row MU+1 of the array.
%    The first superdiagonal in row MU, columns 2 through N.
%    The second superdiagonal in row MU-1, columns 3 through N.
%    The MU-th superdiagonal in row 1, columns MU+1 through N.
%
%  Example:
%
%    N = 5, MU = 2
%
%    A11 A12 A13   0   0
%      0 A22 A23 A24   0
%      0   0 A33 A34 A35
%      0   0   0 A44 A45
%      0   0   0   0 A55
%                --- ---
%                    ---
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of columns of the matrix.
%
%    Input, integer MU, the upper bandwidth.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real A(MU+1,N), the R8BUT matrix.
%
%    Output, integer SEED, an updated seed for the random number generator.
%
  for i = 1 : mu + 1

    for j = 1 : mu + 1 - i
      a(i,j) = 0.0;
    end

    for j = max ( 1, mu + 2 - i ) : n
      [ a(i,j), seed ] = r8_uniform_01 ( seed );
    end

  end

  return
end


text:
function complex_mat = getComplexesFromGeneRules(grRules)
%getComplexesFromGeneRules  Extract enzyme complex info from model grRules.
%
% getComplexesFromGeneRules uses the logic in grRules to identify which
% genes are treated as belonging to an enzyme complex, and returns a binary
% matrix describing all identified complexes and their gene constituents.
%
% USAGE:
%
%   complex_mat = getComplexesFromGeneRules(grRules);
%
% INPUT:
%
%   grRules      grRules from a genome-scale metabolic model structure.
%
% OUTPUT:
%
%   complex_mat  A (GxC) binary matrix, where G is the number of genes
%                appearing in grRules, and C is the number of complexes
%                idenditied. Each column of the matrix represents a unique
%                enzyme complex, with 1s indicating which genes belong to
%                the complex. Although the complexes are unique, they can
%                be subsets of one another.
%


text:
% Author:         Zhiwei (Tony) Qin (zq2107@columbia.edu)
% Date:           April, 2012
% Affiliation:    IEOR, Columbia University, New York

author: andrewssobral

repo: lrslibrary


text:
function SR = superresolve_bicubic(slidingWindow, magFactor)

    SR = imresize(slidingWindow.referenceFrame, magFactor, 'bicubic');

author: thomas-koehler

repo: SupER


text:
function pos =  DECODE_MFSK(Metrics,M)

%give the most likely output, given a set of input Metrics

max_metric = Metrics(1);
pos = 0;

for i = 2:M
    if Metrics(i) > max_metric
        max_metric = Metrics(i);
        pos = i - 1;
    end
end


text:
% given a mask of valid frames in sentences, return the an array of index of last
% valid frames in each sentence.
%
function last_idx = GetLastValidFrameIndex(mask)
T = size(mask,1);
delta = mask(2:end,:) - mask(1:end-1,:);
[max_delta, last_idx] = max(delta);
last_idx(max_delta==0) = T;

end

author: singaxiong

repo: SignalGraph


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
%   Unit step sequence

%  u[n]
 n1=-3:-1;
 n2=0:5;
 n=[n1 n2];
 u1=zeros(size(n1));
 u2=ones(size(n2));
 u=[u1 u2];
 stem(n,u)

%  second way
 figure
 n=-3:5
 n0=0;
 u=(n>=n0)
 stem(n,u)

 %u[n-n0]
 figure
 n1=-3:1;
 n2=2:5;
 n=[n1 n2];
 u1=zeros(size(n1));
 u2=ones(size(n2));
 u2=ones(size(n2));
 u=[u1 u2];
 stem(n,u)

%  second way
 figure
 n=-3:5
 n0=2;
 u=((n-n0)>=0)
 stem(n,u)

repo: mcbench-benchmarks


text:
function bounds = GetBounds(model, vel_lb, vel_ub, T, x0, xf)

    if nargin < 2
        vel_lb = [0.0,0];
        vel_ub = [0.0,0];
    end

    if nargin < 4
        T = 0.4;
    end

    if nargin < 5
        x0 = zeros(model.numState*2,1);
    end

    if nargin < 6
        xf = zeros(model.numState*2,1);
    end

    %% first get the model specific boundary values
    model_bounds = model.getLimits(); % x, dx, ddx, u

    model_bounds.options.enforceVirtualConstraints = true;

    % fixed joint constraint wrench
    model_bounds.inputs.ConstraintWrench.ffourBar.lb = -1000;
    model_bounds.inputs.ConstraintWrench.ffourBar.ub = 1000;

    % fixed joint constraints
    model_bounds.params.pfourBar.lb = zeros(4,1);
    model_bounds.params.pfourBar.ub = zeros(4,1);

    model_bounds.constrBounds.yaw_initial.lb = 0;
    model_bounds.constrBounds.yaw_initial.ub = 0;

    % feedback control gain for virtual constraints
    model_bounds.gains.kp = 100;
    model_bounds.gains.kd = 20;

    %%% Constraint bounds
    %     vel = [0.0, -0.8].';
    %     T   = 0.4;

    %%% Step Duration
    model_bounds.time.duration.lb = T;
    model_bounds.time.duration.ub = T;

    model_bounds.time.t0.lb = 0;
    model_bounds.time.t0.ub = 0;
    model_bounds.time.t0.x0 = 0;
    model_bounds.time.tf.lb = 0;
    model_bounds.time.tf.ub = T;
    model_bounds.time.tf.x0 = T;

    model_bounds.constrBounds.foot_clearance_1.lb = 0.03;
    model_bounds.constrBounds.foot_clearance_1.ub = 0.3;
    model_bounds.constrBounds.foot_clearance_2.lb = 0.15;
    model_bounds.constrBounds.foot_clearance_2.ub = 0.3;
    model_bounds.constrBounds.foot_clearance_3.lb = 0.05;
    model_bounds.constrBounds.foot_clearance_3.ub = 0.3;
    model_bounds.constrBounds.averageVelocity.lb = vel_lb;
    model_bounds.constrBounds.averageVelocity.ub = vel_ub;

    %     model_bounds.states.x.lb(1:2) = vel_lb;
    %     model_bounds.states.x.ub(1:2) = vel_ub;

    %     model_bounds.constrBounds.stepLength.lb = vel(1)*T;
    %     model_bounds.constrBounds.stepLength.ub = vel(1)*T;
    %
    %     model_bounds.constrBounds.stepWidth.lb = vel(1)*T;
    %     model_bounds.constrBounds.stepWidth.ub = vel(1)*T;

    model_bounds.constrBounds.footVelocityBeginning.lb = [-0.15, -0.15, 0]';
    model_bounds.constrBounds.footVelocityBeginning.ub = [0.15, 0.15, 0.3]';
    model_bounds.constrBounds.footVelocityEnd.lb = [-0.15, -0.15, -0.25]';
    model_bounds.constrBounds.footVelocityEnd.ub = [0.15, 0.15, -0.05]';

    %%% Common Virtual Constraints
    model_bounds.params.aoutput.lb = -2*pi;
    model_bounds.params.aoutput.ub = 2*pi;

    model_bounds.params.poutput.lb = [0, 0];
    model_bounds.params.poutput.ub = [0, T];
    model_bounds.params.poutput.x0 = [0, T];

    model_bounds.params.pRightFoot.lb = [-10, -10, 0, 0];
    model_bounds.params.pRightFoot.ub = [10, 10, 0, 0];

    model_bounds.params.pLeftFoot.lb = [-10, -10, 0, 0];
    model_bounds.params.pLeftFoot.ub = [10, 10, 0, 0];

    model_bounds.inputs.ConstraintWrench.fRightFoot.lb = [-1000,-1000,300,-1000]';
    model_bounds.inputs.ConstraintWrench.fRightFoot.ub = [1000,1000,1000,1000]';

    model_bounds.inputs.ConstraintWrench.fLeftFoot.lb = [-1000,-1000,300,-1000]';
    model_bounds.inputs.ConstraintWrench.fLeftFoot.ub = [1000,1000,1000,1000]';

    %% construct the boundary values for each domain
    bounds = struct();

    wt = 0.3662;
    vx_lb = vel_lb(1);
    vy_lb = vel_lb(2);
    vx_ub = vel_ub(1);
    vy_ub = vel_ub(2);

    %%
    bounds.RightStance1 = model_bounds;
    bounds.RightStance1.constrBounds.knee.lb = [deg2rad(50),deg2rad(50)];
    bounds.RightStance1.constrBounds.knee.ub = [deg2rad(60),deg2rad(150)];

    bounds.RightStance1.constrBounds.stepWidth.lb = -wt - vx_ub*T;
    bounds.RightStance1.constrBounds.stepWidth.ub = -wt - vx_lb*T;

    bounds.RightStance1.constrBounds.stepLength.lb = -vy_ub*T;
    bounds.RightStance1.constrBounds.stepLength.ub = -vy_lb*T;

    bounds.RightStance1.constrBounds.full_state.lb = x0;
    bounds.RightStance1.constrBounds.full_state.ub = x0;
    idx = [3:9, 12:14, 19:25, 28:30];
    bounds.RightStance1.constrBounds.actuated_state.lb = x0(idx);
    bounds.RightStance1.constrBounds.actuated_state.ub = x0(idx);
    %%
    bounds.LeftImpact.states.x = model_bounds.states.x;
    bounds.LeftImpact.states.xn = model_bounds.states.x;
    bounds.LeftImpact.states.dx = model_bounds.states.dx;
    bounds.LeftImpact.states.dxn = model_bounds.states.dx;

    bounds.LeftImpact.inputs = struct();
    bounds.LeftImpact.inputs.ConstraintWrench.ffourBar.lb = -20;
    bounds.LeftImpact.inputs.ConstraintWrench.ffourBar.ub = 20;
    bounds.LeftImpact.inputs.ConstraintWrench.fLeftFoot.lb = [-20,-20,0,-10]';
    bounds.LeftImpact.inputs.ConstraintWrench.fLeftFoot.ub = [20,20,150,10]';
    bounds.LeftImpact.params = struct();

    %% Right Stance
    bounds.LeftStance1 = model_bounds;

    bounds.LeftStance1.constrBounds.knee.lb = [deg2rad(50),deg2rad(50)];
    bounds.LeftStance1.constrBounds.knee.ub = [deg2rad(150),deg2rad(60)];

    bounds.LeftStance1.midState_QFit = zeros(6,3);
    bounds.LeftStance1.midState_dQFit = zeros(6,3);

    bounds.LeftStance1.constrBounds.stepWidth.lb = -wt + vx_lb*T;
    bounds.LeftStance1.constrBounds.stepWidth.ub = -wt + vx_ub*T;

    bounds.LeftStance1.constrBounds.stepLength.lb = vy_lb*T;
    bounds.LeftStance1.constrBounds.stepLength.ub = vy_ub*T;

    %%
    bounds.RightImpact.states.x = model_bounds.states.x;
    bounds.RightImpact.states.xn = model_bounds.states.x;
    bounds.RightImpact.states.dx = model_bounds.states.dx;
    bounds.RightImpact.states.dxn = model_bounds.states.dx;
    bounds.RightImpact.inputs = struct();
    bounds.RightImpact.inputs.ConstraintWrench.ffourBar.lb = -20;
    bounds.RightImpact.inputs.ConstraintWrench.ffourBar.ub = 20;
    bounds.RightImpact.inputs.ConstraintWrench.fRightFoot.lb = [-20,-20,0,-10]'';
    bounds.RightImpact.inputs.ConstraintWrench.fRightFoot.ub = [20,20,150,10]';
    bounds.RightImpact.params = struct();

    %%
    bounds.RightStance2 = bounds.RightStance1;

    bounds.RightStance2.midState_QFit = zeros(6,3);
    bounds.RightStance2.midState_dQFit = zeros(6,3);

    bounds.LeftStance2  = bounds.LeftStance1;

    bounds.LeftStance2.constrBounds.full_state.lb = xf;
    bounds.LeftStance2.constrBounds.full_state.ub = xf;

    bounds.LeftStance2.constrBounds.actuated_state.lb = xf(idx);
    bounds.LeftStance2.constrBounds.actuated_state.ub = xf(idx);
end

author: ayonga

repo: frost-dev

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


text:
function movePTPEllipse(tKuka,c,dir,ratio,theta,velocity,accel,TefTool)
%% This funciton is used to
% Move end-effector of the robot on an ellipse.

%% Syntax:
% movePTPEllipse(tKuka,c,dir,ratio,TefTool)

%% Arreguments:
% tKuka: TCP/IP connection obejct.
% c: 3x1 vector, is a vector defining the displacment of the the center
% point of the ellipse in relation to the start point of the ellipse,
% position units in (mm).
% dir: the direction vector of the (a) axis of the ellipse.
% ratio: the radious ratio (a/b) of the ellipse.
% theta: is the angle of the ellipe part , for drawing a complete
% ellipse this angle is equal to 2*pi.
% TefTool: 4X4 Transform matrix of the end-efector in the reference frame
% of the flange of the KUKA iiwa, position units in (mm).
% accel: is the acceleration (mm/sec2).
% velocity: is the motion velocity (mm/sec).

% Copy right: Mohammad SAFEEA
% 12th-Nov-2017

%% Convert inputs to a column vectors
c=colVec(c);
dir=colVec(dir);
ratio=colVec(ratio);
theta=colVec(theta);
velocity=colVec(velocity);
accel=colVec(accel);

%% Check input variables
if(size(c,1)~=3)
fprintf('Error: the vector (c) shall be a 3x1 vector \n');
return;
end

if(size(dir,1)~=3)
fprintf('Error: the vector (dir) shall be a 3x1 vector \n');
return;
end

if(size(ratio,1)~=1)
fprintf('Error: the ratio shall be a scalar  \n');
return;
end

if(size(theta,1)~=1)
fprintf('Error: the variable (theta) shall be a scalar  \n');
return;
end

if(size(TefTool,1)~=4)
fprintf('Error: the transform matrix (TefTool) shall be 4x4  \n');
return;
end

if(size(TefTool,2)~=4)
fprintf('Error: the transform matrix (TefTool) shall be 4x4  \n');
return;
end

if(size(velocity,1)~=1)
fprintf('Error: the velocity on the path shall be a scalar  \n');
return;
end

if(size(accel,1)~=1)
fprintf('Error: the acceleration on the path shall be a scalar  \n');
return;
end

%% Convert the units to (meter)
unitConverter=1000;
c=c/unitConverter;
velocity=velocity/unitConverter;
accel=accel/unitConverter;
TefTool(1:3,4)=TefTool(1:3,4)/unitConverter;

%% start the direct servo
     realTime_startDirectServoJoints(tKuka);

%% get current joints angles of robot
    jPos  = getJointsPos( tKuka );

%% Calculate the current position of TCP point of the robot (p0)
    qs=zeros(7,1);
for i=1:7
    qs(i)=jPos{i};
end
        T0=directKinematics(qs,TefTool); % TCP point transformation matrix
        p0=T0(1:3,4); % Current position of TCP point of the robot
        Tt=T0;

%% Calculate the parameters of the ellipse,
p=p0; % this is the starting point of the ellipse
c=c+p; % calculate the center coordinates from the center displacment vector
[R,theta0,a,b,c,errorFlag]=getEllipseParameters(p,c,ratio,dir);
if errorFlag==true % this happens when the start point is on the (dir) vector direction.
    % in such case the ellipse plane can not be specified in space.
    realTime_stopDirectServoJoints(tKuka);
    beep;
    fprintf('Error executing the ellipse function:\n');
    fprintf('Starting point of the ellipse shall not be on the line aligned with the (dir) vector and passing from the center (c)\n');
    return;
end

%% Get the length of the ellipse arc as a function of the parametric angle (theta)
theta1=theta0+theta;
[ thetaVec,sVec ] = getEllipseLengthVector( a,b, theta0,theta1 );
sizesVec=max(size(sVec));
L=sVec(end); % L is the length of the elliptic curve

%% Calculate the times:
[t0,t1,t2]=calculateInterpolationTimes(L,velocity,accel);

%% Joint space control

    [Ttemp,J]=directKinematics(qs,TefTool);
    vec=Ttemp(1:3,4);

%% dls solver parameters
 numberOfIterationForSolver=100;
 lambda=0.5;
 TefTool=eye(4);

    sCoordinate=0;
    interpolationCounter=1;
    thetaVar=theta0;

%% start of time
        dateVector0=datevec(now);
        time0=dateVector0(6)+dateVector0(5)*60+dateVector0(4)*60*60; % calculate time at this instant
while true

    % Calculate the elapsed time
        dateVector=datevec(now);
        timeNow=dateVector(6)+dateVector(5)*60+dateVector(4)*60*60; % calculate time at this instant
        deltaT=timeNow-time0; % elapsed is zero at first excution

    % calculate position of servo point
    if deltaT<t0
        sCoordinate=0.5*accel*deltaT*deltaT;
    elseif(deltaT<t1)
        s0=0.5*accel*t0*t0;
        sCoordinate=velocity*(deltaT-t0)+s0;
    elseif(deltaT<t2)
           s0=0.5*accel*t0*t0;
           s1=velocity*(t1-t0)+s0;
           sCoordinate=velocity*(deltaT-t1)-0.5*accel*(deltaT-t1)*(deltaT-t1)+s1;
    end
 %% when time ends break the loop
    if deltaT>t2
        break;
    end
%% Interpolate theta from the sVec,thetaVec
    for counter=interpolationCounter:(sizesVec-1)
        if(sCoordinate>sVec(counter))
            tetaRange=thetaVec(counter+1)-thetaVec(counter);
            sRange=sVec(counter+1)-sVec(counter);
            thetaVar=(sCoordinate-sVec(counter))...
                *tetaRange/sRange+thetaVec(counter);
            interpolationCounter=counter;
        end
    end

    pPrime=R*ellipseParametricFunction( a,b,thetaVar );
    p=c+pPrime;

    % calculate target transform
    Tt(1:3,4)=p;

    [ qs ] = kukaDLSSolver( qs, Tt, TefTool,numberOfIterationForSolver,lambda );

    for i=1:7
        jPos{i}=qs(i);
    end

	%% Send joint positions to robot
	sendJointsPositions( tKuka,jPos);

end

    % Stop the reltime motion
    realTime_stopDirectServoJoints(tKuka);

end

function [R,theta0,a,b,c,errorFlag]=getEllipseParameters(p,c,ratio,dir)

%% About
% This function returns the parameters of the ellipse

%% Arreguments:
% p: the starting point of the ellipse
% c: the center of the ellipse
% ratio: the ratio (a/b) of the ellipse
% dir: direction of the

%% Return values:
% R, the rotation matrix from the ellipse frame to the base frame of the
% robot
% theta0: the starting angle of the ellipse (given by the parametric equation of the ellipse)
% a,b: the big and the small radious of the ellipse
% c: is the center of the ellipse
% error flag: returns true if an error happens, this error is when the
% vector (p-c) is coinsident with the vector (dir).

p=colVec(p);
c=colVec(c);
dir=colVec(dir);

% X dirction vector
i=dir/norm(dir);

u=p-c;
%
k=cross(u,i);
normk=norm(k);

if(normk==0)
    errorFlag=true;
else
    errorFlag=false;
end

k=k/norm(k);
j=cross(k,i);

% rotation matrix from ellipse frame to base frame
R=[i,j,k];
% calculate a,b minimum,maximum radious of the ellipse
x=u'*i; % x coordinate of the starting point of the ellipse in the ellipse frame
y=u'*j; % y coordinate of the starting point of the ellipse in the ellipse frame
a=(x*x+y*y/(ratio*ratio))^0.5;
b=ratio*a;
% calculate theta0, begining angle of the ellipse
theta0=atan2(y,x);

end

function y=colVec(x)
    if size(x,2)==1
        y=x;
    else
        y=x';
    end
end

function [ thetaVec,sVec ] = getEllipseLengthVector( a,b, theta0,theta1 )
%% About:
% This function is used to return the length of the ellipse as a function
% of the angle theta.

%% Inputs:
% theta0, theta1: start/end angle of the arc of the ellipse (from the parametric equation of the ellipse)
% a,b: the big and the small radious of the ellipse.

thetaSpan = [theta0 theta1];
s0 = 0;
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[thetaVec,sVec] = ode45(@(theta,s) ((a*a*sin(theta)*sin(theta)+b*b*cos(theta)*cos(theta))^0.5), thetaSpan, s0,opts);

end

function [ vec ] = ellipseParametricFunction( a,b,theta )
%% About:
% Calculate the X and Y coordinates of the ellipse point corrsponding to
% angle (theta) from the parametric equation

%% inputs:
% theta: angle of point on the ellipse (from the parametric equation of the ellipse)
% a,b: the big and the small radious of the ellipse

%% Return value:
% The return value is a column vector, vec=[x;y;0] where:
% x: the x coordiante of the point on the ellipse correspodning to angle
% (theta) as described in the frame of the ellipse
% y: the y coordiante of the point on the ellipse  correspodning to angle
% (theta) as described in the frame of the ellipse

x=a*cos(theta);
y=b*sin(theta);
vec=[x;y;0];

end

function [t0,t1,t2]=calculateInterpolationTimes(L,v,a)
%% About:
% calculates the times of the acceleration, stable motion, deceleration

%% Arreguments:
% L: length of the curve (mm).
% v: linear velocity of eef on the the curve (mm/sec).
% a: linear acceleration of eef on the the curve (mm/sec2).

%% Return values:
% t0; time for acceleration.
% t1: time for acceleration, constant velocity.
% t2: time for acceleration, constant velocity and for deceleration

taw=v/a;

if(L>a*taw*taw)
    t0=taw;
    deltat=(L-a*taw*taw)/v;
    t1=t0+deltat;
    t2=t1+taw;
else
    taw1=(L/a)^0.5;
    t0=taw1;
    t1=taw1;
    t2=2*taw1;
end

end

repo: KST-Kuka-Sunrise-Toolbox


text:
function [afresults, AfAnalysisWindows, AFfile] = PerformAFdetection(subjectID,trr,rr,sqi,HRVparams)
%   PerformAFdetection(subjectID,trr,rr,HRVparams)
%
%	OVERVIEW:
%       Perform Atrial Fibrillation (AF) detection
%
%   INPUT:
%       subjectID : string containing the identifier of the subject to be analyze
%       trr      : (seconds)a single row of time of the rr interval
%                   data
%       rr        : (seconds) a single row of peak to peak interval
%                   data
%       sqi       : Signal Quality Index; Requires a matrix with
%                   at least two columns. Column 1 should be timestamps
%                   should be timestamps of each sqi measure, and Column 2
%                   should be SQI on a scale from 0 to 1.
%       HRVparam  : struct of settings for hrv_toolbox analysis
%
%   OUTPUT:
%       afresults : a single row containing a flag (1) when AF is
%                   dettected in a window and 0 if no AF
%
%   DEPENDENCIES & LIBRARIES:
%       PhysioNet Cardiovascular Signal Toolbox
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%
%   REFERENCE:
%   Vest et al. "An Open Source Benchmarked HRV Toolbox for Cardiovascular
%   Waveform and Interval Analysis" Physiological Measurement (In Press), 2018.
%
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Giulia Da Poian (giulia.dap@gmail.com) on Sep 6, 2017.
%       Dependent scripts written by various authors
%       (see functions for details)
%	COPYRIGHT (C) 2018
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Modified on 02.14.2018 to include SQI check before evaluating AF; if a
%   windows contains low SQI signal do not compute AF and mark the window
%   as NaN

if isempty(sqi)
     sqi(:,1) = trr;
     sqi(:,2) = ones(length(trr),1);
end

% 1. Calculate AF Features
AfAnalysisWindows = CreateWindowRRintervals(trr,rr,HRVparams,'af');
NNsamps = rr .* HRVparams.Fs;

AFtest = nan(length(AfAnalysisWindows),1);

for idx = 1:length(AfAnalysisWindows)
    tstart = AfAnalysisWindows(idx);

    if ~isnan(tstart)

        idxInWin = find(trr >= tstart & trr< tstart + HRVparams.af.windowlength);

        % Added by Giulia: exclude from the Analysis low quality
        % segments (SQI< SQI_threshold)
        sqiWin = sqi(sqi(:,1) >= tstart & sqi(:,1) < tstart + HRVparams.af.windowlength,2);
        LowQualityIdxs = find(sqiWin < HRVparams.sqi.LowQualityThreshold);

        % If enough data has an adequate SQI, perform the calculations
        cond1 = (numel(LowQualityIdxs)/length(sqiWin)) <  HRVparams.RejectionThreshold;
        % RR interval time series must be > 12 ans <60 to extract feautures
        cond2 = (length(NNsamps(idxInWin)) > 12 && length(NNsamps(idxInWin)) < 60);

        if (cond1 && cond2)
            features_af = AF_features(NNsamps(idxInWin),HRVparams.Fs);
            AFtest(idx) = SVM_AFdetection_withoutTrainingModel(features_af,1);
        end
    end
end

% 2. Export AF Data as CSV File

afresults = AFtest(:);
afcol_titles = {'AFtest'};

outputType = 'AF';
AFfile = SaveHRVoutput(subjectID,AfAnalysisWindows,afresults, ...
    afcol_titles, outputType, HRVparams, trr, rr);

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function varargout = pow10(varargin)
%POW10 (overloaded)

varargout{1} = 10.^varargin{1};

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [ y, m, d, f ] = now_to_ymdf_common ( )

%*****************************************************************************80
%
%% NOW_TO_YMDF_COMMON expresses the current date as a Common YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer Y, M, D, real F,
%    the YMDF date.
%
  c = clock ( );

  y = c(1);
  m = c(2);
  d = c(3);
  h = c(4);
  n = c(5);
  s = c(6);
  mu = 0.0;

  f = mu;
  f = s + f / 1000.0;
  f = n + f / 60.0;
  f = h + f / 60.0;
  f = f / 24.0;

  return
end


text:
classdef NoOptimizationVariablesValidator < AbstractLaunchVehicleDataValidator
    %NoOptimizationVariablesValidator Summary of this class goes here
    %   Detailed explanation goes here

    properties
        lvdData LvdData
    end

    methods
        function obj = NoOptimizationVariablesValidator(lvdData)
            obj.lvdData = lvdData;
        end

        function [errors, warnings] = validate(obj)
            errors = LaunchVehicleDataValidationError.empty(0,1);
            warnings = LaunchVehicleDataValidationWarning.empty(0,1);

            varSet = obj.lvdData.optimizer.vars;
            numVars = length(varSet.getTotalScaledXVector());
            if(numVars == 0)
                str = sprintf('No optimization variables enabled on script.');
                warnings(end+1) = LaunchVehicleDataValidationWarning(str);
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function show_per_class_table( soa, soa_ids, res_id, measure )

% Pascal classes names
classes={'aeroplane', 'bicycle'  , 'bird'       , 'boat'       ,...
         'bottle'   , 'bus'      , 'car'        , 'cat'        ,...
         'chair'    , 'cow'      , 'diningtable', 'dog'        ,...
         'horse'    , 'motorbike', 'person'     , 'pottedplant',...
         'sheep'    , 'sofa'     , 'train'      , 'tvmonitor'  };

% Allocate
all_res = zeros(length(soa_ids),length(classes)+1);
ncands = zeros(length(res_id),1);

% Fill
for ii=1:length(soa_ids)
    ncands(ii) = soa.(soa_ids{ii}).mean_n_masks(res_id(ii));
    for kk=1:length(classes)
        if strcmp(measure,'jaccard_object')
            res = soa.(soa_ids{ii}).per_class_results{kk}.meanmax;
        else
            res = soa.(soa_ids{ii}).per_class_results{kk}.global_J;
        end
        all_res(ii,kk) = res(res_id(ii));
    end
    res = soa.(soa_ids{ii}).(measure);
    all_res(ii,end) = res(res_id(ii));
end

% Get the maximums to set in boldface
[~,max_ids] = max(all_res,[],1);
disp(['MCG is the best in ' num2str(sum(max_ids(1:end-1)==1)) ' categories'])

% Show header
to_disp = '    Method &   NCands & ';
for kk=1:length(classes)
    to_disp = [to_disp sprintf('%14s',classes{kk}) ' & ']; %#ok<AGROW>
end
to_disp = [to_disp '        Global \\'];
disp(to_disp)

% Show table
for ii=1:length(soa_ids)
    curr_n_reg = soa.(soa_ids{ii}).mean_n_masks(res_id(ii));
    to_disp = [sprintf('%10s', soa_ids{ii}) ' & ' sprintf('%4s%d',' ', round(curr_n_reg)) ' & '];
    for kk=1:length(classes)
        if (ii==max_ids(kk))
            to_disp = [to_disp sprintf('%13s',['\textbf{' sprintf('%2.1f',100*all_res(ii,kk))]) '} & ']; %#ok<AGROW>
        else
            to_disp = [to_disp sprintf('%13s',sprintf('%2.1f',100*all_res(ii,kk))) '  & ']; %#ok<AGROW>
        end
    end
    if (ii==max_ids(end))
        to_disp = [to_disp sprintf('%13s',['\textbf{' sprintf('%2.1f',100*all_res(ii,end))]) '} \\ ']; %#ok<AGROW>
    else
        to_disp = [to_disp sprintf('%13s',sprintf('%2.1f',100*all_res(ii,end))) '  \\']; %#ok<AGROW>
    end
    disp(to_disp)
end

end

author: jponttuset

repo: mcg

save_path: github-repos/MATLAB/jponttuset-mcg


text:
function o = symmetrise(o,varargin)
% all crystallographically equivalent orientations

if nargin > 1 && isa(varargin{1},'symmetry')
  CS = varargin{1};
  SS = getClass(varargin(2:end),'symmetry',specimenSymmetry);
else
  CS = o.CS;
  SS = o.SS;
end

if check_option(varargin,'proper')
  CS = CS.properGroup;
  SS = SS.properGroup;
end
o = symmetrise@rotation(o,CS,SS);

if o.antipodal, o = [o;inv(o)]; end

if check_option(varargin,'unique')
  [~,ind] = unique(o,'noSymmetry');
  o = subSet(o,ind);
end

author: mtex-toolbox

repo: mtex

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


text:
function [ o, x, w ] = cn_geg_02_xiu ( n, alpha )

%*****************************************************************************80
%
%% CN_GEG_02_XIU implements the Xiu rule for region CN_GEG.
%
%  Discussion:
%
%    The rule has order
%
%      O = N + 1.
%
%    The rule has precision P = 2.
%
%    CN_GEG is the cube [-1,+1]^N with the Gegenbauer weight function
%
%      w(alpha;x) = product ( 1 <= i <= n ) (1-x(i)^2)^alpha.
%
%    with -1.0 < alpha.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dongbin Xiu,
%    Numerical integration formulas of degree two,
%    Applied Numerical Mathematics,
%    Volume 58, 2008, pages 1515-1520.
%
%  Parameters:
%
%    Input, integer N, the spatial dimension.
%
%    Input, real ALPHA, the parameter.
%    -1.0 < ALPHA.
%
%    Input, integer O, the order.
%
%    Output, real X(N,O), the abscissas.
%
%    Output, real W(O), the weights.
%
  if ( alpha <= -1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'CN_GEG_02_XIU - Fatal error!\n' );
    fprintf ( 1, '  ALPHA <= -1.0\n' );
    error ( 'CN_GEG_02_XIU - Fatal error!' );
  end

  o = n + 1;
  x = zeros ( n, o );
  w = zeros ( o, 1 );

  for j = 1 : o

    i = 0;
    for r = 1 : floor ( n / 2 )
      arg = 2 * r * ( j - 1 ) * pi / ( n + 1 );
      i = i + 1;
      x(i,j) = sqrt ( 2.0 ) * cos ( arg );
      i = i + 1;
      x(i,j) = sqrt ( 2.0 ) * sin ( arg );
    end

    if ( i < n )
      i = i + 1;
      x(i,j) = r8_mop ( j - 1 );
    end

  end

  gamma0 = 1.0;
  delta0 = 0.0;
  c1 = 1.0 / ( 2.0 * alpha + 3.0 );

  x(1:n,1:o) = ( sqrt ( gamma0 * c1 ) * x(1:n,1:o) - delta0 ) / gamma0;

  expon = 0;
  volume_1d = c1_geg_monomial_integral ( alpha, expon );
  volume = volume_1d ^ n;

  w(1:o) = volume / o;

  return
end

author: johannesgerer

repo: jburkardt-m


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

%%
clear; close all; clc;

%% Syntax
% |[F,V]=import_obj_geom(fileName);|

%% Description
% This function imports the geometry (faces and nodes0 contained in an OBJ
% file. All texture/material data is ignored.

%% Examples
%

defaultFolder = fileparts(fileparts(mfilename('fullpath')));
loadPath=fullfile(defaultFolder,'data','OBJ');
fileName=fullfile(loadPath,'test.obj');

[F,V]=import_obj_geom(fileName);

%%
% Visualize

cFigure;
gpatch(F,V,'w','k');
axisGeom; camlight headlight;
gdrawnow;

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