text:
function DATA = read_edat_output_2008(fname, varargin)
% Function that creates a structure DATA containing columns
% of the edat file output (saved in text tab delimited "excel" format)
%
% :Usage:
% ::
%
%     DATA = read_edat_output_2008(fname, varargin)
%
% For this code to work on a Mac, you must: 1) export .edat2 file as an Excel file,
% then 2) open this file in Excel on a Mac and save as a .csv, 3) read that
% .csv file
%
% :Examples:
% ::
%
%    fname = 'myfile.txt';
%    DATA = read_edat_output_2008(fname)
%
% :Defaults:
% These are the default formats this function expects:
% tab delimited, 1 header row, then row of column names, then data
%
% You can override some of them by using the following --
% E.g., for zero header rows and comma delimited data:
% ::
%
%    DATA = read_edat_output_2008(fname, 'nheaderrows', 0, 'mydelimiter', ',')
%
% You can force the number of columns to be a certain value by doing the
% following:
%% ::
%
%    DATA = read_edat_output_2008(fname, 'nheaderrows', 1, 'numc', 103);
%
% This could be useful if your last row contains empty cells at the end,
% which will mess up the automatic calculation of number of columns.
%
% ..
%    Tor Wager, Oct 2008
% ..

% ..
%    read the database initially to get all column names, etc.
% ..
DATA = [];

% assume you have one header row
nheaderrows = 1;
mydelimiter = '\t';
numc = [];

for i = 1:length(varargin)
        if ischar(varargin{i})
            switch varargin{i}
                % reserved keywords
                case 'nheaderrows', nheaderrows = varargin{i+1};
                case 'mydelimiter', mydelimiter = varargin{i+1};

                case 'numc', numc = varargin{i+1};

                otherwise, warning(['Unknown input string option:' varargin{i}]);
            end
        end
end

col1 = textread(fname, '%s%*[^\n]', 'delimiter', mydelimiter, 'headerlines', nheaderrows, 'bufsize', 10000000);

nrows = length(col1);

%numc = 148;
[d] = textread(fname,'%s','headerlines', nheaderrows, 'delimiter', mydelimiter, 'bufsize', 10000000);

for i = 1:length(d), if isempty(d{i}), d{i} = 'NaN'; end, end

if isempty(numc)
    % determine from data
    numc = length(d) ./ nrows;
else
    % force numc to be input value
    disp(['Forcing number of columns to be ' num2str(numc)])
    totalel = nrows * numc;
    if length(d) < totalel
        numtoadd = totalel - length(d);
        for i = 1:numtoadd
            d{end+1} = 'NaN';
        end
    else
        d = d(1:totalel);
    end
end

% % if strcmp(d{1},'This file contained edited data.'),
% %     % add an extra input row!
% %     disp('edited data: extra row.')
% %     addone = 1;
% %     [d] = textread(fname,'%s','headerlines',2,'delimiter','\t');
% % end

if mod(length(d), nrows) == 0
    fprintf('This data matrix appears to be OK, with %3.0f rows (including col. names) x %3.0f cols\n', nrows, numc)
else
    fprintf('This data matrix IS NOT OK\n')
    fprintf('i think it has %3.0f header rows, %3.2f data rows (including col. names) and %3.2f cols\n', nheaderrows, nrows, numc)

    disp('Here are the first 10 rows of the first column: ')
    disp(char(col1{1:10}))

    disp('One potential cause is if you have empty spaces in the last row in your file.')
    disp('then matlab will think it''s reached the end of the file and you will have too few rows.')
    disp('Try entering ''numc'' as an input with the number of columns, or delete the last row.')

    return
end

names = d(1:numc)'; disp(['Last name column read is ' names{end}])

disp(' ')
disp('Your column names:')
fprintf('%s ', names{:});
disp(' ')

% ------------------------------------------------------------------------
% replace illegal characters that will cause Matlab problems
% ------------------------------------------------------------------------

for i = 1:length(names), names{i}(findstr('.',names{i})) = '_'; end
for i = 1:length(names), names{i}(findstr('[',names{i})) = ''; end
for i = 1:length(names), names{i}(findstr(']',names{i})) = ''; end
for i = 1:length(names), names{i}(findstr(':',names{i})) = ''; end
for i = 1:length(names), names{i}(findstr('=',names{i})) = ''; end

% ------------------------------------------------------------------------
% Get data matrix in appropriate format
% ------------------------------------------------------------------------

% This only works if there are no blanks
%dmat = reshape(d, numc, nrows)';

% The code below works with blanks

% dmat(1, :) = [];  % get rid of names

% adjust for header rows
nrows = nrows - 1;

fmtstring = '%s%*[^\n]';

for i = 1:numc
    %mycol = dmat(:, i);

    % read in the nth column, skipping header rows and column names
    mycol = textread(fname, fmtstring, 'delimiter', mydelimiter, 'headerlines', nheaderrows + 1, 'bufsize', 100000);

    fmtstring = ['%*s' fmtstring];  % add this so we skip this col next time

    % Get the first non-empty entry and determine if it's numeric or text
    isvalidnumber(i) = 0;
    for j = 1:length(mycol)
        empt = isempty(mycol{j});

        if ~empt
            isvalidnumber(i) = ~isempty(str2num(mycol{j}));
            break
        end
    end

    if isvalidnumber(i)
        % This column has numbers
        for j = 1:nrows
            try
                mynum = str2num(mycol{j});  % sometimes, with 17:40:12 notation, returns more than 1 number
                if isempty(mynum), mynum = NaN; end

            catch
                disp('We seem to be mixing up text and numbers in columns.  Are your rows and cols correct?')
                disp('Here''s the current column (1st 15 rows)')
                mycol(1:15)

                disp('Here''s the first 10 rows and cols:')
                dmat(1:10, 1:10)

                keyboard
            end

            DATA.(names{i})(j, 1) = mynum(1);
        end
    else
        % This column is text
        DATA.(names{i}) = mycol;
    end

end

end

% %
% %
% %
% % % ------------------------------------------------------------------------
% % % create formatting string with number of columns and output var names
% % % ------------------------------------------------------------------------
% %
% % fmt = repmat('%s',1,numc);
% % fmt = [fmt '%*[^\n]'];
% %
% % outs = ['[' names{1}];
% % for i = 2:numc, outs = [outs ',' names{i}]; end
% % outs = [outs ']'];
% %
% % % works, but requires an extra step
% % %outs = ['[out{1}'];
% % %for i = 2:numc, outs = [outs ',out{' num2str(i) '}'];, end
% % %outs = [outs ']'];
% %
% % % ------------------------------------------------------------------------
% % % read the database again with proper formatting and output names
% % % ------------------------------------------------------------------------
% %
% % if addone
% %     str = [outs ' = textread(''' fname ''',fmt,''headerlines'',3,''delimiter'',''\t'');'];
% % else
% %     str = [outs ' = textread(''' fname ''',fmt,''headerlines'',2,''delimiter'',''\t'');'];
% % end
% %
% % eval(str)
% %
% % warning off
% % ww = whos('*RT*'); ww = ww(end).name; eval(['ww = isempty(str2num(' ww '{1}));'])
% % warning on
% % if ww,
% %     str = [outs ' = textread(''' fname ''',fmt,''headerlines'',3,''delimiter'',''\t'');'];
% %     eval(str)
% % end
% %
% % % not necessary if all filenames are OK
% % %for i = 1:length(outs),
% % %    eval([names{i} ' = out{i};'])
% % %end
% %
% %
% % % ------------------------------------------------------------------------
% % % convert blanks to NaN's, to avoid losing placeholders, and
% % % convert columns with numeric information to numeric vectors
% % % ------------------------------------------------------------------------
% % for i = 1:numc
% %
% %     % The replacing blanks part
% %     eval(['wh = strmatch('' '',str2mat(' names{i} '{:}));'])
% %     eval([names{i} '(wh) = {''NaN''};'])
% %
% %     % the conversion to numbers part
% %     eval(['a = str2num(str2mat(' names{i} '{:}));'])
% %     if isempty(a) | sum(isnan(a)) == length(a)
% %         % leave it alone; it's text
% %     else
% %         eval([names{i} ' = a;']),
% %     end
% % end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function coeffs = ellipseCartesianCoefficients(elli)
%ELLIPSECARTESIANCOEFFICIENTS Cartesian coefficients of an ellipse.
%
%   COEFFS = ellipseCartesianCoefficients(ELLI)
%   Computes the cartesian coefficients of the ellipse ELLI, given by:
%     COEFFS = [A B C D E F]
%   such that the points on the ellipse follow:
%     A*X^2 + B*X*Y + C*Y^2 + D*X + E*Y + F = 0
%
%   Example
%     elli = [30 20 40 20 30];
%     coeffs = ellipseCartesianCoefficients(elli)
%     elli2 = createEllipse(coeffs)
%     elli2 =
%        30.0000   20.0000   40.0000   20.0000   30.0000
%
%   See also
%     ellipses2d, createEllipse, equivalentEllipse
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2022-09-05, using Matlab 9.12.0.1884302 (R2022a)
% Copyright 2022 INRAE - BIA Research Unit - BIBS Platform (Nantes)

% retrieve ellipse center and squared radiusses
xc = elli(1);
yc = elli(2);
a2 = elli(3)^2;
b2 = elli(4)^2;

% pre-compute trigonometric functions (angle is in degrees)
cot = cos(elli(5) * pi / 180);
sit = sin(elli(5) * pi / 180);

% identification of each parameter
A = a2 * sit * sit + b2 * cot * cot;
B = 2 * (b2 - a2) * sit * cot;
C = a2 * cot * cot + b2 * sit * sit;
D = - 2 * A * xc - B * yc;
E = - B * xc - 2 * C * yc;
F = A * xc * xc + B * xc * yc + C * yc * yc - a2 * b2;

% concatenate into a single row vector
coeffs = [A B C D E F];

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [cfg, def] = cfg_mlbatch_appcfg(varargin)

% Add BasicIO to applications list of cfg_util. This file is an example how
% to add your own application configuration to cfg_util. To add an
% application, create a file called cfg_mlbatch_appcfg.m in the application
% folder and add this folder to the MATLAB path. cfg_util will look for
% files with the exact name cfg_mlbatch_appcfg.m and run all of them in
% order of their occurence on the path.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_mlbatch_appcfg.m 3179 2009-06-03 12:41:21Z volkmar $

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

if ~isdeployed
    % get path to this file
    p = fileparts(mfilename('fullpath'));
    % in this case, the cfg_mlbatch_appcfg file is not in the application
    % folder. Therefore, the path to the application needs to be added
    % explicitly
    addpath(fullfile(p, 'cfg_basicio'));
end
% these two files are now on MATLABs path
cfg = cfg_cfg_basicio;
def = cfg_cfg_basicio_def;

author: spm


text:
function value = p17_i4 ( action, name, value )

%*****************************************************************************80
%
%% P17_I4 sets or gets I4 parameters for problem 17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION,
%    'D' to set a parameter to its default value.
%    'S' to set a parameter,
%    'G' to get a parameter,
%    'I' to increment a parameter.
%
%    Input, string NAME, the name of the variable.
%    '#' is the number of calls to the integrand routine.
%
%    Input/output, integer VALUE.
%    * If ACTION = 'I', then VALUE is an input quantity, and is the
%      new value to be added to NAME.
%    * If ACTION = 'S', then VALUE is an input quantity, and is the
%      new value to be assigned to NAME.
%    * If ACTION = 'G', then VALUE is an output quantity, and is the
%      current value of NAME.
%
  persistent calls;

  if ( action == 'D' | action == 'd' )

    if ( name == '#' | name == '*' )
      calls = 0;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == '#' )
      value = calls;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P17_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P17_I4 - Fatal error!' );
    end

  elseif ( action == 'I' | action == 'i' )

    if ( name == '#' )
      calls = calls + value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P17_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P17_I4 - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == '#' )
      calls = value;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P17_I4 - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P17_I4 - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P17_I4 - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P17_I4 - Fatal error!' );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [hd_record]=mahdecomp(beta_gibbs,delta_gibbs,D_record,strshocks_record,It,Bu,Y,X,Z,n,m,p,k1,k3,T)

% function [hd_record]=mahdecomp(beta_gibbs,delta_gibbs,sigma_gibbs,D_record,It,Bu,Y,X,Z,n,m,p,k1,k3,T)
% performs algortihm 3.2.1, and returns posterior draws from the historical decomposition of the data sample
% inputs:  - matrix 'beta_gibbs': the matrix recording the post-burn draws of beta
%          - matrix 'delta_gibbs': the matrix recording the post-burn draws of delta
%          - matrix 'D_record': the matrix recording the simulated values of the D matrix
%          - integer 'It': the total number of iterations run by the Gibbs sampler
%          - integer 'Bu': the number of initial iterations discared as burn-in sample
%          - matrix 'Y': the matrix of endogenous variables, defined in (3.5.10)
%          - matrix 'X': the matrix of endogenous regressors, defined in (3.5.10)
%          - matrix 'Z': the matrix of exogenous regressors, defined in (3.5.10)
%          - integer 'n': the number of endogenous variables in the model
%          - integer 'm': the number of exogenous variables in the model
%          - integer 'p': the number of lags in the model
%          - integer 'k1': the number of coefficients related to the endogenous variables for each equation in the model
%          - integer 'k3': the number of coefficients related to the exogenous variables for each equation, in the reformulated model (3.5.5)
%          - integer 'T': the sample size, i.e. the number of time periods used to estimate the model
% outputs: - cell 'hd_record': the cell array containing records of simulated  hisotrical decompositions

% this function implements algorithm 3.2.1, adapted to the mean-adjusted VAR model

% preliminary tasks
% first create the hd_record and temp cells
hd_record=cell(n,n+1);
temp=cell(n,2);

% then initiate the Gibbs algorithm
for ii=1:It-Bu

% step 2: recover parameters
beta=beta_gibbs(:,ii);
B=reshape(beta,k1,n);
delta=delta_gibbs(:,ii);
Delta=reshape(delta,k3,n);
D=reshape(D_record(:,ii),n,n);

% step 3: obtain irfs and orthogonalised irfs
[~,ortirfmatrix]=bear.mairfsim(B,D,p,n,T);

% step 5: compute the historical contribution of each shock
   % fill the Yhd matrices
   % loop over rows of temp
   for jj=1:n
      % loop over columns of temp
      for kk=1:n
      % create the virf and vshocks vectors (shocks correspond to step 4)
         for ll=1:T
         virf(ll,1)=ortirfmatrix(jj,kk,ll);
         end
      vshocks=strshocks_record{kk,1}(ii,:)';
         % loop over sample periods
         for ll=1:T
         hd_record{jj,kk}(ii,ll)=virf(1:ll,1)'*flipud(vshocks(1:ll,1));
         end
      end
   end

% then go for next Gibbs iteration
end

% step 6: compute the contributions of deterministic variables
% loop over rows of temp/hd_record
for ii=1:n
% fill the Ytot matrix in temp
% initial condition
temp{ii,1}=hd_record{ii,1};
   % sum over the remaining columns of hd_record
   for jj=2:n
   temp{ii,1}=temp{ii,1}+hd_record{ii,jj};
   end
% fill the Y matrix in temp
temp{ii,2}=repmat(Y(:,ii)',It-Bu,1);
% fill the Yd matrix in hd_record
hd_record{ii,n+1}=temp{ii,2}-temp{ii,1};
% go for next variable
end

author: european-central-bank

repo: BEAR-toolbox


text:
function model = L1APG_initialize(img, id, x1, y1, x2, y2)

init_pos = [y1, y2, y1; x1 x1 x2];
sz_T = [12 15];
% sz_T = [24 12];

% parameters setting for tracking
para.lambda = [0.2, 0.001, 10]; % lambda 1, lambda 2 for a_T and a_I respectively, lambda 3 for the L2 norm parameter
% set para.lambda = [a,a,0]; then this the old model
para.angle_threshold = 40;
para.Lip = 8;
para.Maxit = 5;
para.nT = 10; % number of templates for the sparse representation
% para.rel_std_afnv = [0.03,0.0005,0.0005,0.03,1,1]; % diviation of the sampling of particle filter
para.rel_std_afnv = [0.003, 0.0005, 0.0005, 0.003, 1, 1];
para.n_sample = 300;		% number of particles
para.sz_T = sz_T;
para.init_pos = init_pos;
para.bDebug = 1;		% debugging indicator

% generate the initial templates for the 1st frame
[T, T_norm, T_mean, T_std] = InitTemplates(sz_T, para.nT, img, init_pos);
norms = T_norm .* T_std; % template norms

% get affine transformation parameters from the corner points in the first frame
aff_obj = corners2affine(init_pos, sz_T);
map_aff = aff_obj.afnv;

dim_T = size(T,1);	% number of elements in one template, sz_T(1)*sz_T(2)=12x15 = 180
A = [T eye(dim_T)]; % data matrix is composed of T, positive trivial T.
% fixT = T(:,1)/para.nT; % first template is used as a fixed template
%Temaplate Matrix
Temp = A;
Dict = Temp'*Temp;
Temp1 = T*pinv(T);

% build model
model.T = T;
model.T_mean = T_mean;
model.norms = norms;
model.occlusionNf = 0;
model.map_aff = map_aff;
model.A = A;
model.Temp = Temp;
model.Dict = Dict;
model.Temp1 = Temp1;
model.Lambda = para.lambda;
model.para = para;
model.id = id;
model.lost = 0;
model.aff_samples = ones(para.n_sample,1)*map_aff;

author: yuxng


text:
function model = assignQuantDir(model)
% Quantitatively assigns reaction directionality based on estimated bounds
% on transformed reaction Gibbs energies
%
% USAGE:
%
%    model = assignQuantDir(model)
%
% INPUT:
%    model:    structure with fields:
%
%                * .SIntRxnBool - `n x 1` boolean of internal reactions
%                * .DrGtMin - `n x 1` array of estimated lower bounds on
%                  transformed reaction Gibbs energies.
%                * .DrGtMax - `n x 1` array of estimated upper bounds on
%                  transformed reaction Gibbs energies.
%
% OUTPUT:
%    model:    structure with fields:
%
%                * .quantDir - `n x 1` array indicating quantitatively assigned
%                  reaction directionality. 1 for reactions that are
%                  irreversible in the forward direction, -1 for
%                  reactions that are irreversible in the reverse
%                  direction, and 0 for reversible reactions.
%
% .. Author: - Hulda SH, Nov. 2012

model.NaNdfG0MetBool=isnan(model.DfGt0);
model.NaNd0GRxnBool = isnan(model.DrGt0);

DrGtMin=model.DrGtMin;
DrGtMax=model.DrGtMax;

[nMet,nRxn]=size(model.S);

% model.quantDir = zeros(size(DrGtMin));
% model.quantDir(DrGtMax < 0) = 1;
% model.quantDir(DrGtMin > 0) = -1;

nEqualDrGt=nnz(DrGtMin==DrGtMax & DrGtMin~=0);
if any(nEqualDrGt)
    fprintf('%s\n',[num2str(nEqualDrGt) '/' num2str(length(DrGtMin)) ' reactions with DrGtMin=DrGtMax~=0' ]);
end

nZeroDrGt=nnz(DrGtMin==0 & DrGtMax==0);
if any(nZeroDrGt)
    fprintf('%s\n',[num2str(nZeroDrGt) '/' num2str(length(DrGtMin)) ' reactions with DrGtMin=DrGtMax=0' ]);
end

if any(model.NaNd0GRxnBool)
    warning('Some DfGt0 are NaN');
end

if any(DrGtMin>DrGtMax)
    error('DrGtMin greater than DrGtMax');
end

%reaction directionality
%keep exchange bounds the same as for the recostruction
model.lb_reconThermo=model.lb;
model.ub_reconThermo=model.ub;

%now set internal reaction directions
for n=1:nRxn
    if model.SIntRxnBool(n)
        if model.NaNd0GRxnBool(n)
            %for the reactions that involve a NaN metabolite standard Gibbs energy of
            %formation, use the directions given by the reconstruction
            if model.lb(n)<0 && model.ub(n)>0
                model.lb_reconThermo(n)=-Inf;
                model.ub_reconThermo(n)=Inf;
            end
            %forward
            if model.lb(n)>=0 && model.ub(n)>0
                model.lb_reconThermo(n)=0;
                model.ub_reconThermo(n)=Inf;
            end
            %reverse
            if model.lb(n)<0 && model.ub(n)<=0
                model.lb_reconThermo(n)=-Inf;
                model.ub_reconThermo(n)=0;
            end
            if model.lb(n)==0 && model.ub(n)==0
                error(['Reaction ' model.rxns{n} ' bounds set to zero'])
            end
            %note that there is no thermodynamic directionality assignment
            %for this reaction
            model.directionalityThermo{n}=NaN;
        else
            if DrGtMax(n)<0
                model.directionalityThermo{n}='forward';
                model.lb_reconThermo(n)=0;
                model.ub_reconThermo(n)=Inf;
            end
            if DrGtMin(n)>0
                model.directionalityThermo{n}='reverse';
                model.lb_reconThermo(n)=-Inf;
                model.ub_reconThermo(n)=0;
            end
            if DrGtMin(n)<0 && DrGtMax(n)>0
                model.directionalityThermo{n}='reversible';
                model.lb_reconThermo(n)=-Inf;
                model.ub_reconThermo(n)=Inf;
            end
            if DrGtMin(n)==DrGtMax(n)
                model.directionalityThermo{n}='equilibrium';
            end
            if model.lb(n)==0 && model.ub(n)==0
                error(['Reaction ' model.rxns{n} ' bounds set to zero'])
            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function [r, pol, res, zer, zj, fj, wj, errvec, wt] = aaatrig(F, varargin)
%AAATRIG   Trigonometric AAA and AAA-Lawson (near-minimax) real or complex
%      rational approximation.
%   R = AAATRIG(F, Z) computes an trigonometric AAA rational approximant R
%   (function handle) to data F on the set of sample points Z.  The rational
%   approximant is periodic with period 2*pi. F may be given by its values at
%   Z, or as a function handle or a chebfun.  R = AAATRIG(F, Z, 'degree', N)
%   computes the minimax approximation of degree N (i.e., rational type
%   (N,N)).
%
%   [R, POL, RES, ZER] = AAATRIG(F, Z) returns vectors of poles POL, residues
%   RES, and zeros ZER of R. The poles and zeros are repeated at intervals
%   of 2*pi.
%
%   [R, POL, RES, ZER, ZJ, FJ, WJ] = AAATRIG(F, Z) also returns the vectors of
%   support points ZJ, approximation values FJ = r(ZJ), and weights WJ of
%   the barycentric representation of R.
%
%   [R, POL, RES, ZER, ZJ, FJ, WJ, ERRVEC] = AAATRIG(F, Z) also returns the
%   vector of errors ||f-r||_infty in successive iteration steps of AAATRIG.
%
%   R = AAATRIG(F,Z,FORM) computes a rational approximant of type FORM.
%   FORM can either be 'odd' (default) or 'even'.
%
%   R = AAATRIG(F, Z, NAME, VALUE) sets the following parameters:
%   - 'tol', TOL: relative tolerance (default TOL = 1e-13),
%   - 'degree', N: maximal degree (default N = 99).
%      Output rational approximant will be at most of type (N,N).
%      Identical to 'mmax', N+1.
%      By default, this will turn on Lawson iteration: see next paragraph.
%   - 'mmax', MMAX: maximal number of terms in the barycentric representation
%       (default MMAX = 100). R will be of degree MMAX-1.
%       Identical to 'degree', MMAX-1. Also turns on Lawson iteration.
%   - 'dom', DOM: domain (default DOM = [0, 2*pi]). No effect if Z is provided.
%   - 'cleanup', 'off' or 0: turns off automatic removal of numerical Froissart
%       doublets
%   - 'cleanuptol', CLEANUPTOL: cleanup tolerance (default CLEANUPTOL = TOL).
%       Poles with residues less than this number times the geometric mean size
%       of F times the minimum distance to Z are deemed spurious by the cleanup
%       procedure. If TOL = 0, then CLEANUPTOL defaults to 1e-13.
%   - 'lawson', NLAWSON: take NLAWSON iteratively reweighted least-squares steps
%       to bring approximation closer to minimax; specifying NLAWSON = 0
%       ensures there is no Lawson iteration.  See next paragraph.
%
%   If 'degree' or equivalently 'mmax' is specified and 'lawson' is not, then
%   AAATRIG attempts to find a minimax approximant of degree N by Lawson iteration.
%   This will generally be successful only if the minimax error is well
%   above machine precision, and is more reliable for complex problems than
%   real ones.  If 'degree' and 'lawson' are both specified, then exactly
%   NLAWSON Lawson steps are taken (so NLAWSON = 0 corresponds to AAA
%   approximation with no Lawson iteration).  The final weight vector WT of
%   the Lawson iteration is available with
%   [R, POL, RES, ZER, ZJ, FJ, WJ, ERRVEC, WT] = AAATRIG(F, Z).
%
%   Note that R may have fewer than N poles and zeros.  This may happen,
%   for example, if N is too large, or if F is even and N is odd, or if F is
%   odd and N is even.
%
%   One can also execute R = AAATRIG(F), with no specification of a set Z.
%   If F is a vector, this is equivalent to R = AAATRIG(F, Z) with
%   Z = LINSPACE(0, 2*pi, LENGTH(F)).  If F is a function handle or a chebfun,
%   AAATRIG attempts to resolve F on its domain, which defaults to [0,2*pi] for
%   a function handle.
%
% Examples:
%
%   f = @(x) exp(cos(x));
%   r = aaatrig(f); xx = linspace(-pi,pi); plot(xx,r(xx)-f(xx))
%
%   r = aaatrig(f,'degree',4); xx = linspace(-pi,pi); plot(xx,r(xx)-f(xx))
%
%   Z = exp(2i*pi*linspace(0,1,500));
%   [r,pol,res] = aaatrig(@tan,Z); disp([pol res])
%
%   X = linspace(0,2*pi,1000); F = gamma(2+sin(X));
%   subplot(1,2,1)
%   r = aaatrig(F,X,'degree',10,'lawson',0); plot(X,F-r(X)), hold on
%   r = aaatrig(F,X,'degree',10); plot(X,F-r(X)), hold off
%
%   Z = exp(1i*linspace(0,2*pi,1000)); G = exp(1i*tan(Z));
%   subplot(1,2,2)
%   r = aaatrig(G,Z,'degree',5,'lawson',0); plot(G-r(Z)), axis equal, hold on
%   r = aaatrig(G,Z,'degree',5); plot(G-r(Z)), axis equal, hold off
%
%   References:
%   [1] Yuji Nakatsukasa, Olivier Sete, Lloyd N. Trefethen, "The AAA algorithm
%   for rational approximation", SIAM J. Sci. Comp. 40 (2018), A1494-A1522.
%
%   [2] Yuji Nakatsukasa and Lloyd N. Trefethen, "An algorithm for real and
%   complex rational minimax approximation", SIAM J. Sci. Comp. (2020).
%
%   [3] Peter J. Baddoo, "The AAAtrig algorithm for rational approximation
%   of periodic functions", SIAM J. Sci. Comp. (2021).
%
% See also AAA, TRIGRATINTERP, CHEBPADE, MINIMAX, PADEAPPROX.

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

% Parse inputs:
[F, Z, M, form, dom, tol, mmax, cleanup_flag, cleanup_tol, needZ, mmax_flag, nlawson] ...
    = parseInputs(F, varargin{:});

if ( needZ )
    % Z was not provided.  Try to resolve F on its domain.
    [r, pol, res, zer, zj, fj, wj, errvec] = ...
        aaatrig_autoZ(F, form, dom, tol, mmax, cleanup_flag, cleanup_tol, mmax_flag, nlawson);
    return
end

% Remove any infinite or NaN function values (avoid SVD failures):
toKeep = ~isinf(F);
F = F(toKeep); Z = Z(toKeep);
toKeep = ~isnan(F);
F = F(toKeep); Z = Z(toKeep);

% Find sample points at infity
infP = find(Z==+1i*Inf);
infM = find(Z==-1i*Inf);
finfP = F(infP); % value at +1i*Inf
finfM = F(infM); % value at -1i*Inf
F([infP,infM])=[]; Z([infP,infM])=[];
if strcmp(form,'even') && numel([finfP,finfM])==2 && finfP~=finfM
    error('The even representation must take the same values at +/-i*Inf.')
end

% Define basis functions
if strcmp(form,'even')
    cst = @(x) cot(x);
elseif strcmp(form,'odd')
    cst = @(x) csc(x);
end

% Project sample points onto a single period window
Z = Z - 2*pi*floor(real(Z/(2*pi)));

% Remove repeated elements of Z and corresponding elements of F:
[Z, uni] = unique(Z,'stable'); F = F(uni);

M = length(Z);

% Relative tolerance:
reltol = tol * norm(F, inf);

% Left scaling matrix:
SF = spdiags(F, 0, M, M);

% Initialization for AAA iteration:
J = 1:M;
zj = [];
fj = [];
C = [];
errvec = [];
R = mean(F);

% AAA iteration:
for m = 1:mmax
    % Select next support point where error is largest:
    [~, jj] = max(abs(F - R));          % Select next support point.
    zj = [zj; Z(jj)];                   % Update support points.
    fj = [fj; F(jj)];                   % Update data values.
    J(J == jj) = [];                    % Update index vector.
    C = [C cst((Z - Z(jj))/2)];
    % Compute weights:
    Sf = diag(fj);                      % Right scaling matrix.
    A = SF*C - C*Sf;                    % Loewner matrix.
    % Add value(s) at infinity to least-squares problem
    Jv = J;
    if ~isempty([finfP,finfM])
        if strcmp(form,'odd')
            if ~isempty(finfP)
                A = [A; (finfP - fj.').*exp(-1i*zj.'/2)];
                Jv = [Jv, size(A,1)];
            end
            if ~isempty(finfM)
                A = [A; (finfM - fj.').*exp(+1i*zj.'/2)];
                Jv = [Jv, size(A,1)];
            end
        elseif strcmp(form,'even')
        A = [A; (finfP-fj.')];
        Jv = [Jv, size(A,1)];
        end
    end
    [~, ~, V] = svd(A(Jv,:), 0); % Reduced SVD. Includes the bottom
    wj = V(:,m);                        % weight vector = min sing vector
    % Rational approximant on Z:
    N = C*(wj.*fj);                     % Numerator
    D = C*wj;                           % Denominator
    R = F;
    R(J) = N(J)./D(J);

    % Error in the sample points and at infinity:
    err = F - R;
    if isempty([finfP,finfM])
    maxerr = norm(1, inf);
    else
        if strcmp(form,'odd')
            if ~isempty(finfP)
            errInfP = finfP - (fj.'.*exp(-1i*zj.'/2)*wj)./(exp(-1i*zj.'/2)*wj); % error at +1i*infinity
            err = [err; errInfP];
            end
            if ~isempty(finfM)
            errInfM = finfM - (fj.'.*exp(+1i*zj.'/2)*wj)./(exp(+1i*zj.'/2)*wj); % error at -1i*infinity
            err = [err; errInfM];
            end
        elseif strcmp(form,'even')
            errInf = finfP - fj.'*wj./sum(wj); % error at +/-1i*infinity
            err = [err; errInf];
        end
    end
    maxerr = norm(err,inf);
    errvec = [errvec; maxerr];
    % Check if converged:
    if ( maxerr <= reltol )
        break
    end
end

maxerrAAA = maxerr;                     % error at end of AAA

% When M == 2, one weight is zero and r is constant.
% To obtain a good approximation, interpolate in both sample points.
if ( M == 2 )
    zj = Z;
    fj = F;
    wj = [1; -1];
    wj = wj/norm(wj);   % Impose norm(w) = 1 for consistency.
    errvec(2) = 0;
    maxerrAAA = 0;
end

% We now enter Lawson iteration: barycentric IRLS = iteratively reweighted
% least-squares if 'lawson' is specified with NLAWSON > 0 or 'mmax' is
% specified and 'lawson' is not.  In the latter case the number of steps
% is chosen adaptively.  Note that the Lawson iteration is unlikely to be
% successful when the errors are close to machine precision.

wj0 = wj; fj0 = fj;     % Save parameters in case Lawson fails
wt = NaN(M,1); wt_new = ones(M,1);
if ( nlawson > 0 )      % Lawson iteration

  if ~isempty([finfP,finfM])
    warning('Specifying the function values at infinity is not currently compatible with Lawson iteration.')
  end

    maxerrold = maxerrAAA;
    maxerr = maxerrold;
    nj = length(zj);
    A = [];
    for j = 1:nj                              % Cauchy/Loewner matrix
        A = [A cst((Z-zj(j))/2) F.*cst((Z-zj(j))/2)];
    end
    for j = 1:nj
        [i,~] = find(Z==zj(j));               % support pt rows are special
        A(i,:) = 0;
        A(i,2*j-1) = 2;
        A(i,2*j) = 2*F(i);
    end
    stepno = 0;
    while ( (nlawson < inf) & (stepno < nlawson) ) |...
          ( (nlawson == inf) & (stepno < 20) ) |...
          ( (nlawson == inf) & (maxerr/maxerrold < .999) & (stepno < 1000) )
        stepno = stepno + 1;
        wt = wt_new;
        W = spdiags(sqrt(wt),0,M,M);
        [U,S,V] = svd(W*A,0);
        c = V(:,end);
        denom = zeros(M,1); num = zeros(M,1);
        for j = 1:nj
            denom = denom + c(2*j).*cst((Z-zj(j))/2);
            num = num - c(2*j-1).*cst((Z-zj(j))/2);
        end
        R = num./denom;
        for j = 1:nj
            [i,~] = find(Z==zj(j));           % support pt rows are special
            R(i) = -c(2*j-1)/c(2*j);
        end
        err = F - R; abserr = abs(err);
        wt_new = wt.*abserr; wt_new = wt_new/norm(wt_new,inf);
        maxerrold = maxerr;
        maxerr = max(abserr);
    end
    wj = c(2:2:end);
    fj = -c(1:2:end)./wj;
    % If Lawson has not reduced the error, return to pre-Lawson values.
    if (maxerr > maxerrAAA) & (nlawson == Inf)
        wj = wj0; fj = fj0;
    end
end
% Remove support points with zero weight:
I = find(wj == 0);
zj(I) = [];
wj(I) = [];
fj(I) = [];

% Construct function handle:
r = @(zz) revaltrig(zz, zj, fj, wj, form);

% Compute poles, residues and zeros:
[pol, res, zer] = prztrig(zj, fj, wj, form);

if ( cleanup_flag & nlawson == 0)       % Remove Froissart doublets
    [r, pol, res, zer, zj, fj, wj] = ...
        cleanuptrig(r, form, finfP, finfM,pol, res, zer, zj, fj, wj, Z, F, cleanup_tol);
end

end % of AAATRIG()

%% parse Inputs:

function [F, Z, M, form, dom, tol, mmax, cleanup_flag, cleanup_tol, ...
    needZ, mmax_flag, nlawson] = parseInputs(F, varargin)
% Input parsing for AAATRIG.

% Check if F is empty:
if ( isempty(F) )
    error('CHEBFUN:aaatrig:emptyF', 'No function given.')
elseif ( isa(F, 'chebfun') )
    if ( size(F, 2) ~= 1 )
        error('CHEBFUN:aaatrig:nColF', 'Input chebfun must have one column.')
    end
end

% Sample points:
if ( ~isempty(varargin) && isfloat(varargin{1}) )
    % Z is given.
    Z = varargin{1};
    if ( isempty(Z) )
        error('CHEBFUN:aaatrig:emptyZ', ...
            'If sample set is provided, it must be nonempty.')
    end
    varargin(1) = [];
end

% Set defaults for other parameters:
tol = 1e-13;         % Relative tolerance.
mmax = 100;          % Maximum number of terms.
cleanup_tol = 1e-13; % Cleanup tolerance.
nlawson = Inf;       % number of Lawson steps (Inf means adaptive)
% Domain:
if ( isa(F, 'chebfun') )
    dom = F.domain([1, end]);
else
    dom = [0, 2*pi];
end
cleanup_flag = 1;   % Cleanup on.
mmax_flag = 0;      % Checks if mmax manually specified.
cleanup_set = 0;    % Checks if cleanup_tol manually specified.
form = 'odd';
% Check if parameters have been provided:
while ( ~isempty(varargin) )

    if  strncmpi(varargin{1},'even',4)
          form = 'even';
          varargin(1) = [];
    elseif strncmpi(varargin{1},'odd',3)
          form = 'odd';
          varargin(1) = [];

    elseif ( strncmpi(varargin{1}, 'tol', 3) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 1]) )
            tol = varargin{2};
            if ~cleanup_set & tol > 0 % If not manually set, set cleanup_tol to tol.
              cleanup_tol = tol;
            end
        end
        varargin([1, 2]) = [];

    elseif ( strncmpi(varargin{1}, 'degree', 6) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 1]) )
            if ( mmax_flag == 1 ) && ( mmax ~= varargin{2}+1 )
                error('CHEBFUN:aaatrig:degmmaxmismatch', ' mmax must equal degree+1.')
            end
            mmax = varargin{2}+1;
            mmax_flag = 1;
        end
        varargin([1, 2]) = [];

    elseif ( strncmpi(varargin{1}, 'mmax', 4) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 1]) )
            if ( mmax_flag == 1 ) && ( mmax ~= varargin{2})
                error('CHEBFUN:aaatrig:degmmaxmismatch', ' mmax must equal degree+1.')
            end
            mmax = varargin{2};
            mmax_flag = 1;
        end
        varargin([1, 2]) = [];

    elseif ( strncmpi(varargin{1}, 'lawson', 6) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 1]) )
            nlawson = varargin{2};
        end
        varargin([1, 2]) = [];

    elseif ( strncmpi(varargin{1}, 'dom', 3) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 2]) )
            dom = varargin{2};
        end
        varargin([1, 2]) = [];
        if ( isa(F, 'chebfun') )
            if ( ~isequal(dom, F.domain([1, end])) )
                warning('CHEBFUN:aaatrig:dom', ...
                    ['Given domain does not match the domain of the chebfun.\n', ...
                    'Results may be inaccurate.'])
            end
        end

    elseif ( strncmpi(varargin{1}, 'cleanuptol', 10) )
        if ( isfloat(varargin{2}) && isequal(size(varargin{2}), [1, 1]) )
          cleanup_tol = varargin{2};
          cleanup_set = 1;
        end
        varargin([1, 2]) = [];

    elseif ( strncmpi(varargin{1}, 'cleanup', 7) )
        if ( strncmpi(varargin{2}, 'off', 3) || ( varargin{2} == 0 ) )
            cleanup_flag = 0;
        end
        varargin([1, 2]) = [];

    else
        error('CHEBFUN:aaatrig:UnknownArg', 'Argument unknown.')
    end
end

% Deal with Z and F:
if ( ~exist('Z', 'var') && isfloat(F) )
    % F is given as data values, pick same number of sample points:
    Z = linspace(dom(1), dom(2), length(F)).';
end

if ( exist('Z', 'var') )
    % Z is given:
    needZ = 0;

    % Work with column vector:
    Z = Z(:);
    M = length(Z);

    % Function values:
    if ( isa(F, 'function_handle') || isa(F, 'chebfun') )
        % Sample F on Z:
        F = F(Z);
    elseif ( isnumeric(F) )
        % Work with column vector and check that it has correct length.
        F = F(:);
        if ( length(F) ~= M )
            error('CHEBFUN:aaatrig:lengthFZ', ...
                'Inputs F and Z must have the same length.')
        end
    elseif ( ischar(F) )
        % F is given as a string input. Convert it to a function handle.
        F = str2op(vectorize(F));
        F = F(Z);
    else
        error('CHEBFUN:aaatrig:UnknownF', 'Input for F not recognized.')
    end

else
    % Z was not given.  Set flag that Z needs to be determined.
    % Also set Z and M since they are needed as output.
    needZ = 1;
    Z = [];
    M = length(Z);
end

if ~mmax_flag & (nlawson == Inf)
    nlawson = 0;
end

end % End of PARSEINPUT().

%% Cleanup.  In June 2022 the residue size test was changed to be relative to
%            the distance to the approximation set Z.

function [r, pol, res, zer, z, f, w] = ...
    cleanuptrig(r, form, finfP, finfM, pol, res, zer, z, f, w, Z, F, cleanup_tol)
% Remove spurious pole-zero pairs.

% Find negligible residues:
if any(F)
   geometric_mean_of_absF = exp(mean(log(abs(F(F~=0)))));
else
   geometric_mean_of_absF = 0;
end
Zdistances = NaN(size(pol));
for j = 1:length(Zdistances);
   Zdistances(j) = min(abs(pol(j)-Z));
end
ii = find(abs(res)./Zdistances < cleanup_tol * geometric_mean_of_absF);
ni = length(ii);
if ( ni == 0 )
    % Nothing to do.
    return
elseif ( ni == 1 )
    warning('CHEBFUN:aaatrig:Froissart','1 Froissart doublet');
else
    warning('CHEBFUN:aaatrig:Froissart',[int2str(ni) ' Froissart doublets']);
end
% For each spurious pole find and remove closest support point.
for j = 1:ni
    % Find the closest support point modulo the period
    np = fix(real((z-pol(ii(j)))/pi));
    azp = abs(z-(pol(ii(j)) + np*2*pi));
    jj = find(azp == min(azp),1);
    % Remove support point(s):
    z(jj) = [];
    f(jj) = [];
end

% Remove support points z from sample set:
for jj = 1:length(z)
    F(Z == z(jj)) = [];
    Z(Z == z(jj)) = [];
end

m = length(z);
M = length(Z);

% Define basis functions
if strcmp(form,'even')
    cst = @(x) cot(x);
elseif strcmp(form,'odd')
    cst = @(x) csc(x);
end

% Build Loewner matrix:
SF = spdiags(F, 0, M, M);
Sf = diag(f);
C = cst(bsxfun(@minus, Z, z.')/2);      % Cauchy matrix.
A = SF*C - C*Sf;                    % Loewner matrix.

% Enforce value(s) at infinity
if ~isempty([finfP,finfM])
    if strcmp(form,'odd')
        if ~isempty(finfP)
            A = [A; (finfP - f.').*exp(-1i*z.'/2)];
        end
        if ~isempty(finfM)
            A = [A; (finfM - f.').*exp(+1i*z.'/2)];
        end
    elseif strcmp(form,'even')
    A = [A; (finfP-f.')];   % Add a row to A to enforce behaviour at infinity.
    end
end

% Solve least-squares problem to obtain weights:
[~, ~, V] = svd(A, 0);
w = V(:,m);

% Build function handle and compute poles, residues and zeros:
r = @(zz) revaltrig(zz, z, f, w, form);
[pol, res, zer] = prztrig(z, f, w, form);

end % End of CLEANUPTRIG().

%% Automated choice of sample set

function [r, pol, res, zer, zj, fj, wj, errvec] = ...
    aaatrig_autoZ(F, form, dom, tol, mmax, cleanup_flag, cleanup_tol, mmax_flag, nlawson)
%

% Flag if function has been resolved:
isResolved = 0;

% Main loop:
for n = 5:14
    % Sample points:
    % Next line enables us to do pretty well near poles
    Z = linspace(dom(1)+1.37e-7*diff(dom), dom(2), 2 + 2^n).'; Z(end) = [];
    [r, pol, res, zer, zj, fj, wj, errvec] = aaatrig(F, Z, form, 'tol', tol, ...
        'mmax', mmax, 'cleanup', cleanup_flag, 'cleanuptol', cleanup_tol, 'lawson', nlawson);

    % Test if rational approximant is accurate:
    reltol = tol * norm(F(Z), inf);

    % On Z(n):
    err(1,1) = norm(F(Z) - r(Z), inf);
    Zrefined = linspace(dom(1)+1.37e-8*diff(dom), dom(2)-3.08e-9*diff(dom), ...
        1 + round(1.5 * (1 + 2^(n+1)))).';
    Zrefined(end) = [];
    err(2,1) = norm(F(Zrefined) - r(Zrefined), inf);

    if ( all(err < reltol) )
        % Final check that the function is resolved, inspired by sampleTest().
        % Pseudo random sample points in [-1, 1]:
        xeval = [-0.357998918959666; 0.036785641195074];
        % Scale to dom:
        xeval = (dom(2) - dom(1))/2 * xeval + (dom(2) + dom(1))/2;

        if ( norm(F(xeval) - r(xeval), inf) < reltol )
            isResolved = 1;
            break
        end
    end
end

if ( ( isResolved == 0 ) && ~mmax_flag )
    warning('CHEBFUN:aaatrig:notResolved', ...
        'Function not resolved using %d pts.', length(Z))
end

end % End of AAATRIG_AUTOZ().

function op = str2op(op)
    % Convert string inputs to either numeric format or function_handles.
    sop = str2num(op);
    if ( ~isempty(sop) )
        op = sop;
    else
        depVar = symvar(op);
        if ( numel(depVar) ~= 1 )
            error('CHEBFUN:CHEBFUN:str2op:indepvars', ...
             'Incorrect number of independent variables in string input.');
        end
        op = eval(['@(' depVar{:} ')', op]);
    end
end % End of STR2OP().

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [DAD,DInv]=equilibrateSymMat(A)
%EQUILIBRATESYMMAT Equilibrate a symmetric matrix maintaining symmetry.
%          This finds a matrix D such that D*A*D is such that the largest
%          magnitude element in each row and in each column is 1. It is
%          assumed that there is no all-zero row in A.
%
%INPUTS: A An nXn symmetric matrix having no all-zero row.
%
%OUTPUTS: DAD The equilibrated matrix. That is A pre and post multiplied by
%             the diagonal matrix D.
%        DInv The inverse of the diagonal matrix D. Invert the diagonal to
%             get D.
%
%This function implements the algorithm of Sectioj 5 of [1].
%
%EXAMPLE 1:
%Here, we create a symmetric matrix and equilibrate it. The equilibration
%makes the largest magnitude values in each row and column equal 1.
%However, it also worsens the condition number, as we shall show:
% V=vander([1;2;3;4]);
% A=V+V';
% DAD=equilibrateSymMat(A);
%
% %The original maximum absolute value elements in the columns and rows.
% max(abs(A),[],1)
% max(abs(A),[],2)
%
% %After equilibration, the maximum mangitude elements in the columns and
% %rows are all 1.
% max(abs(DAD),[],1)
% max(abs(DAD),[],2)
%
% %One can see that the condition number increased.
% cond(A,Inf)
% cond(DAD,Inf)
%
%EXAMPLE 2:
%This is the example given in Section 7 of [1]. AGian, we can see that the
%matrix has been equilibrated, but the condition number got worse.
% A=[1,2;
%    2,1];
% DAD=equilibrateSymMat(A);
%
% %The original maximum absolute value elements in the columns and rows.
% max(abs(A),[],1)
% max(abs(A),[],2)
%
% %After equilibration, the maximum mangitude elements in the columns and
% %rows are all 1.
% max(abs(DAD),[],1)
% max(abs(DAD),[],2)
%
% %One can see that the condition number increased.
% cond(A,Inf)%The inital condition number is 3.
% cond(DAD,Inf)%The equilibrated condition number is 16/3.
%
%REFERENCES:
%[1] J. R. Bunch, "Equilibration of symmetric matrices in the max-norm,"
%    Journal of the Association for Computing Machinery, vol. 18, no. 4,
%    pp. 566-572, Oct. 1971.
%
%September 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=size(A,1);

d=zeros(n,1);
for i=1:n
    d(i)=sqrt(abs(A(i,i)));
    for j=1:(i-1)
        t=abs(A(i,j));
        d(i)=max(t,d(i));
    end
    if(d(i)~=0)
        for j=1:i
            A(i,j)=A(i,j)/d(i);
        end

        for j=i:n
            A(j,i)=A(j,i)/d(i);
        end
    end
end

for i=1:n
    if(d(i)==0)
        for j=(i+1):n
            t=abs(A(j,i));
            d(i)=max(t,d(i));
        end

        if(d(i)==0)
            error('The matrix has an all-zero row.')
        end
        for j=(i+1):n
            A(j,i)=A(j,i)/d(i);
        end
    end
end

L=tril(A);
DAD=L+L'-diag(diag(L));

DInv=diag(d);

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:
%+++ Data load
load T2DM
%+++ Subwindow Permutation Analysis
F=spa(X,y,5,5,'autoscaling',300,0.8,10);
bar(F.COSS);
xlabel('Variable index');
figure;
plotspa(F,11);
vargrouping(F,[1:21],1)

author: viggin

repo: domain-adaptation-toolbox


text:
function [u_norm, coeffs, u] = smallestinconvexhull(M, x, U, tol)
% Computes a minimal norm convex combination of given tangent vectors in Manopt.
%
% function [u_norm, coeffs, u] = smallestinconvexhull(M, x, U)
% function [u_norm, coeffs, u] = smallestinconvexhull(M, x, U, tol)
%
% M is a manifold as returned by a Manopt factory.
% x is a point on this manifold.
% U is a cell containing N tangent vectors U{1} to U{N} at x.
% tol (default: 1e-8): tolerance for solving the quadratic program.
%
% This function computes u, a tangent vector at x contained in the convex
% hull spanned by the N vectors U{i}, with minimal norm (according to the
% Riemannian metric on M). This is obtained by solving a convex quadratic
% program involving the Gram matrix of the given tangent vectors.
% The quadratic program is solved using Matlab's built-in quadprog,
% which requires the optimization toolbox. If this toolbox is not
% available, consider replacing with CVX for example.
%
%
% u_norm is the norm of the smallest vector u.
% coeffs is a vector of length N with entries in [0, 1] summing to 1.
% u is the sought vector: u = coeffs(1)*U{1} + ... + coeffs(N)*U{N}.
%
% Nicolas Boumal, Feb. 19, 2013
% Modified April 6, 2016 to work with Manopt.

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, June 28, 2016.
% Contributors:
% Change log:
%
%   June 28, 2016 (NB):
%       Adapted for Manopt from original code by same author (Feb. 19, 2013)

% Example code: pick a manifold, a point, and a collection of tangent
% vectors at that point, then get the smallest vector in the convex hull
% of those:
%
% M = spherefactory(5);
% x = M.rand();
% N = 3;
% U = cell(N,1);
% for k = 1 : N, U{k} = M.randvec(x); end
% [u_norm, coeffs, u] = smallestinconvexhull(M, x, U)

    % We simply need to solve the following quadratic program:
    % minimize ||u||^2 such that u = sum_i s_i U_i, 0 <= s_i <= 1
    %                            and sum_i s_i = 1
    %
    % This is equivalent to solving:
    %  min s'*G*s s.t. 0 <= s <= 1, s'*ones = 1, with G(i, j) = <U_i, U_j> (Gram matrix)
    % Then our solution is s_1 U_1 + ... + s_N U_N.

    % Compute the Gram matrix of the given tangent vectors
    N = numel(U);
    G = grammatrix(M, x, U);

    % Solve the quadratic program.
    % If the optimization toolbox is not available, consider replacing with
    % CVX.

    if ~exist('tol', 'var') || isempty(tol)
        tol = 1e-8;
    end

    opts = optimset('Display', 'off', 'TolFun', tol);
    [s_opt, cost_opt] ...
          = quadprog(G, zeros(N, 1),     ...  % objective (squared norm)
                     [], [],             ...  % inequalities (none)
                     ones(1, N), 1,      ...  % equality (sum to 1)
                     zeros(N, 1),        ...  % lower bounds (s_i >= 0)
                     ones(N, 1),         ...  % upper bounds (s_i <= 1)
                     [],                 ...  % we do not specify an initial guess
                     opts);

    % Norm of the smallest tangent vector in the convex hull:
    u_norm = real(sqrt(2*cost_opt));

    % Keep track of optimal coefficients
    coeffs = s_opt;

    % If required, construct the vector explicitly.
    if nargout >= 3
        u = lincomb(M, x, U, coeffs);
    end

end

author: NicolasBoumal

repo: manopt


text:
% GCI_creak_postproc.m
% Function to do the post-processing step to remove false positive GCIs
% detected in creaky voice regions
%
% Octave compatible
%
% Description
% Function to do the post-processing step to remove false positive GCIs
% detected in creaky voice regions
%
% Inputs
%  GCI             : [samples] [Mx1] Glottal closure instants
%  creak           : [binary]  [Nx1] Creak decision
%  search_reg      : [integer] [1x1] Search factor
%  rep             : [samples] [Nx1] Resonator output
%  removeThresh    : [integer] [1x1] Threshold for false GCI removal
%  repNum          : [integer] [1x1] number of repititions
%
% Outputs
%  GCI             : [samples] [Mx1] Glottal closure instants
%
% Example
%  GCI = GCI_creak_postproc(GCI,creak,search_reg,rep,removeThresh,repNum)
%
% References
%  [1] Kane, J., Gobl, C., (2013) `Evaluation of glottal closure instant
%       detection in a range of voice qualities', Speech Communication
%       55(2), pp. 295-314.
%
% Copyright (c) 2013 Trinity College Dublin
%
% License
%  This code is a part of the Voice Analysis Toolkit with the following
%  licence:
%  The software product (and any modifications thereof) is distributed under
%  a dual licence: an open source license for individual, non-commercial
%  purposes, and a commercial license. The opensource licence under which
%  the product is distributed is GNU GPL v2. For individual users, this
%  licence suits their use as these are not distributing proprietary
%  modifications, additions to, or derivatives of the product and don't
%  require legal protection of a commercial licence. For commercial users,
%  where open source does not meet their requirements, we offer commercial
%  licensing of the product. A commercial license permits customers to
%  modify, add or produce derivative products without the obligation of
%  making the subsequent code open source. For more information regarding
%  our commercial licence, please contact john.whelan@tcd.ie
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  John Kane kanejo@tcd.ie
%
% $Id <info set by the versioning system> $

function GCI = GCI_creak_postproc(GCI,creak,search_reg,rep,removeThresh,repNum)

% Separate GCIs detected in creaky voice regions from other regions
creak_GCI=GCI(creak(GCI)==1);
GCI(creak(GCI)==1)=[];

for m=1:repNum
    n=2;
    while n < length(creak_GCI)
        cur_rep_max = abs(min(rep(creak_GCI(n)-round(search_reg):creak_GCI(n)+round(search_reg))));

        if mean([abs(rep(creak_GCI(n-1))) abs(rep(creak_GCI(n+1)))])*removeThresh > cur_rep_max
             creak_GCI(n)=NaN;
            n=n+2;
        else n=n+1;
        end
    end
    creak_GCI(isnan(creak_GCI))=[];
end

GCI=sort(unique([GCI(:)' creak_GCI(:)']));

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function varargout = drawPolygon3d(varargin)
%DRAWPOLYGON3D Draw a 3D polygon specified by a list of vertex coords.
%
%   drawPolygon3d(POLY);
%   packs coordinates in a single N-by-3 array.
%
%   drawPolygon3d(PX, PY, PZ);
%   specifies coordinates in separate numeric vectors (either row or
%   columns)
%
%   drawPolygon3d(..., PARAM, VALUE);
%   Specifies style options to draw the polyline, see plot for details.
%
%   H = drawPolygon3d(...);
%   also returns a handle to the list of created line objects.
%
%   Example
%     t = linspace(0, 2*pi, 100)';
%     xt = 10 * cos(t);
%     yt = 5 * sin(t);
%     zt = zeros(1,100);
%     figure; drawPolygon3d(xt, yt, zt, 'b');
%
%   See Also:
%   polygons3d, fillPolygon3d, drawPolyline3d
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2011-08-17 from drawPolyline3d, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

% HISTORY
% 2019-02-02 add support for multiple polygons

%% Process input arguments

% extract handle of axis to draw on
ax = gca;
var1 = varargin{1};
if isAxisHandle (var1)
    ax          = var1;
    varargin(1) = [];
end

% check case we want to draw several curves, stored in a cell array
var1 = varargin{1};
if iscell(var1)
    hold on;
    h = [];
    for i = 1:length(var1(:))
        h = [h; drawPolygon3d(ax, var1{i}, varargin{2:end})]; %#ok<AGROW>
    end
    if nargout > 0
        varargout{1} = h;
    end
    return;
end

%% extract polygon coordinates
if min(size(var1)) == 1
    % if first argument is a vector (either row or column), then assumes
    % first argument contains x coords, second argument contains y coords
    % and third one the z coords
    px = var1;
    if length(varargin) < 3
        error('geom3d:drawPolygon3d:Wrong number of arguments in fillPolygon3d');
    end
    py = varargin{2};
    pz = varargin{3};
    varargin = varargin(4:end);
else
    % first argument contains both coordinate
    px = var1(:, 1);
    py = var1(:, 2);
    pz = var1(:, 3);
    varargin = varargin(2:end);
end

%% draw the polygon

% check that the polygon is closed
if px(1) ~= px(end) || py(1) ~= py(end) || pz(1) ~= pz(end)
    px = [px(:); px(1)];
    py = [py(:); py(1)];
    pz = [pz(:); pz(1)];
end

% draw the closed curve
h = plot3(ax, px, py, pz, varargin{:});

%% Format output

if nargout > 0
    varargout = {h};
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
path = '../traindata/original_segment/';
for i = 1:19
    save_folder = ['../testdata/original_mask/' sprintf('%02d',i) '/'];
    mkdir(save_folder);

    dirs = dir(['../data/original_segment/' sprintf('%02d',i),'/*.mat']);
    dirs = struct2cell(dirs);
    for j = 2976:size(dirs,2)
        delete(['../data/original_segment/' sprintf('%02d',i),'/' dirs{1,j}])
        %load(['../data/original_segment/' sprintf('%02d',i),'/' dirs{1,j}]);
        %save([save_folder dirs{1,j}],'library_mask','library_mask_pole');

    end

end

author: xjqicuhk

repo: SIMS

save_path: github-repos/MATLAB/xjqicuhk-SIMS


text:
function voxel3(T,varargin)
%VOXEL3 Visualize a third-order tensor with voxels.
%   voxel3(T) visualizes the third-order tensor T by plotting its elements
%   as voxels whose color and opacity are proportional to their value. The
%   figure contains two sliders for setting the parameters thresh and
%   degree (press 'h' to hide/show them). Let alpha = (T(i,j,k)-min(T(:))/
%   (max(T(:))-min(T(:))), then the opaqueness of each voxel, where 0 is
%   transparent and 1 is opaque, is computed as
%
%      alpha                            if alpha >= thresh
%      thresh^(1-degree)*alpha^degree   if alpha <  thresh
%
%   voxel3(T,varargin) passes the parameters varargin{:} to the plot.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)

% Check the dimensions.
if ndims(T) ~= 3
    error('voxel3:T','ndims(T) should be 3.');
end

% Check options.
idx = find(strcmpi('thresh',varargin),1,'last');
if idx < length(varargin)
    thresh = varargin{idx+1};
    varargin = [varargin(1:idx-1) varargin(idx+2:end)];
end
idx = find(strcmpi('degree',varargin),1,'last');
if idx < length(varargin)
    degree = varargin{idx+1};
    varargin = [varargin(1:idx-1) varargin(idx+2:end)];
end

% If the data is complex, convert it to the modulus.
% Remove any NaNs and Infs.
if any(~isreal(T(:))), T = abs(T); end
T(isnan(T) | (isinf(T) & T < 0)) = min(T(~isinf(T(:))));
T(isinf(T) & T > 0) = max(T(~isinf(T(:))));

% Compute the vertices.
st = size(T);
sv = st+1;
verts = zeros(prod(sv),3,'single');
verts(:,3) = repmat((.5:1:st(1)+.5).',size(verts,1)/sv(1),1);
verts(:,2) = repmat(kron((.5:1:st(2)+.5).',ones(sv(1),1)),sv(3),1);
verts(:,1) = kron((.5:1:st(3)+.5).',ones(size(verts,1)/sv(3),1));

% Compute the minimal number of colors and faces.
idx = {(1:prod(st+[1 0 0]))',(1:prod(st+[0 1 0]))',(1:prod(st+[0 0 1]))'};
off = [0 cumsum(cellfun(@numel,idx))];
color = zeros(sum(cellfun(@numel,idx)),1,'single');
faces = zeros(sum(cellfun(@numel,idx)),4,'single');
Tint = cat(1,T(1,:,:),max(T(1:end-1,:,:),T(2:end,:,:)),T(end,:,:));
[i,j,k] = ind2sub(size(Tint),idx{1});
color(off(1)+idx{1}) = single(Tint(:));
faces(off(1)+idx{1},:) = [sub2ind(sv,i,j  ,k)   sub2ind(sv,i,j+1,k) ...
                          sub2ind(sv,i,j+1,k+1) sub2ind(sv,i,j  ,k+1)];
Tint = cat(2,T(:,1,:),max(T(:,1:end-1,:),T(:,2:end,:)),T(:,end,:));
[i,j,k] = ind2sub(size(Tint),idx{2});
color(off(2)+idx{2}) = single(Tint(:));
faces(off(2)+idx{2},:) = [sub2ind(sv,i  ,j,k)   sub2ind(sv,i+1,j,k) ...
                          sub2ind(sv,i+1,j,k+1) sub2ind(sv,i  ,j,k+1)];
Tint = cat(3,T(:,:,1),max(T(:,:,1:end-1),T(:,:,2:end)),T(:,:,end));
[i,j,k] = ind2sub(size(Tint),idx{3});
color(off(3)+idx{3}) = single(Tint(:));
faces(off(3)+idx{3},:) = [sub2ind(sv,i  ,j  ,k) sub2ind(sv,i+1,j  ,k) ...
                          sub2ind(sv,i+1,j+1,k) sub2ind(sv,i  ,j+1,k)];

% Compute the opacity for each voxel.
mn = min(T(:));
mx = max(T(:));
alpha = (color-mn)/(mx-mn);
cutoff = 1e-4;

% Set up the figure.
cax = newplot;
if exist('degree','var') || exist('thresh','var')
    if exist('degree','var'), k = degree; else k = 1.0; end
    if exist('thresh','var'), t = thresh; else t = 0.5; end
else
    k = 1.0; t = 0.5;
    set(gcf,'Toolbar','figure');
    zoom off; pan off; rotate3d off; datacursormode off;
    lbl1 = uicontrol('Style','text','Position',[5 25 65 15], ...
                     'HorizontalAlignment','left');
    sld1 = uicontrol('Style','slider','Position',[75 25 120 15], ...
                     'Min',0,'Max',1,'Value',0.5, ...
                     'SliderStep',[0.1 0.25], ...
                     'Callback',{@redraw,'t'}, ...
                     'KeyPressFcn',@(obj,evt)toggle(evt.Key));
    lbl2 = uicontrol('Style','text','Position',[5 5 65 15], ...
                     'HorizontalAlignment','left');
    sld2 = uicontrol('Style','slider','Position',[75 5 120 15], ...
                     'Min',1,'Max',5,'Value',1, ...
                     'SliderStep',[0.25 0.5], ...
                     'Callback',{@redraw,'k'}, ...
                     'KeyPressFcn',@(obj,evt)toggle(evt.Key));
    isVisible = true;
    set(gcf,'KeyPressFcn',@(obj,evt)toggle(evt.Key));
end
redraw();

function toggle(key)
    % Toggle show or hide of uicontrols.
    if strcmpi(key,'h')
        onoff = 'on'; if isVisible, onoff = 'off'; end
        cellfun(@(c)set(c,'Visible',onoff),{lbl1,sld1,lbl2,sld2});
        isVisible = ~isVisible;
    end
end

function redraw(hobj,~,param)

    % Update the parameter.
    if nargin > 1
        switch param
            case 't', t = get(hobj,'Value');
            case 'k', k = get(hobj,'Value');
        end
    end
    if exist('lbl1','var')
        set(lbl1,'String',sprintf('thresh = %g',t));
        set(lbl2,'String',sprintf('degree = %g',k));
    end

    % Speed up rendering a bit.
    plot3(cax,nan,nan,nan,varargin{:});
    set(gcf,'DoubleBuffer','off');
    set(cax,'XLimMode','manual','YLimMode','manual', ...
            'ZLimMode','manual','CLimMode','manual','ALimMode','manual');

    % Set axis properties.
    xlabel('k');
    ylabel('j');
    zlabel('i');
    xlim([.5 st(3)+.5]);
    ylim([.5 st(2)+.5]);
    zlim([.5 st(1)+.5]);
    set(cax,'YDir','reverse');
    set(cax,'ZDir','reverse');
    caxis([mn mx]);

    % Display grid.
    step = max(1,round(st/6));
    set(cax,'XTickMode','manual','YTickMode','manual','ZTickMode','manual');
    tk = [1:step(3):st(3)-1 st(3)];
    if tk(end)-tk(end-1) < .3*step(3), tk = [tk(1:end-2) tk(end)]; end
    set(cax,'XTick',tk);
    tk = [1:step(2):st(2)-1 st(2)];
    if tk(end)-tk(end-1) < .3*step(2), tk = [tk(1:end-2) tk(end)]; end
    set(cax,'YTick',tk);
    tk = [1:step(1):st(1)-1 st(1)];
    if tk(end)-tk(end-1) < .3*step(1), tk = [tk(1:end-2) tk(end)]; end
    set(cax,'ZTick',tk);
    grid on;

    % Draw the voxels.
    opaq = alpha;
    opaq(alpha<t) = t^(1-k)*alpha(alpha<t).^k;
    opaq(opaq > 0.9) = 1; % Fix Matlab render bug.
    patch('Vertices',verts, ...
          'Faces',faces(opaq > cutoff,:), ...
          'EdgeColor','none','FaceAlpha','flat','FaceColor','flat', ...
          'FaceVertexAlphaData',opaq(opaq > cutoff), ...
          'FaceVertexCData',color(opaq > cutoff));

end

end

author: andrewssobral

repo: lrslibrary


text:
function value = calculateValue(s0,price,inFlow,turbFlow,spillFlow,C2A,n)
% Calculate the revenue generated at each time step for given turbine and
% spill flows.
%
% Copyright (c) 2012, MathWorks, Inc.
%%
% Initialize Value
value = zeros(n,1);

% Intialize Storage Vector
s = zeros(n,1);
s(1) = s0;

% Loop through time, calculate storage at each time step and calculate
% revenue from storage and turbine flow
for ii = 2:n
    s(ii) = s(ii-1)+(inFlow(ii-1) - turbFlow(ii-1) - spillFlow(ii-1))*C2A;
    value(ii) = price(ii)*turbFlow(ii)*(0.00001*(s(ii)+s(ii-1))/2+10)/1000;
end

author: Sable

repo: mcbench-benchmarks


text:
function x=nsst_rec1(dst,lpfilt)
% This function performs the inverse (local) nonsubsampled shearlet transform as given
% in G. Easley, D. Labate and W. Lim, "Sparse Directional Image Representations
% using the Discrete Shearlet Transform", Appl. Comput. Harmon. Anal. 25 pp.
% 25-46, (2008).
%
% Input:
%
% dst          - the nonsubsampled shearlet coefficients
%
% lpfilt       - the filter to be used for the Laplacian
%                Pyramid/ATrous decomposition using the codes
%                written by Arthur L. Cunha
%
% Output
%
% x         - the reconstructed image
%
% Code contributors: Glenn R. Easley, Demetrio Labate, and Wang-Q Lim.
% Copyright 2011 by Glenn R. Easley. All Rights Reserved.
%

level=length(dst)-1;
y{1}=dst{1};
for i=1:level,
      y{i+1} = real(sum(dst{i+1},3));
end

x=real(atrousrec(y,lpfilt));

author: thomaskuestner


text:
function [biasFieldCorrected, varargout] = segment(this, varargin)
% Segments brain images using SPM's unified segmentation approach.
% This warps the brain into a standard space and segment it there using tissue
% probability maps in this standard space.
%
% Since good warping and good segmentation are interdependent, this is done
% iteratively until a good tissue segmentation is given by probality maps
% that store how likely a voxel is of a certain tissue type
% (either in native or standard space).
% Furthermore, a deformation field from native to standard space (or back)
% has then been found for warping other images of the same native space.
%
%   Y = MrImage()
%   [biasFieldCorrected, tissueProbMaps, deformationFields, biasField] = ...
%   Y.segment(...
%       'representationIndexArray', representationIndexArray, ...
%       'spmParameterName1', spmParameterValue1, ...
%       ...
%       'spmParameterNameN', spmParameterValueN)
%
% This is a method of class MrImageSpm4D.
%
% NOTE: If an nD image is given, then all dimension > 3  will be treated as
% channels.
%
% IN
%   tissueTypes         cell(1, nTissues) of strings to specify which
%                       tissue types shall be written out:
%                       'GM'    grey matter
%                       'WM'    white matter
%                       'CSF'   cerebrospinal fluid
%                       'bone'  skull and surrounding bones
%                       'fat'   fat and muscle tissue
%                       'air'   air surrounding head
%
%                       default: {'GM', 'WM', 'CSF'}
%
%   mapOutputSpace      'native' (default) or 'warped'/'mni'/'standard'
%                       defines coordinate system in which images shall be
%                       written out;
%                       'native' same space as image that was segmented
%                       'warped' standard Montreal Neurological Institute
%                                (MNI) space used by SPM for unified segmentation
%  deformationFieldDirection determines which deformation field shall be
%                       written out,if any
%                       'none' (default) no deformation fields are stored
%                       'forward' subject => mni (standard) space
%                       'backward'/'inverse' mni => subject space
%                       'both'/'all' = 'forward' and 'backward'
%  saveBiasField        0 (default) or 1
%  biasRegularisation   describes the amount of expected bias field
%                       default: 0.001 (light)
%                       no: 0; extremely heavy: 10
%  biasFWHM             full-width-at-half-maximum of the Gaussian
%                       non-uniformity bias field (in mm)
%                       default: 60 (mm)
%  fileTPM              tissue probablity maps for each tissue class
%                       default: SPM' TPMs in spm/tpm
%  mrfParameter         strenght of the Markov Random Field cleanup
%                       performed on the tissue class images
%                       default: 1
%  cleanUp              crude routine for extracting the brain from
%                       segmented images ('no', 'light', 'thorough')
%                       default: 'light'
%  warpingRegularization regularization for the different terms of the
%                       registration
%                       default: [0 0.001 0.5 0.05 0.2]
%  affineRegularisation regularisation for the initial affine registration
%                       of the image to the tissue probability maps (i.e.
%                       into standard space)
%                       for example, the default ICBM template are slighlty
%                       larger than typical brains, so greater zooms are
%                       likely to be needed
%                       default: ICBM spase template - European brains
%  smoothnessFwhm       fudge factor to account for correlation between
%                       neighbouring voxels (in mm)
%                       default: 0 (for MRI)
%  samplingDistance     approximate distance between sampled points when
%                       estimating the model parameters (in mm)
%                       default: 3
%
%   Parameters for high-dim application:
%
%   representationIndexArray:   a selection (e.g. {'t', 1}) which is then
%                               applied to obtain one nD image, where all
%                               dimensions > 3 are treated as additional
%                               channels
%                               default representationIndexArray: all
%                               dimensions not the imageSpaceDims
%   imageSpaceDims              cell array of three dimLabels defining the
%                               dimensions that define the physical space
%                               the image is in
%                               default imageSpaceDims: {'x', 'y', 'z'}
%   splitComplex                'ri' or 'mp'
%                               If the data are complex numbers, real and
%                               imaginary or magnitude and phase are
%                               realigned separately.
%                               default: mp (magnitude and p)
%                               Typically, realigning the magnitude and
%                               applying it to the phase data makes most
%                               sense; otherwise, using real and imaginary
%                               part, more global phase changes would
%                               impact on estimation
%
% OUT
%   biasCorrected       bias corrected images
%   tissueProbMaps      (optional) cell(nTissues,1) of 3D MrImages
%                       containing the tissue probability maps in the
%                       respective order as volumes,
%   deformationFields   (optional) cell(nDeformationFieldDirections,1)
%                       if deformationFieldDirection is 'both', this cell
%                       contains the forward deformation field in the first
%                       entry, and the backward deformation field in the
%                       second cell entry; otherwise, a cell with only one
%                       element is returned
%   biasField           (optional) bias field
%
% EXAMPLE
% [biasFieldCorrected, tissueProbMaps, deformationFields, biasField] =
%   Y.segment();
%
% for 7T images stronger non-uniformity expected
% [biasFieldCorrected, tissueProbMaps, deformationFields, biasField] = ...
%   m.segment('biasRegularisation', 1e-4, 'biasFWHM', 18, ...
%   'cleanUp', 2, 'samplingDistance', 2);
%
%   See also MrImage spm_preproc MrImageSpm4D.segment
% Author:   Saskia Bollmann & Lars Kasper
% Created:  2019-12-23
% Copyright (C) 2019 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

% spm parameters (details above)
spmDefaults.tissueTypes = {'WM', 'GM', 'CSF'};
spmDefaults.mapOutputSpace = 'native';
spmDefaults.deformationFieldDirection = 'none';
spmDefaults.saveBiasField = 0;
spmDefaults.biasRegularisation = 0.001;
spmDefaults.biasFWHM = 60;
spmDefaults.fileTPM = [];
spmDefaults.mrfParameter = 1;
spmDefaults.cleanUp = 'light';
spmDefaults.warpingRegularization = [0 0.001 0.5 0.05 0.2];
spmDefaults.affineRegularisation = 'mni';
spmDefaults.smoothnessFwhm = 0;
spmDefaults.samplingDistance = 3;

[spmParameters, unusedVarargin] = tapas_uniqc_propval(varargin, spmDefaults);

% for split/apply functionality
methodParameters = {spmParameters};

% use cases: abs of complex, single index on many!
defaults.representationIndexArray   = {{}}; % default: use all
defaults.imageSpaceDims             = {};
defaults.splitComplex               = 'mp';

args = tapas_uniqc_propval(unusedVarargin, defaults);
tapas_uniqc_strip_fields(args);

% set imageSpaceDims
if isempty(imageSpaceDims)
    imageSpaceDims = {'x','y','z'};
end

% check whether real/complex
isReal = isreal(this);

if isReal
    splitComplexImage = this.copyobj();
else
    splitComplexImage = this.split_complex(splitComplex);
end

% prepare output container with right size
varargoutForMrImageSpm4D = cell(numel(representationIndexArray),nargout-1);
hasArgOut = nargout > 1;
hasBiasField = nargout > 3;
for iRepresentation = 1:numel(representationIndexArray)

    % apply selection
    inputSegment = splitComplexImage.select(...
        representationIndexArray{iRepresentation});

    % Merge all n>3 dims, which are not part of the representationIndexArray, into 4D array
    mergeDimLabels = setdiff(inputSegment.dimInfo.dimLabels, imageSpaceDims);
    % additional channels need to be in the t dimensions so they become part of
    % the same nifti file
    % empty mergeDimLabels just return the original object, e.g. for true 3D
    % images
    [mergedImage, newDimLabel] = ...
        inputSegment.merge(mergeDimLabels, 'dimLabels', 't');

    if hasArgOut
        [biasFieldCorrected{iRepresentation}, varargoutForMrImageSpm4D{iRepresentation, :}] = ...
            mergedImage.apply_spm_method_per_4d_split(@segment, ...
            'methodParameters', methodParameters);
    else
        biasFieldCorrected{iRepresentation} = mergedImage.apply_spm_method_per_4d_split(@segment, ...
            'methodParameters', methodParameters);
    end

    % un-do merge operation using combine
    if ~isempty(mergeDimLabels)
        % not necessary for 4D images - just reset dimInfo
        if numel(mergeDimLabels) == 1
            origDimInfo = inputSegment.dimInfo;
            biasFieldCorrected{iRepresentation}.dimInfo = origDimInfo;
            % also for the bias fields
            if hasBiasField
                varargoutForMrImageSpm4D{iRepresentation, 3}.dimInfo = origDimInfo;
            end
        else
            % created original dimInfo per split
            origDimInfo = inputSegment.dimInfo.split(mergeDimLabels);
            % un-do reshape
            split_array = biasFieldCorrected{iRepresentation}.split('splitDims', newDimLabel);
            split_array = reshape(split_array, size(origDimInfo));
            % add original dimInfo
            for nSplits = 1:numel(split_array)
                split_array{nSplits}.dimInfo = origDimInfo{nSplits};
            end
            % and combine
            biasFieldCorrected{iRepresentation} = split_array{1}.combine(split_array);

            % same for the bias fields
            if hasBiasField
                % un-do reshape
                split_array = varargoutForMrImageSpm4D{iRepresentation, 3}{1}.split('splitDims', newDimLabel);
                split_array = reshape(split_array, size(origDimInfo));
                % add original dimInfo
                for nSplits = 1:numel(split_array)
                    split_array{nSplits}.dimInfo = origDimInfo{nSplits};
                end
                varargoutForMrImageSpm4D{iRepresentation, 3} = {split_array{1}.combine(split_array)};
            end
        end
    end

    if ~isReal
        % un-do complex split
        biasFieldCorrected{iRepresentation} = biasFieldCorrected{iRepresentation}.combine_complex();
    end

    % add representation index back to TPMs and deformation field to
    % combine them later
    if ~isempty(representationIndexArray{iRepresentation})
        for nOut = 1:nargout-2
            addDim = varargoutForMrImageSpm4D{iRepresentation, nOut}{1}.dimInfo.nDims + 1;
            dimLabels = representationIndexArray{iRepresentation}{1};
            % only pick the first one
            samplingPoints = representationIndexArray{iRepresentation}{2}(1);
            for nClasses = 1:numel(varargoutForMrImageSpm4D{iRepresentation, nOut})
                varargoutForMrImageSpm4D{iRepresentation, nOut}{nClasses}.dimInfo.add_dims(...
                    addDim, 'dimLabels', dimLabels, ...
                    'samplingPoints', samplingPoints);
            end
        end
    end
end

% combine bias field corrected
biasFieldCorrected = biasFieldCorrected{1}.combine(biasFieldCorrected);

% combine varargout
for nOut = 1:nargout-1
    toBeCombined = varargoutForMrImageSpm4D(:, nOut);
    % the TPMs are cell-arrays of images, so we need to combine them per
    % tissue class
    for nClasses = 1:numel(toBeCombined{1})
        for nCells = 1:size(toBeCombined, 1)
            thisCombine(nCells) = toBeCombined{nCells}(nClasses);
        end
        combinedImage{nClasses, 1} = thisCombine{1}.combine(thisCombine);
    end
    varargout{1, nOut} = combinedImage;
    clear toBeCombined thisCombine combinedImage;
end

end

author: translationalneuromodeling

repo: tapas


text:
% Jiao Xianjun (putaoshu@msn.com; putaoshu@gmail.com)
% multi_rtl_sdr_split_scanner.m
% Frequency band scanning via multiple rtl-sdr dongles. Each dongle for each sub-band to speedup!
% A script of project: https://github.com/JiaoXianjun/multi-rtl-sdr-calibration

% Assume that you have installed rtl-sdr
% (http://sdr.osmocom.org/trac/wiki/rtl-sdr) and have those native utilities run correctly already.

% For example, you have multiple dongles, please run multiple rtl_tcp in multiple shell respectively as:
% rtl_tcp -p 1234 -d 0
% rtl_tcp -p 1235 -d 1
% rtl_tcp -p 1236 -d 2
% ...

% Then run this script in MATLAB.

% ATTENTION! In some computer, every time before you run this script, maybe you need to terminate multiple rtl_tcp and re-launch them again.
% ATTENTION! Please reduce number of inspected points by reducing frequency range or increasing step size, if your computer hasn't enough memory installed. Because all signlas are stored firstly before processing.

% Change following parameters as you need:

% Number of dongles you have connected to your computer
num_dongle = 1; % more dongles, much faster.

% Beginning of the band you are interested in
start_freq = 502e6; % for test
% start_freq = 935e6; % Beginning of Primary GSM-900 Band downlink
% start_freq = (1575.42-15)*1e6; % GPS L1
% start_freq = (1207.14-30)*1e6; % COMPASS B2I
% start_freq = (1561.098-30)*1e6; % COMPASS B1I
% start_freq = 1.14e9; % Beginning of GNSS(GPS/GLONASS/COMPASS/Galileo) Band

% End of the band you are interested in
end_freq = 702e6; % for test
% end_freq = 960e6; % End of Primary GSM-900 Band downlink
% end_freq = (1575.42+30)*1e6; % GPS L1
% end_freq = (1207.14+30)*1e6; % COMPASS B2I
% end_freq = (1561.098+30)*1e6; % COMPASS B1I
% end_freq = 1.63e9; % Beginning of GNSS(GPS/GLONASS/COMPASS/Galileo) Band

freq_step = 0.5e6; % less step, higher resolution, narrower FIR bandwidth, slower speed

observe_time = 0.2; % observation time at each frequency point. ensure it can capture your signal!

RBW = freq_step; % Resolution Bandwidth each time we inspect

gain = 20; % If this is larger than 0, the fixed gain will be set to dongles

% use high sampling rate and FIR to improve estimation accuracy
sample_rate = 2.048e6; % sampling rate of dongles

coef_order = (2^(ceil(log2(sample_rate/RBW))))-1;
coef_order = min(coef_order, 127);
coef_order = max(coef_order, 31);
coef = fir1(coef_order, RBW/sample_rate);
% freqz(coef, 1, 1024);

num_samples = observe_time*sample_rate;

clf;
close all;

% construct freq set for each dongle
freq_orig = start_freq:freq_step:end_freq;
num_freq_per_sub_band = ceil(length(freq_orig)/num_dongle);
num_pad = num_freq_per_sub_band*num_dongle - length(freq_orig);
freq = [freq_orig freq_orig(end)+(1:num_pad).*freq_step];
freq = vec2mat(freq, num_freq_per_sub_band);

real_count = zeros(1, num_dongle);
s_all = uint8( zeros(2*num_samples, num_freq_per_sub_band*num_dongle) );
decimate_ratio = floor(sample_rate/(2*RBW));

% check if previous tce objects existed. if so clear them
if ~isempty(who('tcp_obj'))
    for i=1:length(tcp_obj)
        fclose(tcp_obj{i});
        delete(tcp_obj{i});
    end
    clear tcp_obj;
end

% construct tcp objects
tcp_obj = cell(1, num_dongle);
for i=1:num_dongle
    tcp_obj{i} = tcpip('127.0.0.1', 1233+i); % for dongle i
end

% set some parameters to tcp objects, and open them.
for i=1:num_dongle
    set(tcp_obj{i}, 'InputBufferSize', 8*2*num_samples);
    set(tcp_obj{i}, 'Timeout', 60);
end
for i=1:num_dongle
    fopen(tcp_obj{i});
end

% set gain
for i=1:num_dongle
    set_gain_tcp(tcp_obj{i}, gain*10); %be careful, in rtl_sdr the 10x is done inside C program, but in rtl_tcp the 10x has to be done here.
end

% set sampling rate
for i=1:num_dongle
    set_rate_tcp(tcp_obj{i}, sample_rate);
end

% set different start freq to different dongle
for i=1:num_dongle
    set_freq_tcp(tcp_obj{i}, freq(i,1));
end

% read and discard to flush
for i=1:num_dongle
    fread(tcp_obj{i}, 8*2*num_samples, 'uint8');
end

% capture samples of all frequencies firstly!
tic;
for freq_idx = 1:num_freq_per_sub_band
    while 1 % read data at current frequency until success
        for i=1:num_dongle
            set_freq_tcp(tcp_obj{i}, freq(i, freq_idx)); % set different frequency to different dongle
        end
        for i=1:num_dongle
            [s_all(:, freq_idx + (i-1)*num_freq_per_sub_band), real_count(i)] = fread(tcp_obj{i}, 2*num_samples, 'uint8'); % gather data from different dongles to s_all
        end

        if sum(real_count-(2*num_samples)) ~= 0
            disp(num2str([idx 2*num_samples, real_count]));
        else
            break;
        end
    end
end
e = toc;
ideal_time_cost = observe_time*num_freq_per_sub_band;

% close TCP
for i=1:num_dongle
    fclose(tcp_obj{i});
end
for i=1:num_dongle
    delete(tcp_obj{i});
end
clear tcp_obj;

disp('Scanning done!');
disp(['actual time cost ' num2str(e) ' ideal cost ' num2str(ideal_time_cost) ' efficiency(ideal/actual) ' num2str(ideal_time_cost/e)]);
disp(' ');
disp('Begin process ...');

% generate power spectrum
tic;
r = raw2iq( double( s_all ) ); % remove DC. complex number constructed.
r_flt = filter(coef, 1, r);% filter target band out
power_spectrum = mean(abs(r_flt(1:decimate_ratio:end, :)).^2, 1);% get averaged power
e1 = toc;
disp(['time cost ' num2str(e1) ' scan/process ' num2str(e/e1)]);
disp(['total time cost ' num2str(e1+e)]);

% plot power spectrum (converted to dB)
freq_linear = freq';
freq_linear = freq_linear(:)';
figure;
format_string = {'b.-', 'r.-', 'k.-', 'm.-'};
for i=1:num_dongle
    plot(freq(i,:).*1e-6, 10.*log10(power_spectrum( (i-1)*num_freq_per_sub_band+1: i*num_freq_per_sub_band)), format_string{i}); hold on;
end

legend_string = cell(1, num_dongle);
for i=1:num_dongle
    legend_string{i} = ['dongle ' num2str(i)];
end
legend(legend_string);

filename = ['split_scan_' num2str(start_freq) '_' num2str(end_freq) '_gain' num2str(gain) '_' num2str(num_dongle) 'dongles.mat'];
save(filename, 'power_spectrum', 'start_freq', 'end_freq', 'freq_step', 'observe_time', 'RBW', 'gain', 'sample_rate', 'coef', 'freq');

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
classdef MimImportNewManualSegmentation < MimGuiPlugin
    % MimImportNewManualSegmentation. Gui Plugin
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimImportNewManualSegmentation 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 = 'Import new segmentation'
        SelectedText = 'Import new segmentation'
        ToolTip = 'Imports a new segmentation from an external file'
        Category = 'Segment region'
        Visibility = 'Dataset'
        Mode = 'Segment'

        Icon = 'add_patch.png'
        Location = 26

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImportManualSegmentation();
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded;
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


% pivot_input.m // Jon Lee
% data for pivoting example

A = [1 2 1 0 0 0; 3 1 0 1 0 0; 1.5 1.5 0 0 1 0; 0 1 0 0 0 1];
c = [6 7 -2 0 4 4.5]';
b = [7 9 6 3.3]';
beta = [1,2,4,6];
[m,n] = size(A);
eta = setdiff(1:n,beta); % lazy eta initialization


text:
function y = nanvar(x,w,dim)
%Replacement for Matlab NANVAR Variance function, ignoring NaNs.
%
%
%
%

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

sz = size(x);
if nargin < 3 || isempty(dim)
    % The output size for [] is a special case when DIM is not given.
    if isequal(x,[]), y = NaN(class(x)); return; end

    % Figure out which dimension sum will work along.
    dim = find(sz ~= 1, 1);
    if isempty(dim), dim = 1; end
elseif dim > length(sz)
    sz(end+1:dim) = 1;
end

% Need to tile the mean of X to center it.
tile = ones(size(sz));
tile(dim) = sz(dim);

if isequal(w,0) || isequal(w,1)
    % Count up non-NaNs.
    n = sum(~isnan(x),dim);

    if w == 0
        % The unbiased estimator: divide by (n-1).  Can't do this when
        % n == 0 or 1, so n==1 => we'll return zeros
        denom = max(n-1, 1);
    else
        % The biased estimator: divide by n.
        denom = n; % n==1 => we'll return zeros
    end
    denom(n==0) = NaN; % Make all NaNs return NaN, without a divideByZero warning

    x0 = x - repmat(nanmean(x, dim), tile);
    y = nansum(abs(x0).^2, dim) ./ denom; % abs guarantees a real result

% Weighted variance
elseif numel(w) ~= sz(dim)
    error('MATLAB:nanvar:InvalidSizeWgts','The length of W must be compatible with X.');
elseif ~(isvector(w) && all(w(~isnan(w)) >= 0))
    error('MATLAB:nanvar:InvalidWgts','W must be a vector of nonnegative weights, or a scalar 0 or 1.');
else
    % Embed W in the right number of dims.  Then replicate it out along the
    % non-working dims to match X's size.
    wresize = ones(size(sz)); wresize(dim) = sz(dim);
    wtile = sz; wtile(dim) = 1;
    w = repmat(reshape(w, wresize), wtile);

    % Count up non-NaNs.
    n = nansum(~isnan(x).*w,dim);

    x0 = x - repmat(nansum(w.*x, dim) ./ n, tile);
    y = nansum(w .* abs(x0).^2, dim) ./ n; % abs guarantees a real result
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function data = h5loadgroup(loc)

data = struct();

num_objs = H5G.get_num_objs(loc);

%
% Load groups and datasets
%
for j_item = 0:num_objs-1
    objtype = H5G.get_objtype_by_idx(loc, j_item);
    objname = H5G.get_objname_by_idx(loc, j_item);

    if objtype == 0
        % Group
        name = regexprep(objname, '.*/', '');

        if ~isempty(regexp(name,'^[a-zA-Z].*'))
            group_loc = H5G.open(loc, name);
            try
                sub_data = h5loadgroup(group_loc);
                H5G.close(group_loc);
            catch exc
                H5G.close(group_loc);
                rethrow(exc);
            end
            if isempty(path_parts)
                data = setfield(data, name, sub_data);
            else
                data = sub_data;
                return
            end
        end

    elseif objtype == 1
        % Dataset
        name = regexprep(objname, '.*/', '');
        if ~isempty(regexp(name,'^[a-zA-Z].*'))
            dataset_loc = H5D.open(loc, name);
            try
                % NOTE: HDF5 stores contiguous muti-dimensional arrays in row-major order.
                % Matlab stores them in row-major order. We want to transpose the loaded data
                % it back to Matlab's column-major storage order and thus get back the
                % original array.
                sub_data = H5D.read(dataset_loc, 'H5ML_DEFAULT', 'H5S_ALL','H5S_ALL','H5P_DEFAULT');
                sub_data = HDF5_PostProcessing(sub_data);

                H5D.close(dataset_loc);
            catch exc
                H5D.close(dataset_loc);
                rethrow(exc);
            end

            sub_data = fix_data(sub_data);

            data = setfield(data, name, sub_data);
        end
    end

end

% ---------------------------------------------------------------------
function data = fix_data(data)
% Fix some common types of data to more friendly form.

if isstruct(data)
    fields = fieldnames(data);
    if length(fields) == 2 & strcmp(fields{1}, 'r') & strcmp(fields{2}, 'i')
        if isnumeric(data.r) & isnumeric(data.i)
            data = data.r + 1j*data.i;
        end
    end
end

if isnumeric(data) & ndims(data) > 1
    % permute dimensions
    data = permute(data, fliplr(1:ndims(data)));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% load dataset
load('SimTest.mat');
% current data
steering_req_vec = debug.debug_mvdc_curvvel_tracking_debug_Req_Delta_rad.Data;
steering_vec = debug.debug_VehicleSensorData_Delta_Wheel_rad.Data;
kappa_vec = debug.debug_mvdc_state_estimation_debug_StateEstimate_kappa_radpm.Data;
velocity_vec = debug.debug_mvdc_state_estimation_debug_StateEstimate_v_mps.Data;
curv_req = debug.debug_mvdc_curvvel_tracking_debug_Curv_TargetCurv_radpm.Data;

% parameters
filter_coeff_virtual_ss = 0.07;
P_VDC_SCLearn_Sigma = 0.01;
P_VDC_SCLearn_LengthScales = [0.04, 5];
l_front_m = 1.51;
l_rear_m = 1.388;
P_VDC_SCLearn_FilterCoeff = 0.999;
P_VDC_SCLearn_KappaMax_radpm = 0.13;
P_VDC_SCLearn_vMax_mps = 60;
P_VDC_MinVelSlipCalc_mps = 3;

% iterate via all data points
steering_steady_virtual = zeros(length(steering_req_vec)+1, 1);
steering_comp = zeros(length(steering_req_vec), 1);
% clear persistent variaables
clear learnSteeringCharacteristic;
for i = 1:1:length(steering_req_vec)
    % call steering characteristic update
    [SC_Kappa_Out, SC_Vel_Out, SC_DeltaNeutral_Out, SC_meas] =...
        learnSteeringCharacteristic(steering_steady_virtual(i), ...
                                        kappa_vec(i), velocity_vec(i), P_VDC_SCLearn_Sigma, ...
                            P_VDC_SCLearn_LengthScales, l_front_m, l_rear_m, ...
                            P_VDC_SCLearn_FilterCoeff, P_VDC_SCLearn_KappaMax_radpm, ...
                            P_VDC_SCLearn_vMax_mps, P_VDC_MinVelSlipCalc_mps);
    steering_steady_virtual(i+1) = (1-filter_coeff_virtual_ss)*steering_steady_virtual(i) + ...
                                    filter_coeff_virtual_ss*steering_req_vec(i);
    steering_comp(i) = interp2(SC_Vel_Out, SC_Kappa_Out, SC_DeltaNeutral_Out,...
                            velocity_vec(i), curv_req(i),...
                            'linear', 0);
end

%% visualize Results
close all;

figure;
surface(SC_Vel_Out, SC_Kappa_Out, SC_DeltaNeutral_Out);
colorbar;
xlabel('Velocity in mps');
ylabel('Kappa in radpm');
zlabel('Steering Angle Correction in rad');
title('Smoothed Understeer Compensation Profile');

figure;
surface(SC_Vel_Out, SC_Kappa_Out, SC_meas);
colorbar;
xlabel('Velocity in mps');
ylabel('Kappa in radpm');
zlabel('Steering Angle Correction in rad');
title('Raw Understeer Compensation Profile');

% visualize steering time values
figure;
plot(steering_vec);
grid on; hold on;
plot(kappa_vec*2.898);
plot(steering_req_vec);
plot(steering_steady_virtual);
plot(steering_comp);
legend('steering', 'neutral', 'req', 'virtual_ss', 'ff');

author: TUMFTM


text:
function c = l1pp_cholesky ( n, h )

%*****************************************************************************80
%
%% L1PP_CHOLESKY computes the Cholesky factor of the 1D PP Laplacian.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%    N must be at least 3.
%
%    Input, real H, the spacing between points.
%
%    Output, real C(N,N), the Cholesky factor.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'L1PP_CHOLESKY - Fatal error!\n' );
    fprintf ( 1, '  N < 3.\n' );
    error ( 'L1PP_CHOLESKY - Fatal error!' );
  end

  c = zeros ( n, n );

  for i = 1 : n - 1
    c(i,i) = sqrt ( i + 1.0 ) / sqrt ( i );
  end

  for i = 1 : n - 2
    c(i,i+1) = - i / ( i + 1.0 ) * sqrt ( i + 1.0 ) / sqrt ( i );
  end

  for i = 1 : n - 2
    c(i,n) = - 1.0 / ( i + 1.0 ) * sqrt ( i + 1.0 ) / sqrt ( i );
  end

  i = n - 1;
  c(i,n) = - n / ( i + 1.0 ) * sqrt ( i + 1.0 ) / sqrt ( i );

  c(1:n,1:n) = c(1:n,1:n) / h;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ R ] = vl_nnr2R( r,dCdR )
%r2R Axis angle to rotation matrix layer
% Forwards mode:
% R = vl_nnr2R(r);
%   r is of size 1 x 1 x 3 x nbatch containing axis angle vectors
%   R is of size 1 x 3 x 3 x nbatch containing rotation matrices
% Backwards mode:
% dCdr = vl_nnr2R(r,dCdR);
%   dCdR is derivative of the cost function with respect to R and is of
%   size 1 x 3 x 3 x nbatch
%   dCdr is derivative of the cost function with respect to r and is of
%   size 1 x 1 x 3 x nbatch
%
% Useful reference for formula:
% [1] Gallego, Guillermo, and Anthony Yezzi. "A compact formula for the
%     derivative of a 3-D rotation in exponential coordinates." Journal
%     of Mathematical Imaging and Vision 51.3 (2015): 378-384.

nbatch = size(r,4);
if nargin<2
    % Forwards

    % Vectors with zero length are a special case - they just have an
    % identity rotation matrix
    idmask = squeeze(all(r==0,3));

    R = zeros(1,3,3,nbatch,'single');

    % Fill in the identity cases
    R(1,1,1,idmask)=1;
    R(1,2,2,idmask)=1;
    R(1,3,3,idmask)=1;

    theta = sqrt(sum(r(:,:,:,~idmask).^2,3)); % 1 x 1 x 1 x non-zero nbatch rotation angles
    k1 = r(:,:,1,~idmask)./theta; % rotation axis unit vectors
    k2 = r(:,:,2,~idmask)./theta; % rotation axis unit vectors
    k3 = r(:,:,3,~idmask)./theta; % rotation axis unit vectors

    % Fill in the non-identity cases using the following formula (see [1]):
    % K = [0 -k(3) k(2); k(3) 0 -k(1); -k(2) k(1) 0];
    % R = eye(3) + sin(theta).*K + (1-cos(theta)).*K*K;

    R(1,1,1,~idmask) = (cos(theta) - 1).*k2.^2 + (cos(theta) - 1).*k3.^2 + 1;
    R(1,1,2,~idmask) =  -k3.*sin(theta) - k1.*k2.*(cos(theta) - 1);
    R(1,1,3,~idmask) = k2.*sin(theta) - k1.*k3.*(cos(theta) - 1);
    R(1,2,1,~idmask) = k3.*sin(theta) - k1.*k2.*(cos(theta) - 1);
    R(1,2,2,~idmask) = (cos(theta) - 1).*k1.^2 + (cos(theta) - 1).*k3.^2 + 1;
    R(1,2,3,~idmask) = -k1.*sin(theta) - k2.*k3.*(cos(theta) - 1);
    R(1,3,1,~idmask) = -k2.*sin(theta) - k1.*k3.*(cos(theta) - 1);
    R(1,3,2,~idmask) = k1.*sin(theta) - k2.*k3.*(cos(theta) - 1);
    R(1,3,3,~idmask) = (cos(theta) - 1).*k1.^2 + (cos(theta) - 1).*k2.^2 + 1;
else
    % Backwards

	idmask = squeeze(all(r==0,3));

    % Handle identity case
    dRdr1 = zeros(size(dCdR));
    dRdr1(1,2,3,idmask)=-1;
    dRdr1(1,3,2,idmask)=1;
    dRdr2 = zeros(size(dCdR));
    dRdr2(1,1,3,idmask)=1;
    dRdr2(1,3,1,idmask)=-1;
    dRdr3 = zeros(size(dCdR));
    dRdr3(1,1,2,idmask)=-1;
    dRdr3(1,2,1,idmask)=1;

    % Call forwards version to get R
    [ R ] = vl_nnr2R( r );

    % Subselect and precompute some values
    R11 = R(1,1,1,~idmask); R12 = R(1,1,2,~idmask); R13 = R(1,1,3,~idmask);
    R21 = R(1,2,1,~idmask); R22 = R(1,2,2,~idmask); R23 = R(1,2,3,~idmask);
    R31 = R(1,3,1,~idmask); R32 = R(1,3,2,~idmask); R33 = R(1,3,3,~idmask);
    r1 = r(1,1,1,~idmask); r2 = r(1,1,2,~idmask); r3 = r(1,1,3,~idmask);
    sqnorm = (r1.^2 + r2.^2 + r3.^2);

    % For all other cases, find 3x3 derivatives of R with respect to each
    % element of r
    dRdr1(1,1,1,~idmask) = (R31.*(R31.*r1 + r1.*r2 - r3.*(R11 - 1)))./sqnorm - (R21.*(r1.*r3 - R21.*r1 + r2.*(R11 - 1)))./sqnorm;
    dRdr1(1,1,2,~idmask) = (R32.*(R31.*r1 + r1.*r2 - r3.*(R11 - 1)))./sqnorm - (R22.*(r1.*r3 - R21.*r1 + r2.*(R11 - 1)))./sqnorm;
    dRdr1(1,1,3,~idmask) = (R33.*(R31.*r1 + r1.*r2 - r3.*(R11 - 1)))./sqnorm - (R23.*(r1.*r3 - R21.*r1 + r2.*(R11 - 1)))./sqnorm;
    dRdr1(1,2,1,~idmask) = (R11.*(r1.*r3 - R21.*r1 + r2.*(R11 - 1)))./sqnorm - (R31.*(r1.^2 + R21.*r3 - R31.*r2))./sqnorm;
    dRdr1(1,2,2,~idmask) = (R12.*(r1.*r3 - R21.*r1 + r2.*(R11 - 1)))./sqnorm - (R32.*(r1.^2 + R21.*r3 - R31.*r2))./sqnorm;
    dRdr1(1,2,3,~idmask) = (R13.*(r1.*r3 - R21.*r1 + r2.*(R11 - 1)))./sqnorm - (R33.*(r1.^2 + R21.*r3 - R31.*r2))./sqnorm;
    dRdr1(1,3,1,~idmask) = (R21.*(r1.^2 + R21.*r3 - R31.*r2))./sqnorm - (R11.*(R31.*r1 + r1.*r2 - r3.*(R11 - 1)))./sqnorm;
    dRdr1(1,3,2,~idmask) = (R22.*(r1.^2 + R21.*r3 - R31.*r2))./sqnorm - (R12.*(R31.*r1 + r1.*r2 - r3.*(R11 - 1)))./sqnorm;
    dRdr1(1,3,3,~idmask) = (R23.*(r1.^2 + R21.*r3 - R31.*r2))./sqnorm - (R13.*(R31.*r1 + r1.*r2 - r3.*(R11 - 1)))./sqnorm;

    dRdr2(1,1,1,~idmask) = (R31.*(r2.^2 - R12.*r3 + R32.*r1))./sqnorm - (R21.*(R12.*r2 + r2.*r3 - r1.*(R22 - 1)))./sqnorm;
    dRdr2(1,1,2,~idmask) = (R32.*(r2.^2 - R12.*r3 + R32.*r1))./sqnorm - (R22.*(R12.*r2 + r2.*r3 - r1.*(R22 - 1)))./sqnorm;
    dRdr2(1,1,3,~idmask) = (R33.*(r2.^2 - R12.*r3 + R32.*r1))./sqnorm - (R23.*(R12.*r2 + r2.*r3 - r1.*(R22 - 1)))./sqnorm;
    dRdr2(1,2,1,~idmask) = (R11.*(R12.*r2 + r2.*r3 - r1.*(R22 - 1)))./sqnorm - (R31.*(r1.*r2 - R32.*r2 + r3.*(R22 - 1)))./sqnorm;
    dRdr2(1,2,2,~idmask) = (R12.*(R12.*r2 + r2.*r3 - r1.*(R22 - 1)))./sqnorm - (R32.*(r1.*r2 - R32.*r2 + r3.*(R22 - 1)))./sqnorm;
    dRdr2(1,2,3,~idmask) = (R13.*(R12.*r2 + r2.*r3 - r1.*(R22 - 1)))./sqnorm - (R33.*(r1.*r2 - R32.*r2 + r3.*(R22 - 1)))./sqnorm;
    dRdr2(1,3,1,~idmask) = (R21.*(r1.*r2 - R32.*r2 + r3.*(R22 - 1)))./sqnorm - (R11.*(r2.^2 - R12.*r3 + R32.*r1))./sqnorm;
    dRdr2(1,3,2,~idmask) = (R22.*(r1.*r2 - R32.*r2 + r3.*(R22 - 1)))./sqnorm - (R12.*(r2.^2 - R12.*r3 + R32.*r1))./sqnorm;
    dRdr2(1,3,3,~idmask) = (R23.*(r1.*r2 - R32.*r2 + r3.*(R22 - 1)))./sqnorm - (R13.*(r2.^2 - R12.*r3 + R32.*r1))./sqnorm;

    dRdr3(1,1,1,~idmask) = (R31.*(r2.*r3 - R13.*r3 + r1.*(R33 - 1)))./sqnorm - (R21.*(r3.^2 + R13.*r2 - R23.*r1))./sqnorm;
    dRdr3(1,1,2,~idmask) = (R32.*(r2.*r3 - R13.*r3 + r1.*(R33 - 1)))./sqnorm - (R22.*(r3.^2 + R13.*r2 - R23.*r1))./sqnorm;
    dRdr3(1,1,3,~idmask) = (R33.*(r2.*r3 - R13.*r3 + r1.*(R33 - 1)))./sqnorm - (R23.*(r3.^2 + R13.*r2 - R23.*r1))./sqnorm;
	dRdr3(1,2,1,~idmask) = (R11.*(r3.^2 + R13.*r2 - R23.*r1))./sqnorm - (R31.*(R23.*r3 + r1.*r3 - r2.*(R33 - 1)))./sqnorm;
    dRdr3(1,2,2,~idmask) = (R12.*(r3.^2 + R13.*r2 - R23.*r1))./sqnorm - (R32.*(R23.*r3 + r1.*r3 - r2.*(R33 - 1)))./sqnorm;
    dRdr3(1,2,3,~idmask) = (R13.*(r3.^2 + R13.*r2 - R23.*r1))./sqnorm - (R33.*(R23.*r3 + r1.*r3 - r2.*(R33 - 1)))./sqnorm;
    dRdr3(1,3,1,~idmask) = (R21.*(R23.*r3 + r1.*r3 - r2.*(R33 - 1)))./sqnorm - (R11.*(r2.*r3 - R13.*r3 + r1.*(R33 - 1)))./sqnorm;
    dRdr3(1,3,2,~idmask) = (R22.*(R23.*r3 + r1.*r3 - r2.*(R33 - 1)))./sqnorm - (R12.*(r2.*r3 - R13.*r3 + r1.*(R33 - 1)))./sqnorm;
    dRdr3(1,3,3,~idmask) = (R23.*(R23.*r3 + r1.*r3 - r2.*(R33 - 1)))./sqnorm - (R13.*(r2.*r3 - R13.*r3 + r1.*(R33 - 1)))./sqnorm;

    dCdr(1,1,1,:) = sum(sum(dCdR.*dRdr1,2),3);
    dCdr(1,1,2,:) = sum(sum(dCdR.*dRdr2,2),3);
    dCdr(1,1,3,:) = sum(sum(dCdR.*dRdr3,2),3);

    R = dCdr;
end

end

author: anilbas