text:
function [rpav]=realized_preaveraged_variance(price,time,timeType,samplingType,samplingInterval,options)
% Estimated quadratic variation using Preaveraged Realized Variance
%
% USAGE:
%   [RPAV] = realized_preaveraged_variance(PRICE)
%   [RPAV] = realized_preaveraged_variance(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINTERVAL)
%   [RPAV] = realized_preaveraged_variance(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINTERVAL,OPTIONS)
%
% INPUTS:
%   PRICE            - m by 1 vector of high frequency prices
%   TIME             - [OPTIONAL] m by 1 vector of times where TIME(i) corresponds to PRICE(i).
%   TIMETYPE         - [OPTIONAL] String describing the way times are measured
%                       'wall'    24-hour clock of the form HHMMSS.mmm, e.g. 101543 or 153217
%                       'seconds' Time measured in seconds past midnight
%                       'unit'  Unit normalized date format, e.g. .1, .234, .9
%                         Unit normalized times are more general than the other types and can be
%                         applied to data from more than one calendar day
%   SAMPLINGTYPE     - [OPTIONAL] String describing the type of sampling to use when
%                        filtering PRICE
%                        'CalendarTime' - Sample in calendar time using observations separated by
%                          SAMPLINGINTERVAL seconds
%                        'CalendarUniform' - Sample in calendar time using SAMPLINGINTERVAL
%                          observations spread uniformly between TIME(1) and TIME(m)
%                        'BusinessTime' - Sample in business (tick) time using observation separated
%                          by SAMPLINGINTERVAL ticks
%                        'BusinessUniform' - Sample in business (tick) time using observations
%                          uniformly spaced in business time.
%                        'Fixed' - Sample at specific points in time. When using fixed,
%                          SAMPLINGINTERVAL must be a n by 1 vector of times with the same TIMETYPE
%                          as TIME (i.e. seconds if TIME is in seconds)
%   SAMPLINGINTERVAL  - [OPTIONAL] Scalar integer or n by 1 vector whose meaning depends on the
%                         selected SAMPLINGTYPE
%   OPTIONS           - [OPTIONAL] Preaveraged Realized Variance option structure initialized by calling
%                         realized_options('Preaveraging'). See help realized_options for a description of
%                         available options.
%
% OUTPUTS:
%   RK          - Preaveraged realized variance estimate
%
% COMMENTS:
%  Follows Christensen, Oomen and Podolski (2014) mostl closely, with the
%  noise variance estimator usef in Hautsch and Podolski (2013)
%
% EXAMPLES:
%
%  See also REALIZED_OPTIONS, REALIZED_KERNEL, REALIZED_NOISE_ESTIMATE, REALIZED_VARIANCE,
%  REALIZED_VARIANCE_OPTIMAL_SAMPLING, REALIZED_RANGE, REALIZED_QUANTILE_VARIANCE

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 2/27/2014

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 1
        m = length(price);
        time = linspace(9.5*3600,16*3600,m)';
        timeType = 'seconds';
        samplingType = 'businesstime';
        samplingInterval = 1;
        options = realized_options('preaveraging');
    case 5
        options = realized_options('preaveraging');
    case 6
        % Nothing
    otherwise
        error('One, five or six inputs required.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Need method to estimate omega, and options.

% Inserted to protect against inputing integer times
time = double(time);

% 1. Filter the rpice
filteredPrice = realized_price_filter(price,time,timeType,samplingType,samplingInterval);
returns = diff(log(filteredPrice));
% 2. Compute rpreaveraged returns
theta = options.theta;
m = length(returns);
K = ceil(theta * sqrt(m));
g = @(x) min(x,1-x);
w = g((1:(K-1))/K);

preav_returns = nan(m-K+2,1);
for i=1:m-K+2
    preav_returns(i) = w*returns(i:i+K-2);
end

% 3. Compute constants
psi_1 = K*sum((g((1:K)/K) - g((0:(K-1))/K)).^2);
psi_2 = 1/K*sum((g((1:(K-1))/K)).^2);
% 4. Estimate noise
[noiseVariance, ~, ~, noiseEstimateOomen] = realized_noise_estimate(price, time, timeType, options);

omega = noiseEstimateOomen;
if omega<0
    omega = noiseVariance;
end
% 5. Compute Preaveraed Variance
const1 = m/(m-K+2);
const2 = 1/(K*psi_2);
bias = psi_1/(theta.^2 * psi_2) * omega^2;

rpav = const1 * const2 * sum(preav_returns.^2) - bias;

author: bashtage

repo: mfe-toolbox


text:
function x = vec(X)
x = X(:);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function success=testClassifierImport()

testDataDirName='/groups/branson/bransonlab/projects/JAABA/test_data';
jabFileName=fullfile(testDataDirName,'larping.jab');
nameOfJabFileWithClassifierToImport=fullfile(testDataDirName,'larping_different_features.jab');
nameOfExpDir='GMR_71G01_AE_01_TrpA_Rig1Plate15BowlA_20120316T144027';
gtMode=false;
data=JLabelData('setstatusfn',@(str)(fprintf('%s\n',str)), ...
                'clearstatusfn',@()(nop()), ...
                'isInteractive',true);
data.openJabFile(jabFileName,gtMode);

iExp=whichstr(nameOfExpDir,data.expnames);
if isempty(iExp)
  error('The experiment seems to be missing.');
elseif length(iExp)>1
  error('The experiment seems to be present more than once in the JLabelData object.');
end
iFly=1;
data.setCurrentTarget(iExp,iFly);

% The prediction should be larping for these frames
ts=3900:3940;
larpingBehavior=whichstri('larping',data.labelnames);
data.Predict(iExp,iFly,ts(1),ts(end));  % Make sure the prediction has been done
pred=data.GetPredictedIdx(iExp,iFly,ts(1),ts(end));
if all(pred.predictedidx==larpingBehavior)
  fprintf('Predictions are correct prior to importing classifier!\n');
else
  error('Predictions are wrong prior to importing classifier');
end

% Do the import
data.importClassifier(nameOfJabFileWithClassifierToImport);

% The prediction should be none for these frames
noneBehavior=whichstri('none',data.labelnames);
data.Predict(iExp,iFly,ts(1),ts(end));  % Make sure the prediction has been done
pred=data.GetPredictedIdx(iExp,iFly,ts(1),ts(end));
if all(pred.predictedidx==noneBehavior)
  fprintf('Predictions are correct after importing classifier!\n');
else
  error('Predictions are wrong after importing classifier');
end

data.closeJabFile();
data=[];  %#ok
success=true;

end

author: kristinbranson

repo: JAABA


text:
classdef opticheckval
    %OPTICHECKVAL Contains Static Methods for Checking Option Values

    methods(Static)

        function err = checkScalarSet(value,field,set)
            %Scalar in set
            err = [];
            if(~isscalar(value) || ~isnumeric(value) || ~isa(value,'double') || ~any(value==set))
                str = sprintf('\n');
                for i = 1:length(set)
                    str = sprintf('%s %g\n',str,set(i));
                end
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should a double scalar from one of the following supported values:\n%s',field,str);
            end
        end

        function err = checkScalar01(value,field)
            %Scalar 0 or 1
            err = [];
            if(~isscalar(value) || ~isnumeric(value) || (value ~= 0 && value ~=1) || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be an integer scalar with the value 0 or 1',field);
            end
        end

        function err = checkVector01(value,field)
            %Vector 0 or 1
            err =[];
            if(~isnumeric(value) || any((value == 0 | value == 1) ~= 1) || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a vector with the values 0 or 1',field);
            end
        end

        function err = checkScalarGrtZ(value,field)
            %Scalar double > Zero
            err = [];
            if(~isscalar(value) || ~isnumeric(value) || value <= 0 || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a double scalar greater than 0',field);
            end
        end

        function err = checkScalarIntGrtZ(value,field)
            %Scalar integer > zero
            err = [];
            if(~isscalar(value) || ~isnumeric(value) || value < 1 || (round(value) ~= value) || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be an integer scalar greater than 0',field);
            end
        end

        function err = checkScalarNonNeg(value,field)
            %Scalar double >= Zero
            err = [];
            if(~isscalar(value) || ~isnumeric(value) || value < 0 || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a scalar greater or equal to 0',field);
            end
        end

        function err = checkScalarIntNonNeg(value,field)
            %Scalar integer >= zero
            err = [];
            if(~isscalar(value) || ~isnumeric(value) || value < 0 || (round(value) ~= value) || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be an integer scalar greater than 0',field);
            end
        end

        function err = checkScalarDbl(value,field)
            %Scalar double
            err = [];
            if(~isscalar(value) || ~isnumeric(value) || ~isreal(value) || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a scalar double',field);
            end
        end

        function err = checkScalarBoundLELE(value,field,lb,ub)
            %Scalar double lb <= value <= ub
            err = [];
            if(~isscalar(value) || ~isnumeric(value) || ~isreal(value) || value < lb || value > ub || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a scalar double %g <= value <= %g',field,lb,ub);
            end
        end

        function err = checkScalarBoundLLE(value,field,lb,ub)
            %Scalar double lb < value <= ub
            err = [];
            if(~isscalar(value) || ~isnumeric(value) || ~isreal(value) || value <= lb || value > ub || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a scalar double %g < value <= %g',field,lb,ub);
            end
        end

        function err = checkScalarBoundLEL(value,field,lb,ub)
            %Scalar double lb <= value < ub
            err = [];
            if(~isscalar(value) || ~isnumeric(value) || ~isreal(value) || value < lb || value >= ub || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a scalar double %g <= value < %g',field,lb,ub);
            end
        end

        function err = checkScalarBoundLL(value,field,lb,ub)
            %Scalar double lb < value < ub
            err = [];
            if(~isscalar(value) || ~isnumeric(value) || ~isreal(value) || value <= lb || value >= ub || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a scalar double %g < value < %g',field,lb,ub);
            end
        end

        function err = checkScalarIntBoundLEL(value,field,lb,ub)
            %Scalar integer lb <= value < ub
            err = [];
            if(~isscalar(value) || ~isnumeric(value) || ~isreal(value) || (round(value) ~= value) || value < lb || value >= ub || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a scalar integer %g <= value < %g',field,lb,ub);
            end
        end

        function err = checkScalarIntBoundLELE(value,field,lb,ub)
            %Scalar integer lb <= value <= ub
            err = [];
            if(~isscalar(value) || ~isnumeric(value) || ~isreal(value) || (round(value) ~= value) || value < lb || value > ub || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a scalar integer %g <= value <= %g',field,lb,ub);
            end
        end

        function err = checkVectorIntGrtZ(value,field)
            %Vector of integers > 0
            err = [];
            if(~isnumeric(value) || any(value < 1) || any(round(value) ~= value) || ~any(isa(value,'double')))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be an integer vector with all elements greater than 0',field);
            end
        end

        function err = checkDualCol(value,field)
            %Double matrix with two columns
            err = [];
            if(~isnumeric(value) || ~isreal(value) || size(value,2) ~= 2 || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a double matrix with two columns',field);
            end
        end

        function err = checkDblVec(value,field)
            %Double vector (not sparse)
            err = [];
            if(~isnumeric(value) || ~isreal(value) || ~isa(value,'double') || (size(value,1) > 1 && size(value,2) > 1) || issparse(value))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a real double dense vector (1D)',field);
            end
        end

        function err = checkDblMat(value,field)
            %Double Matrix or Vector or Scalar
            err = [];
            if(~isnumeric(value) || ~isreal(value) || ~isa(value,'double'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a real double matrix or vector (2D or 1D)',field);
            end
        end

        function err = checkDblVecOrCell(value,field)
            %Double Vector or Cell Array of Double Vectors
            err = [];
            if(~iscell(value)), value = {value}; end
            for i = 1:length(value)
                val = value{i};
                if(~isnumeric(val) || ~isreal(val) || ~isa(val,'double') || (size(val,1) > 1 && size(val,2) > 1) || issparse(val))
                    err = MException('OPTI:SetFieldError','Parameter ''%s'' (cell %d) should be a real double dense vector or scalar (1D)',field,i);
                end
            end
        end

        function err = checkDblVecOrLogVec(value,field)
            %Double Vector or Logical Vectors
            err = [];
            if(islogical(value))
                if(size(value,1) > 1 && size(value,2) > 1)
                    err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a logical vector OR real double dense vector (1D)',field);
                end
            else
                if(~isnumeric(value) || ~isreal(value) || ~isa(value,'double') || (size(value,1) > 1 && size(value,2) > 1) || issparse(value))
                    err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a real double dense vector (1D) or logical vector',field);
                end
            end
        end

        function err = checkDblMatOrCell(value,field)
            %Double Matrix or Cell Array of Double Matrices
            err = [];
            if(~iscell(value)), value = {value}; end
            for i = 1:length(value)
                val = value{i};
                if(~isnumeric(val) || ~isreal(val) || ~isa(val,'double'))
                    err = MException('OPTI:SetFieldError','Parameter ''%s'' (cell %d) should be a real double matrix or vector (2D or 1D)',field,i);
                end
            end
        end

        function err = checkChar(value,field)
            %Char array
            err = [];
            if(~ischar(value))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a char array (string)',field);
            end
        end

        function err = checkNLOPTAlg(value,field)
            %Char array and valid NLOPT algorithm
            err = [];
            if(~ischar(value))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a char array (string)',field);
            else
                nloptSolver(value); %will throw an error if not found
            end
        end

        function err = checkOptiSolver(value,field)
            %Char array and valid OPTI solver
            err = [];
            if(~ischar(value))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a char array (string)',field);
            else
                rub = checkSolver(value,1); %#ok<NASGU>
            end
        end

        function err = checkYesNo(value,field)
            %Char array and yes or no
            err = opticheckval.checkValidString(value,field,{'yes','no'});
        end

        function err = checkYesCrashNo(value,field)
            %Char array and no, or yes will crash
            err = [];
            if(~ischar(value) || ~strcmpi(value,'no'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' must be set as no, otherwise OPTI may crash (sorry!).\n',field);
            end
        end

        function err = checkEnDis(value,field)
            %Char array and enable or disable
            err = opticheckval.checkValidString(value,field,{'enable','disable'});
        end

        function err = checkOnOff(value,field)
            %Char array and on or off
            err = opticheckval.checkValidString(value,field,{'on','off'});
        end

        function err = checkValidString(value,field,valid)
            %Char array and valid string
            err = [];
            if(~ischar(value) || all(strcmpi(value,valid)==0))
                str = sprintf('\n');
                for i = 1:length(valid)
                    str = sprintf('%s ''%s''\n',str,valid{i});
                end
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should a string from the following supported options:\n%s',field,str);
            end
        end

        function err = checkValidStringArr(value,field,valid)
            %Char array and of valid characters
            err = [];
            if(~ischar(value))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a char array (string)',field);
            else
                if(~iscell(valid)), valid = {valid}; end
                for i = 1:length(value)
                   if(~any(ismember(valid,value(i))))
                       str = sprintf('\n');
                        for j = 1:length(valid)
                            str = sprintf('%s ''%s''\n',str,valid{j});
                        end
                        err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a char array (string) with the following characters only:\n%s',field,str);
                        return;
                    end
                end
            end
        end

        function err = checkStruct(value,field)
            %Structure
            err = [];
            if(~isstruct(value))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a structure\n',field);
            end
        end

        function err = checkStructFields(value,field,reqfields)
            %Structure with required fields
            err = [];
            if(~isstruct(value))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a structure\n',field);
            else
                flds = fieldnames(value);
                if(~iscell(reqfields)), reqfields = {reqfields}; end
                for i = 1:length(reqfields)
                    if(~any(ismember(flds,reqfields{i})))
                        str = sprintf('\n');
                        for j = 1:length(reqfields)
                            str = sprintf('%s ''%s''\n',str,reqfields{j});
                        end
                        err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a structure with the following fields:\n%s',field,str);
                        return;
                    end
                end
            end
        end

        function err = checkFunHandle(value,field)
            %Function Handle
            err = [];
            if(~isa(value,'function_handle'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a function handle.\n\nFor existing m-file functions, pass as @function_name (not a string, using the ''@'' operator).\n',field);
            end
        end

        function err = checkNumorFHndl(value,field)
            %Function handle or numeric matrix
            err = [];
            if(isnumeric(value))
                if(~isreal(value) || ~isa(value,'double'))
                    err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a real, double, matrix/vector or function handle\n',field);
                end
            elseif(~isa(value,'function_handle'))
                err = MException('OPTI:SetFieldError','Parameter ''%s'' should be a function handle OR real, double, matrix\n',field);
            end
        end

    end

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%% Ke Ma, Christopher Bodden
% CS 766 - Project 2

%% simple callback for enlarging graphs/images
function enlargeImage(handles, falseColor)
im = get( gcbo,'cdata' );
imtool(im, [min(im(:)) max(im(:))] )

if exist('falseColor')
    colormap(jet(256));
end
end

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [ x, y, z, w ] = ld0770 ( )

%*****************************************************************************80
%
%% LD0770 computes the 770 point Lebedev angular grid.
%
%  Modified:
%
%    14 September 2010
%
%  Author:
%
%    Dmitri Laikov
%
%  Reference:
%
%    Vyacheslav Lebedev, Dmitri Laikov,
%    A quadrature formula for the sphere of the 131st
%    algebraic order of accuracy,
%    Russian Academy of Sciences Doklady Mathematics,
%    Volume 59, Number 3, 1999, pages 477-481.
%
%  Parameters:
%
%    Output, real X(N), Y(N), Z(N), W(N), the coordinates
%    and weights of the points.
%
  n = 0;
  x = zeros(770,1);
  y = zeros(770,1);
  z = zeros(770,1);
  w = zeros(770,1);
  a = 0.0;
  b = 0.0;
  v = 0.2192942088181184E-03;
  [ n, x, y, z, w ] = gen_oh ( 1, n, a, b, v, x, y, z, w );
  v = 0.1436433617319080E-02;
  [ n, x, y, z, w ] = gen_oh ( 2, n, a, b, v, x, y, z, w );
  v = 0.1421940344335877E-02;
  [ n, x, y, z, w ] = gen_oh ( 3, n, a, b, v, x, y, z, w );
  a = 0.5087204410502360E-01;
  v = 0.6798123511050502E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.1228198790178831;
  v = 0.9913184235294912E-03;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.2026890814408786;
  v = 0.1180207833238949E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.2847745156464294;
  v = 0.1296599602080921E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.3656719078978026;
  v = 0.1365871427428316E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.4428264886713469;
  v = 0.1402988604775325E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.5140619627249735;
  v = 0.1418645563595609E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6306401219166803;
  v = 0.1421376741851662E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6716883332022612;
  v = 0.1423996475490962E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.6979792685336881;
  v = 0.1431554042178567E-02;
  [ n, x, y, z, w ] = gen_oh ( 4, n, a, b, v, x, y, z, w );
  a = 0.1446865674195309;
  v = 0.9254401499865368E-03;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.3390263475411216;
  v = 0.1250239995053509E-02;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.5335804651263506;
  v = 0.1394365843329230E-02;
  [ n, x, y, z, w ] = gen_oh ( 5, n, a, b, v, x, y, z, w );
  a = 0.6944024393349413E-01;
  b = 0.2355187894242326;
  v = 0.1127089094671749E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2269004109529460;
  b = 0.4102182474045730;
  v = 0.1345753760910670E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.8025574607775339E-01;
  b = 0.6214302417481605;
  v = 0.1424957283316783E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.1467999527896572;
  b = 0.3245284345717394;
  v = 0.1261523341237750E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.1571507769824727;
  b = 0.5224482189696630;
  v = 0.1392547106052696E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.2365702993157246;
  b = 0.6017546634089558;
  v = 0.1418761677877656E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.7714815866765732E-01;
  b = 0.4346575516141163;
  v = 0.1338366684479554E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3062936666210730;
  b = 0.4908826589037616;
  v = 0.1393700862676131E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );
  a = 0.3822477379524787;
  b = 0.5648768149099500;
  v = 0.1415914757466932E-02;
  [ n, x, y, z, w ] = gen_oh ( 6, n, a, b, v, x, y, z, w );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function parent = connectedComponents(A)
% label connected componentes in an graph (adjacency matrix)

%elimination tree
parent = etree(A);

isleaf = parent ~= 0;
parent(~isleaf) = 1:nnz(~isleaf);
i = find(isleaf);
for i = i(end:-1:1)
  parent(i) = parent(parent(i));
end

author: mtex-toolbox

repo: mtex

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


text:
function res = gt(a,b)
%GT           Implements  a > b  for gradients, compares only a.x and b.x
%

% written  10/16/98     S.M. Rump
% modified 12/18/02     S.M. Rump  complex comparison
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if ~isa(a,'gradient')
    if isreal(a) & isreal(b.x)
      res = ( a>b.x );
    else
      res = ( real(a)>real(b.x) ) & ( imag(a)>imag(b.x) );
    end
  elseif ~isa(b,'gradient')
    if isreal(a.x) & isreal(b)
      res = ( a.x>b );
    else
      res = ( real(a.x)>real(b) ) & ( imag(a.x)>imag(b) );
    end
  else
    if isreal(a.x) & isreal(b.x)
      res = ( a.x>b.x );
    else
      res = ( real(a.x)>real(b.x) ) & ( imag(a.x)>imag(b.x) );
    end
  end

repo: OpenMAS


text:
function result = p00_even ( prob, int_num )

%*****************************************************************************80
%
%% P00_EVEN uses evenly spaced points to integrate a function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer PROB, the problem index.
%
%    Input, integer INT_NUM, the number of sample points.
%
%    Output, real RESULT, the approximate integral.
%
  [ a, b ] = p00_lim ( prob );

  if ( int_num == 1 )
    x(1) = ( a + b ) / 2.0;
  else
    x(1:int_num) = ( ( int_num - (1:int_num)     ) * a    ...
                   + (           (1:int_num) - 1 ) * b  ) ...
                   / ( int_num               - 1 );
  end

  result = sum ( p00_fun ( prob, int_num, x(1:int_num) ) );

  result = ( b - a ) * result / int_num;

  return
end


text:
function im = read_tiny_big_binary(words_in,ind)
%
%%% Routine to read in tiny images from giant binary file holding all images.
%
%%% There are two modes of operation:
%
% 1. raw indexing mode
% --------------------
% Inputs:
% 1. words_in - 1 x nImages vector, each element being an integer
% in the range 1 to 73,777,893 (or whatever the total number of images in the dataset is).
% Outputs:
% 1. im - 3072 x nImages uint8 data of images.
%
% ++++++++++++++++++++++++++++++++++++++++++++++++++++
%
% 2. word based indexing.
% -----------------------
% Inputs:
% 1. word_in - 1 x nWords cell array, each entry holding one noun
% 2. ind - 1 x nWords cell array, each holding indices of images
% for the corresponding noun
% Outputs:
% 1. im - 3072 x nImages uint8 data of images.
%
% Requires: read_tiny_binary_big_core.mex file
%
% Notes: (i) The code automatically sorts the indices so that they are in
%            order, to improve speed.
%
%       (ii) The access time per image drops as the total number of images to be
%            gathered increases. e.g. for 500,000 images load time is
%            1ms/image. For 100,000 images it is 2ms/image. So try and pass in big
%            blocks of indices.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SET PATHS BELOW TO tiny_images.bin AND tiny_index.mat

% path to giant binary file
binary_fname = '/mit/tiny/data/tiny_images.bin';
% indexing info .mat file
binary_data_fname = '/mit/tiny/data/tiny_index.mat';

% END
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% total # imgs
TOTAL_NUM_IMGS = 79302017;

if (nargin==1)
    RAW_INDEX_MODE = 1;
else
    RAW_INDEX_MODE = 0;
end

if (nargin==2)
    if iscell(ind)
        % do nothing
    else
        % convert to cell arrays
        words_in = {words_in};
        ind = {ind};
    end

    nWords = length(words_in);
end

if RAW_INDEX_MODE
    %% just go and call binary with raw indices
    %% check that all indices are >0 and <TOTAL_NUM_IMGS
    if ((sum(words_in<1)==0) & (sum(words_in>TOTAL_NUM_IMGS)==0))

        % now sort for speed
        [sorted_words,sort_ind] = sort(words_in);

        %% use MEX for sppedy reading in
        im = read_tiny_binary_big_core(binary_fname, uint64(sorted_words));

        % sort back
        im(:,sort_ind) = im;
    else
        error('indices <0 or greater than total number of images');
    end

else % words_in and indices passed in

    % load up index data
    load(binary_data_fname);

    % compute hashes for all words_in
    hashes = compute_hash_function(words_in,hash_key);

    img_offsets = [];

    for a=1:nWords

        % find word using hash
        match_ind = find(hashes(a) == hash);

        if any(match_ind)
            % check all ind are >=1 and <=total number of images

            if ((sum(ind{a}<1)==0) & (sum(ind{a}>num_imgs(match_ind))==0))
                word_offset = offset(match_ind);
                img_offsets = [img_offsets , uint64(double(word_offset) + [ind{a}])];
            else
                error('indices <0 or greater than number per image');
            end
        else
          %im = zeros(3072,1);
          %return;
          error('word not found');
        end

        % now sort for speed
        [sorted_offsets,sort_ind] = sort(img_offsets);

        im = read_tiny_binary_big_core(binary_fname, uint64(sorted_offsets));

        % sort back
        im(:,sort_ind) = im;
    end
end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function [sys,x0,str,ts]=NonlinearPD(t,x,u,flag,Beta1,Beta2,Delta,A1,A2,SampleTime)

if flag==0

    sys = [0;1;1;2;0;1;1];
     x0 = [0 ];
    str = [ ];
     ts = [SampleTime 0];

elseif flag==2

    Bili=fal(u(1),A1,Delta);
    Weifen=fal(u(2),A2,Delta);
    sys=Beta1*Bili+Beta2*Weifen;

elseif flag==3
    sys=x;

elseif flag==4
    sys=sys+SampleTime;

else
    sys=[];
end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function DispCoord1(Mode)

%--------------------------------------------------------------------------------------------
%DispCoord1(MODE) show the coordinate of the point on the 2D or 3D in 2 mode:
%   Mode=1, show in format (x,y) or (x,y)=z whenever you click a mouse button at one point;
%   Mode=2, show in format (x,y) or (x,y)=z in the mouse pointer trace when you move it.
%   Values are shown with 3 decimals.
%
%   Usage Examples,
%   x = -pi:.1:pi;
%   y = sin(x);
%   plot(x,y)
%   DispCoord1(2);
%
%   [X,Y] = meshgrid(-3:.125:3);
%   Z = peaks(X,Y);
%   meshc(X,Y,Z);
%   DispCoord1(1);
%
%   Zhenhai Wang <zhenhai@ieee.org>
%   Version 1.00
%   April, 2002
%--------------------------------------------------------------------------------------------

if Mode==1
    set(gcf,'Pointer','crosshair','WindowButtonDownFcn','DispCoord2');
elseif Mode==2
    set(gcf,'Pointer','crosshair','WindowButtonMotionFcn','DispCoord2');
end

author: Sable

repo: mcbench-benchmarks


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

%*****************************************************************************80
%
%% STUDENT_CDF_VALUES returns some values of the Student CDF.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      Needs["Statistics`ContinuousDistributions`"]
%      dist = StudentTDistribution [ df ]
%      CDF [ dist, x ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real C, is usually called the number of
%    degrees of freedom of the distribution.  C is typically an
%    integer, but that is not essential.  It is required that
%    C be strictly positive.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 13;

  c_vec = [ ...
    1.0, 2.0, 3.0, 4.0, ...
    5.0, 2.0, 5.0, 2.0, ...
    5.0, 2.0, 3.0, 4.0, ...
    5.0 ];

  fx_vec = [ ...
     0.6000231200328521, ...
     0.6001080279134390, ...
     0.6001150934648930, ...
     0.6000995134721354, ...
     0.5999341989834830, ...
     0.7498859393137811, ...
     0.7500879487671045, ...
     0.9500004222186464, ...
     0.9499969138365968, ...
     0.9900012348724744, ...
     0.9900017619355059, ...
     0.9900004567580596, ...
     0.9900007637471291 ];

  x_vec = [ ...
     0.325, ...
     0.289, ...
     0.277, ...
     0.271, ...
     0.267, ...
     0.816, ...
     0.727, ...
     2.920, ...
     2.015, ...
     6.965, ...
     4.541, ...
     3.747, ...
     3.365 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

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

  return
end


text:
function tests = test_isPointInEllipsoid
% Test suite for the file isPointInEllipsoid.
%
%   Test suite for the file isPointInEllipsoid
%
%   Example
%   test_isPointInEllipsoid
%
%   See also
%     isPointInEllipsoid

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2020-11-20,    using Matlab 9.8.0.1323502 (R2020a)
% Copyright 2020 INRAE - BIA-BIBS.

tests = functiontests(localfunctions);

function test_Inside(testCase) %#ok<*DEFNU>
% Test call of function without argument.

elli = [10 20 30   50 30 10   5 10 0];
p = [20 30 35];

b = isPointInEllipsoid(p, elli);

assertTrue(testCase, b);

function test_Outside(testCase) %#ok<*DEFNU>
% Test call of function without argument.

elli = [10 20 30   50 30 10   5 10 0];
p = [-20 10 25];

b = isPointInEllipsoid(p, elli);

assertFalse(testCase, b);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function G=gaussFilter(segma,kSize)
% Creates a 1-D Gaussian kernel of a standard deviation 'segma' and a size
% of 'kSize'.
%
% In theory, the Gaussian distribution is non-zero everywhere. In practice,
% it's effectively zero at places further away from about three standard
% deviations. Hence the reason why the kernel is suggested to be truncated
% at that point.
%
% The 2D Gaussian filter is a complete circular symmetric operator. It can be
% seperated into x and y components. The 2D convolution can be performed by
% first convolving with 1D Gaussian in the x direction and the same in the
% y direction.
%
% Author: Mohd Kharbat at Cranfield Defence and Security
% mkharbat(at)ieee(dot)org , http://mohd.kharbat.com
% Published under a Creative Commons Attribution-Non-Commercial-Share Alike
% 3.0 Unported Licence http://creativecommons.org/licenses/by-nc-sa/3.0/
%
% October 2008

if nargin<1
    segma=1;
end
if nargin<2
    kSize=2*(segma*3);
end

x=-(kSize/2):(1+1/kSize):(kSize/2);
G=(1/(sqrt(2*pi)*segma)) * exp (-(x.^2)/(2*segma^2));


text:
function J = computeCostMulti(X, y, theta)
% J = COMPUTECOSTMULTI(X, y, theta) computes the cost of using theta as the
% parameter for linear regression to fit the data points in X and y

m = length(y); % number of training examples

J = (1/(2*m)) * (X * theta - y)' * (X * theta - y); % Vectorized

end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
% M=1 scattering, cv parameters

run_name = 'DSS_Table2_TIMIT_m1';

src = phone_src('/path/to/timit');

[train_set,test_set,valid_set] = phone_partition(src);

N = 2^13;
T_s = 2560;

filt1_opt.filter_type = {'gabor_1d','morlet_1d'};
filt1_opt.Q = [8 1];
filt1_opt.J = T_to_J(512,filt1_opt);

sc1_opt.M = 1;

Wop = wavelet_factory_1d(N, filt1_opt, sc1_opt);

scatt_fun = @(x)(format_scat(log_scat(renorm_scat(scat(x,Wop)))));

duration_fun = @(x,obj)(32*duration_feature(x,obj));

features = {scatt_fun, duration_fun};

for k = 1:length(features)
	fprintf('testing feature #%d...',k);
	tic;
	if nargin(features{k}) == 1
		sz = size(features{k}(randn(N,1)));
	else
		sz = size(features{k}(randn(N,1),struct('u1',1,'u2',N)));
	end
	aa = toc;
	fprintf('OK (%.2fs) (size [%d,%d])\n',aa,sz(1),sz(2));
end

database_opt.input_sz = N;
database_opt.output_sz = T_s;
database_opt.obj_normalize = 2;
database_opt.collapse = 1;

db = prepare_database(src,features,database_opt);
db.features = single(db.features);
db = svm_calc_kernel(db,'gaussian','triangle',[db.indices{train_set}]);

optt.kernel_type = 'gaussian';
optt.gamma = 2.^[-14:2:-10];
optt.C = 2.^[2:2:6];
optt.search_depth = 2;
optt.full_test_kernel = 1;

[dev_err_grid,C_grid,gamma_grid] = ...
	svm_adaptive_param_search(db,train_set,valid_set,optt);

[dev_err,ind] = min(dev_err_grid{end});
C = C_grid{end}(ind);
gamma = gamma_grid{end}(ind);

optt1 = optt;
optt1.C = C;
optt1.gamma = gamma;

model = svm_train(db,train_set,optt1);
labels = svm_test(db,model,test_set);
err = classif_err(labels,test_set,db.src);

save([run_name '.mat'],'labels','err','C','gamma');

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function M = euclideanfactory(m, n)
% Returns a manifold struct to optimize over real matrices.
%
% function M = euclideanfactory(m)
% function M = euclideanfactory(m, n)
% function M = euclideanfactory([n1, n2, ...])
%
% Returns M, a structure describing the Euclidean space of real matrices,
% equipped with the standard Frobenius distance and associated trace inner
% product, as a manifold for Manopt.
%
% m and n in general can be vectors to handle multidimensional arrays.
% If either of m or n is a vector, they are concatenated as [m, n].
%
% Using this simple linear manifold, Manopt can be used to solve standard
% unconstrained optimization problems, for example in replacement of
% Matlab's fminunc.
%
% See also: euclideancomplexfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors: Bamdev Mishra, May 4, 2015.
% Change log:
%
%   July 5, 2013 (NB):
%       Added egred2rgrad, ehess2rhess, mat, vec, tangent.
%   May 4, 2015 (BM):
%       Added functionality to handle multidimensional arrays.

    % The size can be defined using both m and n, or simply with m.
    % If m is a scalar, then n is implicitly 1.
    % This mimics the use of built-in Matlab functions such as zeros(...).
    if ~exist('n', 'var') || isempty(n)
        if numel(m) == 1
            n = 1;
        else
            n = [];
        end
    end

    dimensions_vec = [m(:)', n(:)']; % We have a row vector.

    M.size = @() dimensions_vec;

    M.name = @() sprintf('Euclidean space R^(%s)', num2str(dimensions_vec));

    M.dim = @() prod(dimensions_vec);

    M.inner = @(x, d1, d2) d1(:).'*d2(:);

    M.norm = @(x, d) norm(d(:), 'fro');

    M.dist = @(x, y) norm(x(:) - y(:), 'fro');

    M.typicaldist = @() sqrt(prod(dimensions_vec));

    M.proj = @(x, d) d;

    M.egrad2rgrad = @(x, g) g;

    M.ehess2rhess = @(x, eg, eh, d) eh;

    M.tangent = M.proj;

    M.exp = @exp;
    function y = exp(x, d, t)
        if nargin == 3
            y = x + t*d;
        else
            y = x + d;
        end
    end

    M.retr = M.exp;

    M.log = @(x, y) y-x;

    M.hash = @(x) ['z' hashmd5(x(:))];

    M.rand = @() randn(dimensions_vec);

    M.randvec = @randvec;
    function u = randvec(x) %#ok<INUSD>
        u = randn(dimensions_vec);
        u = u / norm(u(:), 'fro');
    end

    M.lincomb = @matrixlincomb;

    M.zerovec = @(x) zeros(dimensions_vec);

    M.transp = @(x1, x2, d) d;
    M.isotransp = M.transp; % the transport is isometric

    M.pairmean = @(x1, x2) .5*(x1+x2);

    M.vec = @(x, u_mat) u_mat(:);
    M.mat = @(x, u_vec) reshape(u_vec, dimensions_vec);
    M.vecmatareisometries = @() true;
    M.lie_identity = @() zeros(dimensions_vec);

end

author: NicolasBoumal

repo: manopt


text:
function [SO3G, W, M2] = quadratureSO3Grid(bandwidth, varargin)
% Calculate the nodes and weights of specific quadrature grids on SO(3)
%   'ClenshawCurtis'  - combinate the Gauss quadrature in 1st and 3rd Euler angle
%                       (alpha, gamma) with Clenshaw Curtis quadrature in 2nd
%                       Euler angle (beta)
%
% Syntax
%   [SO3G, W, M2] = quadratureSO3rid(M, 'ClenshawCurtis') quadrature grid of type Clenshaw Curtis
%

% persistent SO3G_p;
% persistent W_p;
% persistent M2_p;

if check_option(varargin, 'gauss') || check_option(varargin, 'chebyshev')
  warning('Quadrature of types gauss and chebychev are not implemented yet.')
  return
end

% if check_option(varargin, 'ClenshawCurtis')
%   M2 = bandwidth;                           % bandwidth
%   % grid has to be full sized
%   N =  (2*bandwidth+2)^2*(2*bandwidth+1);   % number of nodes
% end
%
% if ~isempty(M2_p) && M2_p == M2 && length(SO3G_p) == N
%   SO3G = SO3G_p;
%   if check_option(varargin,'Euler') && isa(SO3G,'rotation')
%     SO3G = Euler(SO3G);
%   elseif ~check_option(varargin,'Euler') && isnumeric(SO3G)
%     [CS,SS] = extractSym(varargin);
%     SO3G = orientation.byEuler(SO3G(:,:,:,1),SO3G(:,:,:,2),SO3G(:,:,:,3),'nfft',CS,SS);
%   end
%   W = W_p;
%   M2 = M2_p;
%
% else
  if check_option(varargin,'ClenshawCurtis')
    SO3G = regularSO3Grid('ClenshawCurtis','bandwidth',bandwidth,varargin{:});

    w = fclencurt2(bandwidth+1);
    if size(SO3G,2) == 1+bandwidth/2
      w = w(1:size(SO3G,2));
    end
    % W = normalized volume * 2xGauß-quadrature-weights * Clenshaw-Curtis-quadrature-weights
    % W =     1/sqrt(8*pi^2)    *      (2*pi/(2*N+2))^2     *              w_b^N
    W = 1/(2*size(SO3G,1)*size(SO3G,3))*repmat(w',size(SO3G,1),1,size(SO3G,3));
  end
%   % store the data
%   SO3G_p = SO3G;
%   W_p = W;
%   M2_p = M2;
% end

end

function w = fclencurt2(n)

  c = zeros(n,2);
  c(1:2:n,1) = (2./[1 1-(2:2:n-1).^2 ])';
  c(2,2)=1;
  f = real(ifft([c(1:n,:);c(n-1:-1:2,:)]));
  w = 2*([f(1,1); 2*f(2:n-1,1); f(n,1)])/2;
  %x = 0.5*((b+a)+n*bma*f(1:n1,2));

end


text:
function fem = laplace_smooth(fem_struct,nodenums,iter);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LAPLACE_SMOOTH performs a Laplacian smoothing, aka "springing", on the
%  mesh. This is achieved by moving each non-boundary point to the average
%  location of the x,y (,z) coordinates of the nodes connected to the given
%  point. The process is repeated according to the 'iter' input parameter.
%  The smoothing may also be performed on a limited portion of the mesh
%  using the 'nodenums' input parameter.
%
% Usage -- fem = laplace_smooth(fem_struct,nodenums,iter);
%                ** nodenums & iter are optional **
%
% Variables
%  fem = updated finite element mesh
%  fem_struct = original finite element mesh
%  nodenums = node numbers of the nodes to be smoothed
%              (default is all nodes)
%  iter = number of iterations of smoothing to perform
%              (default is 5 iterations)
%
% Filename: laplace_smooth.m
% Created by: Ben Holladay
% Date: May 8, 2007
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parse input parameters. Throw exceptions for mal-formed input.
%
% Desired input variables: fem_struct, (nodenums), (iter)
%
% Defaults and valid flags
iterDefault = 5;
% nodenumsDefault = all non-boundary nodes
%
% Switch on number of input parameters.
switch nargin
   case 0
      error('Not enough input parameters: laplace_smooth() must have a fem_struct.');
   case 1
      iter = iterDefault;
      nodenums = 1:length(fem_struct.x);
   case 2
      iter = iterDefault;
      test = (floor(nodenums) == nodenums) & (nodenums > 0) & ...
         (nodenums <= length(fem_struct.x));
      if ~all(test)
         error('Input parameter ''nodenums'' must be a vector of valid node numbers.');
      end
   case 3
      iterValid = (floor(iter) == iter) && (iter > 0);
      if ~((floor(iter) == iter) && (iter > 0))
         error('Input parameter ''iter'' must be a positive integer');
      end
      test = (floor(nodenums) == nodenums) & (nodenums > 0) & ...
         (nodenums <= length(fem_struct.x));
      if ~all(test)
         error('Input parameter ''nodenums'' must be a vector of valid node numbers.');
      end
   otherwise
      error('Too many input parameters.');
end

% Set flag if their is no z variable.
try
   z = fem_struct.z;
   zflag = true;
catch
   zflag = false;
end

% Initialize basic fem_struct variables.
enodes = fem_struct.e;
x = fem_struct.x;
y = fem_struct.y;
bnd = unique(fem_struct.bnd);
np = size(x,1);

% Nodes to be moved.
moveable = setdiff(1:length(x),bnd);
move = intersect(moveable,nodenums);
conflict = setdiff(nodenums,moveable);
if ~isempty(conflict)
   display('Boundary nodes will not be moved.');
end

% Spring the number of times given by iter.
for k = 1:iter
   % Create sparse matrix of neighbor points x,y,&z coordinates.
   tempx = sparse(enodes,enodes(:,[2,3,1]),x(enodes),np,np);
   tempy = sparse(enodes,enodes(:,[2,3,1]),y(enodes),np,np);
   if zflag
      tempz = sparse(enodes,enodes(:,[2,3,1]),z(enodes),np,np);
   end

   % Calculate the number of neighbors.
   numnghb = sparse(enodes,enodes(:,[2,3,1]),1,np,np);
   numnghb = sum(numnghb);

   % Find the new location.
   tempx = sum(tempx) ./ numnghb;
   tempy = sum(tempy) ./ numnghb;

   % Move nodes.
   x(move) = tempx(move);
   y(move) = tempy(move);

   % Find new location move z coordinate.
   if zflag
      tempz = sum(tempz) ./ numnghb;
      z(move) = tempz(move);
   end
end

% Create output fem_struct.
fem = fem_struct;
fem.x = x;
fem.y = y;
if zflag
   fem.z = z;
end
if is_valid_struct(fem_struct)
    fem = el_areas(fem);
end

return

author: CHLNDDEV


text:
% Upgrades the network object to the most recent format.
% The motivation is the same as the motivation behind vl_simplenn_tidy
% (which is also executed during the call) as it enables NetVLAD code
% to evolve while not breaking backward compatibility

function net= relja_simplenn_tidy(net)

    % --- NetVLAD save values which will be deleted by MatConvNet tidy

    toSave= struct();

    % from v100-v101
    fieldNames= {'netID', 'sessionID', 'epoch'};

    for fieldName_= fieldNames
        fieldName= fieldName_{1};
        if isfield(net, fieldName)
            toSave.(fieldName)= net.(fieldName);
            net= rmfield(net, fieldName);
        end
    end

    % --- MatConvNet tidy

    net= vl_simplenn_tidy(net);

    % --- NetVLAD tidy

    for fieldName_= fieldNames
        fieldName= fieldName_{1};
        if isfield(toSave, fieldName)
            net.meta.(fieldName)= toSave.(fieldName);
            toSave= rmfield(toSave, fieldName);
        end
    end
    assert(isempty(fieldnames(toSave)));

end

author: Relja

repo: netvlad

save_path: github-repos/MATLAB/Relja-netvlad


text:
function duhatdt = rhsHeat(t,uhat,kappa,a)
duhatdt = -a^2*(kappa.^2)'.*uhat;  % Linear and diagonal

author: dynamicslab


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code is to replicate the results of RUSBoost model               %
% Our results are generated using Matlab R2020b on Windows 10           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

diary("results_rusboost.txt");
for year_test = 2003:2014
    rng(0,'twister'); % fix random seed for reproducing the results
    % read training data
    fprintf('==> Running RUSBoost (training period: %d-%d, testing period: %d, with %d-year gap)...\n',1991,year_test-2,year_test,2);
    data_train = data_reader('data_FraudDetection_JAR2020.csv','data_default',1991,year_test-2);
    y_train = data_train.labels;
    X_train = data_train.features;
    paaer_train = data_train.paaers;
    % read testing data
    data_test = data_reader('data_FraudDetection_JAR2020.csv','data_default',year_test,year_test);
    y_test = data_test.labels;
    X_test = data_test.features;
    paaer_test = unique(data_test.paaers(data_test.labels~=0));

    % handle serial frauds using PAAER
    y_train(ismember(paaer_train,paaer_test)) = 0;

    % train model
    t1 = tic;
    t = templateTree('MinLeafSize',5); % base model
    % see "tune_RUSBoost.m" for the tuning of the number of trees
    rusboost = fitensemble(X_train,y_train,'RUSBoost',300,t,'LearnRate',0.1,'RatioToSmallest',[1 1]);
    t_train = toc(t1);

    % test model
    t2 = tic;
    [label_predict,dec_values] = predict(rusboost,X_test);
    dec_values = dec_values(:,2);
    t_test = toc(t2);

    % print performance results
    fprintf('Training time: %g seconds | Testing time %g seconds \n', t_train, t_test);
    for topN = [0.01, 0.02, 0.03, 0.04, 0.05]
        metrics = evaluate(y_test,label_predict,dec_values,topN);
        fprintf('Performance (top%d%% as cut-off thresh): \n',topN*100);
        fprintf('AUC: %.4f \n', metrics.auc);
        fprintf('NCDG@k: %.4f \n', metrics.ndcg_at_k);
        fprintf('Sensitivity: %.2f%% \n', metrics.sensitivity_topk*100);
        fprintf('Precision: %.2f%% \n', metrics.precision_topk*100);
    end
end
diary off;

author: JarFraud

repo: FraudDetection


text:
function A = imresize3d_trans(V,scale,tsize,ntype,npad)
    A1 = imresize3d(V(:,:,:,1), scale, tsize, ntype, npad);
    A2 = imresize3d(V(:,:,:,2), scale, tsize, ntype, npad);
    A3 = imresize3d(V(:,:,:,3), scale, tsize, ntype, npad);
    A = cat(4, A1, A2, A3);
end

repo: pTVreg


text:
function r = zerofun(varargin)
r = 0;

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function hlpText = hlp_getFcnPreambleText(functionName)
% Return the 'preamble' text for a given function name.
% This is the initial part of the help text, and is terminated by a
% horizontal dashed line at least 5 characters long. The preamble text
% cannot contain any uncommented line breaks as these will be interpreted
% as indicating the end of the block of help text.
%
% Trailing whitespace is removed from the preamble.
%
% Author: Tim Mullen June 2012, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

% make sure function name is valid
try
    utl_whichfile(functionName);
catch
    error('SIFT:hlp_getFcnPreambleText:badFcnName','Unable to find function %s.m',functionName);
end

% Next, read in the help text
if isdeployed
    hlpText = hlp_readHelpInDeployedMode(functionName);
else
    hlpText = help(functionName);
end

if isempty(hlpText)
    error('Error parsing %s.m: I could not find a block of help text for the function',functionName);
end

% parse the help text. We only want the text preceeding the first
% sufficiently long horizontal dashed line
endpoint = strfind(hlpText,'-----');

if isempty(endpoint)
    disp(hlpText)
    % error('Error parsing %s.m: I could not find the horizontal dashed line separating the preamble from remaining help text',functionName);
    warning('Error parsing %s.m: I could not find the horizontal dashed line separating the preamble from remaining help text',functionName);
    hlpText = '';
end

% return the preamble
hlpText = deblank(hlpText(1:endpoint-1));

repo: REST


text:
function lab_nodes_local = extract(nodes_local, lab_begin, lab_end)

% Special routine for locating the nodes in node_local that correspond to
% columns on the lab. We rely on the fact that lab_cols are consecutive
% integers from lab_begin to lab_end (inclusive).

  lab_nodes_local = [];
  for kk = 1:length(nodes_local)
      if nodes_local(kk) >= lab_begin && nodes_local(kk) <= lab_end
          lab_nodes_local = [lab_nodes_local; kk nodes_local(kk)];
      end
  end

  return
end


text:
function [fea_train, train_targets, fea_test, test_chunk_id] = features(data, prediction_offset)

time_back = 8;

fea_train = zeros(40000, 3 + 89*time_back);
fea_test = zeros(500, 3 + 89*time_back);

train_targets = zeros(40000, 39);

test_chunk_id = [];

fea_cnt = 0;
test_cnt = 0;
for i=1:size(data,1)-time_back-prediction_offset+1
    if data(i,2)==data(i+time_back+prediction_offset-1,2)
        fea_cnt = fea_cnt + 1;
        fea_train(fea_cnt,1:3) = data(i, 4:6);
        this_fea = data(i:i+time_back-1,7:95);
        fea_train(fea_cnt,4:end) = this_fea(:)';

        train_targets(fea_cnt, :) = data(i+time_back+prediction_offset-1, 95-39+1:95);
    end

    if data(i,2) ~= data(i+1,2)
        test_cnt = test_cnt + 1;
        i_back = i - time_back + 1;
        fea_test(test_cnt,1:3) = data(i_back, 4:6);

        this_fea = data(i_back:i_back+time_back-1,7:95);
        fea_test(test_cnt,4:end) = this_fea(:)';
        test_chunk_id(end+1) = data(i_back,2);
    end
end

test_cnt = test_cnt + 1;
i_back = size(data,1) - time_back + 1;
fea_test(test_cnt,1:3) = data(i_back, 4:6);

this_fea = data(i_back:i_back+time_back-1,7:95);
fea_test(test_cnt,4:end) = this_fea(:)';
test_chunk_id(end+1) = data(i_back,2);

train_targets = train_targets(1:fea_cnt,:);
fea_train = fea_train(1:fea_cnt,:);
fea_test = fea_test(1:test_cnt, :);

author: benhamner


text:
function fe = spline_pchip_val ( n, x, f, d, ne, xe )

%*****************************************************************************80
%
%% SPLINE_PCHIP_VAL evaluates a piecewise cubic Hermite function.
%
%  Description:
%
%    This routine may be used by itself for Hermite interpolation, or as an
%    evaluator for SPLINE_PCHIP_SET.
%
%    This routine evaluates the cubic Hermite function at the points XE.
%
%    Most of the coding between the call to CHFEV and the end of
%    the IR loop could be eliminated if it were permissible to
%    assume that XE is ordered relative to X.
%
%    CHFEV does not assume that X1 is less than X2.  Thus, it would
%    be possible to write a version of SPLINE_PCHIP_VAL that assumes a strictly
%    decreasing X array by simply running the IR loop backwards
%    and reversing the order of appropriate tests.
%
%    The present code has a minor bug, which I have decided is not
%    worth the effort that would be required to fix it.
%    If XE contains points in [X(N-1),X(N)], followed by points less than
%    X(N-1), followed by points greater than X(N), the extrapolation points
%    will be counted (at least) twice in the total returned in IERR.
%
%    The evaluation will be most efficient if the elements of XE are
%    increasing relative to X; that is, for all J <= K,
%      X(I) <= XE(J)
%    implies
%      X(I) <= XE(K).
%
%    If any of the XE are outside the interval [X(1),X(N)],
%    values are extrapolated from the nearest extreme cubic,
%    and a warning error is returned.
%
%    This routine was originally named "PCHFE".
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2005
%
%  Author:
%
%    Original FORTRAN77 version by Fred Fritsch.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Fred Fritsch, Ralph Carlson,
%    Monotone Piecewise Cubic Interpolation,
%    SIAM Journal on Numerical Analysis,
%    Volume 17, Number 2, April 1980, pages 238-246.
%
%  Parameters:
%
%    Input, integer N, the number of data points.  N must be at least 2.
%
%    Input, real X(N), the strictly increasing independent
%    variable values.
%
%    Input, real F(N), the function values.
%
%    Input, real D(N), the derivative values.
%
%    Input, integer NE, the number of evaluation points.
%
%    Input, real XE(NE), points at which the function is to
%    be evaluated.
%
%    Output, real FE(NE), the values of the cubic Hermite
%    function at XE.
%

%
%  Check arguments.
%
  if ( n < 2 )
    ierr = -1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SPLINE_PCHIP_VAL - Fatal error!\n' );
    fprintf ( 1, '  Number of data points less than 2.\n' );
    error ( 'SPLINE_PCHIP_VAL - Fatal error!' );
  end

  for i = 2 : n
    if ( x(i) <= x(i-1) )
      ierr = -3;
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SPLINE_PCHIP_VAL - Fatal error!\n' );
      fprintf ( 1, '  X array not strictly increasing.\n' );
      error ( 'SPLINE_PCHIP_VAL - Fatal error!' );
    end
  end

  if ( ne < 1 )
    ierr = -4;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SPLINE_PCHIP_VAL - Fatal error!\n' );
    fprintf ( 1, '  Number of evaluation points less than 1.\n' );
    fe = [];
    return
  end

  ierr = 0;
%
%  Loop over intervals.
%  The interval index is IL = IR-1.
%  The interval is X(IL) <= X < X(IR).
%
  j_first = 1;
  ir = 2;

  while ( 1 )
%
%  Skip out of the loop if have processed all evaluation points.
%
    if ( ne < j_first )
      break
    end
%
%  Locate all points in the interval.
%
    j_save = ne + 1;

    for j = j_first : ne
      if ( x(ir) <= xe(j) )
        j_save = j;
        if ( ir == n )
          j_save = ne + 1;
        end
        break
      end
    end
%
%  Have located first point beyond interval.
%
    j = j_save;

    nj = j - j_first;
%
%  Skip evaluation if no points in interval.
%
    if ( nj ~= 0 )
%
%  Evaluate cubic at XE(J_FIRST:J-1).
%
      [ fe(j_first:j-1), next, ierc ] = chfev ( x(ir-1), x(ir), f(ir-1), ...
        f(ir), d(ir-1), d(ir),  nj, xe(j_first:j-1) );

      if ( ierc < 0 )
        ierr = -5;
        fprintf ( 1, '\n' );
        fprintf ( 1, 'SPLINE_PCHIP_VAL - Fatal error!\n' );
        fprintf ( 1, '  Error return from CHFEV.\n' );
        error ( 'SPLINE_PCHIP_VAL - Fatal error!' );
      end
%
%  In the current set of XE points, there are NEXT(2) to the right of X(IR).
%
      if ( next(2) ~= 0 )

        if ( ir < n )
          ierr = -5;
          fprintf ( 1, '\n' );
          fprintf ( 1, 'SPLINE_PCHIP_VAL - Fatal error!\n' );
          fprintf ( 1, '  IR < N.\n' );
          error ( 'SPLINE_PCHIP_VAL - Fatal error!' );
        end
%
%  These are actually extrapolation points.
%
        ierr = ierr + next(2);

      end
%
%  In the current set of XE points, there are NEXT(1) to the left of X(IR-1).
%
      if ( next(1) ~= 0 )
%
%  These are actually extrapolation points.
%
        if ( ir <= 2 )
          ierr = ierr + next(1);
        else

          j_new = -1;

          for i = j_first : j-1
            if ( xe(i) < x(ir-1) )
              j_new = i;
              break
            end
          end

          if ( j_new == -1 )
            ierr = -5;
            fprintf ( 1, '\n' );
            fprintf ( 1, 'SPLINE_PCHIP_VAL - Fatal error!\n' );
            fprintf ( 1, '  Could not bracket the data point.\n' );
            error ( 'SPLINE_PCHIP_VAL - Fatal error!' );
          end
%
%  Reset J.  This will be the new J_FIRST.
%
          j = j_new;
%
%  Now find out how far to back up in the X array.
%
          for i = 1 : ir-1
            if ( xe(j) < x(i) )
              break
            end
          end
%
%  At this point, either XE(J) < X(1) or X(i-1) <= XE(J) < X(I) .
%
%  Reset IR, recognizing that it will be incremented before cycling.
%
          ir = max ( 1, i-1 );

        end

      end

      j_first = j;

    end

    ir = ir + 1;

    if ( n < ir )
      break
    end

  end

  return
end


text:
function ttt=im2col_inds(a, block)
%IM2COL Rearrange image blocks into columns.
%   B = IM2COL(A,[M N],'distinct') rearranges each distinct
%   M-by-N block in the image A into a column of B. IM2COL pads A
%   with zeros, if necessary, so its size is an integer multiple
%   of M-by-N. If A = [A11 A12; A21 A22], where each Aij is
%   M-by-N, then B = [A11(:) A21(:) A12(:) A22(:)].
%
%   B = IM2COL(A,[M N],'sliding') converts each sliding M-by-N
%   block of A into a column of B, with no zero padding. B has
%   M*N rows and will contain as many columns as there are M-by-N
%   neighborhoods in A. If the size of A is [MM NN], then the
%   size of B is (M*N)-by-((MM-M+1)*(NN-N+1). Each column of B
%   contains the neighborhoods of A reshaped as NHOOD(:), where
%   NHOOD is a matrix containing an M-by-N neighborhood of
%   A. IM2COL orders the columns of B so that they can be
%   reshaped to form a matrix in the normal way. For example,
%   suppose you use a function, such as SUM(B), that returns a
%   scalar for each column of B. You can directly store the
%   result in a matrix of size (MM-M+1)-by-(NN-N+1) using these
%   calls:
%
%        B = im2col(A,[M N],'sliding');
%        C = reshape(sum(B),MM-M+1,NN-N+1);
%
%   B = IM2COL(A,[M N]) uses the default block type of
%   'sliding'.
%
%   B = IM2COL(A,'indexed',...) processes A as an indexed image,
%   padding with zeros if the class of A is uint8 or uint16, or
%   ones if the class of A is double.
%
%   Class Support
%   -------------
%   The input image A can be numeric or logical. The output matrix
%   B is of the same class as the input image.
%
%   Example
%   -------
%   Calculate the local mean using a [2 2] neighborhood with zero padding.
%
%       A = reshape(linspace(0,1,16),[4 4])'
%       B = im2col(A,[2 2])
%       M = mean(B)
%       newA = col2im(M,[1 1],[3 3])
%
%   See also BLOCKPROC, COL2IM, COLFILT, NLFILTER.

%   Copyright 1993-2016 The MathWorks, Inc.

[ma,na] = size(a);
m = block(1); n = block(2);

if any([ma na] < [m n]) % if neighborhood is larger than image
    b = zeros(m*n,0);
    return
end

% Create Hankel-like indexing sub matrix.
mc = block(1); nc = ma-m+1; nn = na-n+1;
cidx = (0:mc-1)'; ridx = 1:nc;
t = cidx(:,ones(nc,1)) + ridx(ones(mc,1),:);    % Hankel Subscripts
tt = zeros(mc*n,nc);
rows = 1:mc;
for i=0:n-1,
    tt(i*mc+rows,:) = t+ma*i;
end
ttt = zeros(mc*n,nc*nn);
cols = 1:nc;
for j=0:nn-1,
    ttt(:,j*nc+cols) = tt+ma*j;
end

%%%
%%% Function parse_inputs
%%%
function [a, block, kind, padval] = parse_inputs(varargin)

narginchk(2,4);

switch nargin
    case 2
        if (strcmp(varargin{2},'indexed'))
            error(message('images:im2col:tooFewInputs'))
        else
            % IM2COL(A, [M N])
            a = varargin{1};
            block = varargin{2};
            kind = 'sliding';
            padval = 0;
        end

    case 3
        if (strcmp(varargin{2},'indexed'))
            % IM2COL(A, 'indexed', [M N])
            a = varargin{1};
            block = varargin{3};
            kind = 'sliding';
            padval = 1;
        else
            % IM2COL(A, [M N], 'kind')
            a = varargin{1};
            block = varargin{2};
            kind = validatestring(varargin{3},{'sliding','distinct'},mfilename,'kind',3);
            padval = 0;
        end

    case 4
        % IM2COL(A, 'indexed', [M N], 'kind')
        a = varargin{1};
        block = varargin{3};
        kind = validatestring(varargin{4},{'sliding','distinct'},mfilename,'kind',4);
        padval = 1;

end

if (isa(a,'uint8') || isa(a, 'uint16'))
    padval = 0;
end

author: TadasBaltrusaitis

repo: OpenFace


text:
%isivector - Test if parameter is a vector of integers satisfying an optional list of tests.
%
%  USAGE
%
%    test = isivector(x,test1,test2,...)
%
%    x              parameter to test
%    test1...       optional list of additional tests (see examples below)
%
%  EXAMPLES
%
%    % Test if x is a vector of doubles
%    isivector(x)
%
%    % Test if x is a vector of strictly positive doubles
%    isivector(x,'>0')
%
%    % Test if x is a vector of doubles included in [2,3]
%    isivector(x,'>=2','<=3')
%
%    % Special test: test if x is a vector of doubles of length 3
%    isivector(x,'#3')
%
%    % Special test: test if x is a vector of strictly ordered doubles
%    isivector(x,'>')
%
%  NOTE
%
%    The tests ignore NaNs, e.g. isivector([500 nan]), isivector([1 nan 3],'>0') and
%    isivector([nan -7],'<=0') all return 1.
%
%  SEE ALSO
%
%    See also isdmatrix, isdvector, isdscalar, isimatrix, isiscalar, isstring,
%    islscalar, islvector, islmatrix.
%

% Copyright (C) 2010 by Michaël Zugaro
%
% 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.

function test = isivector(x,varargin)

% Check number of parameters
if nargin < 1,
  error('Incorrect number of parameters (type ''help isivector'' for details).');
end

% Test: double, vector
test = isa(x,'double') & isvector(x);

% Ignore NaNs
x = x(~isnan(x));

% Test: integers?
test = test & all(round(x)==x);

% Optional tests
for i = 1:length(varargin),
	try
		if varargin{i}(1) == '#',
			if length(x) ~= str2num(varargin{i}(2:end)), test = false; return; end
		elseif isstring(varargin{i},'>','>=','<','<='),
			dx = diff(x);
			if ~eval(['all(0' varargin{i} 'dx);']), test = false; return; end
		else
			if ~eval(['all(x' varargin{i} ');']), test = false; return; end
		end
	catch err
		error(['Incorrect test ''' varargin{i} ''' (type ''help isivector'' for details).']);
	end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function value = i8_btest ( i, pos )

%*****************************************************************************80
%
%% I8_BTEST returns TRUE if bit #POS of I is 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    MIL-STD 1753
%
%  Parameters:
%
%    Input, integer I, the integer to be tested.
%
%    Input, integer POS, the bit position, between 0 and 63.
%
%    Output, logical VALUE, is TRUE if the POS-th bit of I is 1.
%
  if ( pos < 0 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'I8_BTEST - Fatal error!\n' );
    fprintf ( 1, '  POS < 0.\n' );
    error ( 'I8_BTEST - Fatal error!' );

  elseif ( pos < 63 )

    j = i;
    for k = 1 : pos
      j = floor ( j / 2 );
    end

    if ( mod ( j, 2 ) == 0 )
      value = 0;
    else
      value = 1;
    end

  elseif ( pos == 63 )

    if ( i < 0 )
      value = 1;
    else
      value = 0;
    end

  elseif ( 63 < pos )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'I8_BTEST - Fatal error!\n' );
    fprintf ( 1, '  63 < POS.\n' );
    error ( 'I8_BTEST - Fatal error!' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% ENVPROJ - plot envelopes of projections of selected ICA component
%             projections against envelope of the original data
%
% Usage:    >> [envdata] = envproj(data,weights,compnums);
%           >> [envdata] = envproj(data,weights,compnums, ...
%                                  title,limits,chanlist,compnames,colors);
% Inputs:
%   data      = RUNICA input data (chans,frames) <- best one epoch only!
%   weights   = unmixing weight matrix (RUNICA weights*sphere)
%   compnums  = list of component numbers to project and plot
%
% Optional inputs:
%   title     = 'fairly short plot title' (in quotes) (0 -> none)
%   limits    = [xmin xmax ymin ymax] (x's in msec)
%                       (0 or both y's 0 -> [min max])
%   chanlist  = list of data channels to use for max/min (0 -> all)
%   compnames = file of component labels (exactly 5 chars per line,
%               trailing '.' = space) (default|0 -> compnums)
%   colors    = file of color codes, 3 chars per line  ('.' = space)
%                       (0 -> default color order (black/white first))
%
% Output:
%   envdata   = [2,(frames*length(compnums)+1)] envelopes of data, comps.
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 1/1997
%
% See also: ENVTOPO

% Copyright (C) 01-23-97 Scott Makeig, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 03-19-97 use datamean instead of frames/baseframes, diag(COV) for VAR -sm
% 04-03-97 changed name to ENVPROJ -sm
% 05-20-97 used sum-of-squares for var instead of DIAG to allow long data
%          read MAXPLOTDATACHANS from icadefs.m -sm
% 06-07-97 changed order of args to conform to runica -sm
% 06-18-97 read MAXENVPLOTCHANS instead of MAXPLOTDATACHANS -sm
% 07-23-97 dropped datamean from args; mean is distributed among components -sm
% 10-11-97 range of max/min for default ylimits over chanlist only! -sm
% 11-05-97 disallowed white lines unless axis color is white -sm & ch
% 11-13-97 rm'd errorcode variable to improve limits handling -sm
% 11-18-97 added legend to plot; varied line types -sm
% 11-30-97 added f=figure below to keep previous plot from being altered -sm
% 12-01-97 added thicker lines for data envelope -sm
% 12-08-97 debugged thicker lines for data envelope -sm
% 12-10-97 changed f=figure below to f=gcf to avoid matlab bug (gcf confusion) -sm
% 12-10-97 implemented compnames arg -sm
% 12-13-97 took out compnames=compnames'; bug -sm
% 02-06-98 added 'show-largest' feature to legend when numcomps > 7 -sm
% 02-09-98 fixed numcomps bug -sm
% 02-12-98 test for outofbounds compnums early -sm
% 04-28-98 made to work for rectangular weights -sm
% 06-05-98 made NEG_UP optional -sm
% 02-22-99 added FILL mode default when numcomps==1 -sm
% 03-14-99 made envelope finding code more efficient -sm
% 12-19-00 updated ICAPROJ args -sm
% 01-25-02 reformated help & license, added links -ad

function [envdata] = envproj(data,weights,compnums,titl,limits,chanlist,compnames,colors)

FILL = 1; % 1;  % 1 = fill in component envelope unless ncomps>1
DATA_LINEWIDTH = 1.5; % 2
COMP_LINEWIDTH = 1.5; % 1
MAXENVPLOTCHANS = 256;

if nargin < 3,
    help envproj
    fprintf('envproj(): must have at least three arguments.\n');
    return
end

icadefs        % load ENVCOLORS & MAXENVPLOTCHANS default variables
MAX_LEG = 7;   % maximum number of component line types to label in legend
NEG_UP = 0;    % 1 = plot negative up

%
%%%%%%%%%%%% Substitute for omitted arguments %%%%%%%%%%%%%%%%%%%%%%%%
%
if nargin < 8,
    colors = 'envproj.col';
elseif colors(1)==0,
    colors = 'envproj.col';
end

[chans,frames] = size(data);

if nargin < 7
   compnames = 0;
end

if nargin < 6
    chanlist = 0;
end
if nargin < 5,
    limits = 0;
end
if nargin < 4,
    titl = 0;
end
%
%%%%%%%%%%%%%%%%%%%%%%%% Test data size %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[wr,wc]       = size(weights);
%
%%%%%%%%%%%%%%%% Substitute for zero arguments %%%%%%%%%%%%%%%%%%%%%%%%%
%
if chanlist == 0,
    chanlist = [1:chans];
end
if compnums == 0,
    compnums = [1:wr];
end
if size(compnums,1)>1,        % handle column of compnums !
    compnums = compnums';
end
numcomps = length(compnums);
if numcomps > MAXENVPLOTCHANS,
    fprintf(...
'envproj(): cannot plot more than %d channels of data at once.\n',...
                  MAXENVPLOTCHANS);
    return
end

if max(compnums)>wr
 fprintf('\nenvproj(): Component index %d out of bounds (1:%d).\n',...
                                             max(compnums),wr);
 return
end

if min(compnums)<1,
 fprintf('\nenvproj(): Component index %d out of bounds (1:%d).\n',...
                                             min(compnums),wr);
 return
end

if FILL>0 && numcomps == 1
   FILL = 1;
else
   FILL = 0;
end
%
%%%%%%%%%%%%%%%%%%%% Read and adjust limits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
  if limits==0,
    xmin=0;xmax=1;
    ymin=min(min(data(chanlist,:)));
    ymax=max(max(data(chanlist,:)));
  else
    if length(limits)~=4,
        fprintf( ...
 'envproj(): limits should be 0 or an array [xmin xmax ymin ymax].\n');
        return
      end
        if limits(1,1) == 0 && limits(1,2) ==0,
            xmin=0;
            xmax=0;
        else
            xmin = limits(1,1);
            xmax = limits(1,2);
        end
         if limits(1,3) == 0 && limits(1,4) ==0,
            ymin=0;
            ymax=0;
        else
            ymin = limits(1,3);
            ymax = limits(1,4);
        end
  end

  if xmax == 0 && xmin == 0,
    x = (0:1:frames-1);
    xmin = 0;
    xmax = frames-1;
  else
    dx = (xmax-xmin)/(frames-1);
    x=xmin*ones(1,frames)+dx*(0:frames-1); % construct x-values
  end

  if ymax == 0 && ymin == 0,
    ymax=max(max(data(chanlist,:)));
    ymin=min(min(data(chanlist,:)));
  end

%
%%%%%%%%%%%%%%%%%%%% Read the color names %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    if ~ischar(colors)
        fprintf('envproj(): color file name must be a string.\n');
        return
    end
    cid = fopen(colors,'r');
    if cid <3,
        fprintf('envproj(): cannot open file %s.\n',colors);
        return
    else
        colors = fscanf(cid,'%s',[3 MAXENVPLOTCHANS]);
        colors = colors';
        [r c] = size(colors);
        for i=1:r
            for j=1:c
                if colors(i,j)=='.',
                    colors(i,j)=' ';
                end
            end
        end
    end
[rr cc] = size(colors);

%
%%%%%%%%%%%% Compute projected data for single components %%%%%%%%%%%%%
%
projdata = [data(chanlist,:)];
envdata = [min(projdata); max(projdata)]; % begin with envelope of data
fprintf('envproj(): Projecting component(s) ');
maxvars = zeros(1,length(compnums));
n=1;
for c=compnums,            % for each component
  fprintf('%d ',c);
  % [icaproj] = icaproj(data,weights,compindex); % new arg order 12/00
  proj = icaproj(data,weights,c); % let offsets be
                                         % distributed among comps.
  [val,i] = max(sum(proj.*proj)); % find max variance
  maxvars(n) = val;
  proj = proj(chanlist,:);
  projdata = [projdata proj];
  envtmp = [ min(proj) ; max(proj)];
  envdata = [envdata envtmp];
                 % append envelope of projected data sets onto envdata
  % Note: size(envdata) = [length(chanlist)  frames*(numcomps+1)]
  n = n+1;
end
fprintf('\n');
%
%%%%%%%%%%%%%%%%%%%%%%%% Make the plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if ~ischar(titl)
  if titl==0,
    titl = ' ';
  else
    fprintf('envproj(): titl unrecognized.\n');
    return
  end
end

set(gcf,'Color',BACKCOLOR); % set the background color
set(gca,'Color','none'); % set the axis color to the background color
if NEG_UP
  set(gca,'Ydir','reverse');
end
n = 1;
H = zeros(1,numcomps+1);

for c=1:numcomps+1;
  %
  % Plot envelopes
  %
  if n <= 6
    if n == 1              % thick lines, store H for legend()
      H(n)=plot(x,envdata(1,(c-1)*frames+1:c*frames)',...
                                     colors(c),'Linewidth',DATA_LINEWIDTH);
      hold on;
           plot(x,envdata(2,(c-1)*frames+1:c*frames)',...
                                     colors(c),'Linewidth',DATA_LINEWIDTH);
    else
       if FILL
            fillx = [x x(frames:-1:1)];
            filly = [envdata(1,(c-1)*frames+1:c*frames) ...
                            envdata(2,c*frames:-1:(c-1)*frames+1)];
            i = find(isnan(filly)==1);
            if ~isempty(i)
                fprintf('Replacing %d NaNs in envelope with 0s.\n',length(i)/2);
                filly(i)=0;
            end
            H(n) = fill(fillx,filly,colors(c));
       else
            H(n)= plot(x,envdata(1,(c-1)*frames+1:c*frames)',...
                                     colors(c),'Linewidth',COMP_LINEWIDTH);
             hold on;
            plot(x,envdata(2,(c-1)*frames+1:c*frames)',...
                                     colors(c),'Linewidth',COMP_LINEWIDTH);
       end
    end
  else
    H(n)= plot(x,envdata(1,(c-1)*frames+1:c*frames)',...
                                     colors(c),'LineStyle',':','LineWidth',DATA_LINEWIDTH);
    hold on;
          plot(x,envdata(2,(c-1)*frames+1:c*frames)',...
                                     colors(c),'LineStyle',':','LineWidth',DATA_LINEWIDTH);
  end
  n = n+1;
end
set(gca,'Color','none'); % set the axis color to the background color

l=xlabel('Time (ms)');      % xaxis label
set(l,'FontSize',14);

l=ylabel('Potential (uV)'); % yaxis label
set(l,'FontSize',14);

if xmax > 0  && xmin < 0
  plot([0 0],[ymin ymax],'k','linewidth',1.5); % plot vertical line at time 0
end

if ymax ~= 0 |ymin~= 0,
    axis([min(x),max(x),ymin,ymax]);
else
    axis([min(x),max(x),min(envdata(1,:)'),max(envdata(2,:)')]);
end
legtext = [ 'Data ' ];
if compnames(1) ~= 0
    pid = fopen(compnames,'r');
    if pid < 3
       fprintf('envproj(): file %s not found.\n',compnames);
       return
    end
    compnames = fscanf(pid,'%s',[5 numcomps]);
    compnames = [ ['Data.']' compnames];
    compnames = compnames';
    if size(compnames,1) ~= numcomps+1
      fprintf('envproj(): no. of compnames must equal no. of compnums.\n');
      return
    end
    for c=1:5 % remove padding with .'s
      for r = 1:size(compnames,1)
         if compnames(r,c) == '.'
           compnames(r,c) = ' ';
         end
      end
    end
    legtext = compnames;
else
  legtext = ['Data'];
  for c = compnums
      legtext = strvcat(legtext,int2str(c));
  end
end

if numcomps<MAX_LEG+1
  if FILL <= 0  % no legend in FILL mode
       l=legend(H,legtext,0);      % plot legend
  end
else
  fprintf('Finding largest components for legend...\n');
  [maxvars,i] = sort(maxvars);
  maxvars = maxvars(length(maxvars):-1:1); % sort large-to-small
  i= i(length(i):-1:1);
  H = [H(1) H(i+1)];
  H = H(1:MAX_LEG+1);
  legtext = [legtext(1,:);legtext(i+1,:)];
  legtext = legtext(1:MAX_LEG+1,:);
  legend(H,legtext,0); % MAX_LEG largest-peaked comps
end

set(gca,'FontSize',12);
fprintf('If desired, use mouse to move legend or >> delete(legend).\n');

%
%%%%%%%%%%%%% Compute percentage of variance accounted for %%%%%%%%%%%%
%
sumdata = zeros(length(chanlist),frames);
for e=2:numcomps+1,              % sum the component activities
    sumdata = sumdata + projdata(:,(e-1)*frames+1:e*frames);
end

sigssqr = sum(data(chanlist,:).*data(chanlist,:))/(length(chanlist)-1);
dif    = data(chanlist,:)-sumdata;
difssqr = sum(dif.*dif)/(length(chanlist)-1);
pvaf   = round(100.0*(1.0-difssqr/sigssqr)); % % variance accounted for
rtitl = ['(' int2str(pvaf) '%)'];
%titl = [titl '   ' rtitl];

t=title(titl);              % plot title
set(t,'FontSize',14);

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [z,n] = phiflower(p)
%% PHIFLOWER level set function of a flower
%
% Author: Huayi Wei <weihuayi@xtu.edu.cn>

petal = 5;
xc = 0;
yc = 0;
r0 = 0.5;
r1 = 0.2;
theta = atan2((p(:,2) - yc),(p(:,1) - xc));
r = r0 + r1*sin(petal*theta);

z = (p(:,2)-yc).^2 + (p(:,1)-xc).^2 - r.^2;

if nargout == 2
    rp = petal*r1*cos(petal*theta);
    n1 = rp.*cos(theta) - r.*sin(theta);
    n2 = rp.*sin(theta) + r.*cos(theta);
    l = sqrt(n1.^2+n2.^2);
    n = [n2./l, -n1./l];
end
end

repo: ifem


text:
import matlab.unittest.TestSuite
import matlab.unittest.TestRunner
import matlab.unittest.plugins.CodeCoveragePlugin

clear; clc;
suite = TestSuite.fromFile('replicationTests.m');
runner = TestRunner.withTextOutput;
runner.addPlugin(CodeCoveragePlugin.forFolder(bearroot(), 'IncludingSubfolders', true))
result = runner.run(suite);
tb = table(result);
disp(tb);

author: european-central-bank

repo: BEAR-toolbox


text:
function a = cosh(a)
%COSH         Gradient hyperbolic cosine cosh(a)
%

% written  10/16/98     S.M. Rump
% modified 10/14/00     S.M. Rump  use Tony's trick
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 10/08/08     S.M. Rump  improved sparse multiplication: not using intval data type
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  N = getappdata(0,'INTLAB_GRADIENT_NUMVAR');

  ax = sinh(a.x(:));
  a.x = cosh(a.x);
  if issparse(a.dx)
    sizeax = size(a.dx,1);
    [ia,ja,sa] = find(a.dx);
    % take care of scalar a.x: cures Matlab V6.0 bug
    % a=7; i=[1 1]; x=a(i), b=sparse(a); y=b(i)  yields row vector x but column vector y
    ax = ax(ia);
    if isa(a.x,'intval')
      adx = times(ax(:),sa(:),0);
      if adx.complex
        a.dx = intval( sparse(ia,ja,adx.mid,sizeax,N) , sparse(ia,ja,adx.rad,sizeax,N) , 'midrad' );
      else
        a.dx = intval( sparse(ia,ja,adx.inf,sizeax,N) , sparse(ia,ja,adx.sup,sizeax,N) , 'infsup' );
      end
    else
      a.dx = sparse(ia,ja,ax(:).*sa(:),sizeax,N);
    end
  else
    a.dx = a.dx .* ax(:,ones(1,N));
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function rtk=update_stat(rtk,obs,sol_stat)

global glc
opt=rtk.opt;
nobs=size(obs,1);

rtk.sol.time=timeadd(obs(1).time,-rtk.x(rtk.ic+1)/glc.CLIGHT);
rtk.sol.ns=0;
for i=1:nobs
    for j=1:opt.nf
        sat=obs(i).sat;
        if rtk.sat(sat).vsat(j)==0,continue;end
        rtk.sat(sat).lock(j)=rtk.sat(sat).lock(j)+1;
        rtk.sat(sat).outc(j)=0;
        if j==1,rtk.sol.ns=rtk.sol.ns+1;end
    end
end

if rtk.sol.ns<4
    rtk.sol.stat=0;
else
    rtk.sol.stat=sol_stat;
end

if rtk.sol.stat==glc.SOLQ_FIX

else
    rtk.sol.pos=rtk.x(1:3)';
    rtk.sol.posP(1)=rtk.P(1,1);
    rtk.sol.posP(2)=rtk.P(2,2);
    rtk.sol.posP(3)=rtk.P(3,3);
    rtk.sol.posP(4)=rtk.P(1,2);
    rtk.sol.posP(5)=rtk.P(2,3);
    rtk.sol.posP(6)=rtk.P(1,3);
end

% clk
rtk.sol.dtr(1)=rtk.x(rtk.ic+1)/glc.CLIGHT;
% isb
rtk.sol.dtr(2)=rtk.x(rtk.ic+2)/glc.CLIGHT;
rtk.sol.dtr(3)=rtk.x(rtk.ic+3)/glc.CLIGHT;
rtk.sol.dtr(4)=rtk.x(rtk.ic+4)/glc.CLIGHT;
rtk.sol.dtr(5)=rtk.x(rtk.ic+5)/glc.CLIGHT;

for i=1:nobs
    for j=1:rtk.opt.nf
        rtk.sat(obs(i).sat).snr(j)=obs(i).S(j);
    end
end

for i=1:glc.MAXSAT
    for j=1:rtk.opt.nf
        if bitand(rtk.sat(i).slip(j),3)
            rtk.sat(i).slip(j)=rtk.sat(sat).slip(j)+1;
        end
        if rtk.sat(i).fix(j)==2&&sol_stat~=glc.SOLQ_FIX
            rtk.sat(i).fix(j)=1;
        end
    end
end

return

repo: GINav