text:
function options = IRset(varargin)
%IRset Set options for IR Tools functions
%
% options = IRset('param1',value1,'param2',value2,...)
%
% Create/alter options structure for methods in IR Tools, in which the
% named parameters have the specified values.  Any unspecified parameters
% are set to [] (they indicate to use the default value for that parameter
% when passed to an IR function). It is sufficient to type only the
% leading characters that uniquely identify the parameter, and  case is
% ignored for parameter names.
% NOTE: For values that are strings, the complete string is required.
%
% options = IRset(OLDOPTS,'param1',value1,...) creates a copy of OLDOPTS
%   with the named parameters altered with the specified values.
%
% options = IRset(OLDOPTS,NEWOPTS) combines an existing options structure
%   OLDOPTS with a new structure NEWOPTS.  Any parameters in NEWOPTS with
%   non-empty values overwrite the corresponding old parameters in OLDOPTS.
%
% options = IRset('IRmethod') creates an options structure with all the
% parameter names and default values relevant to an IR method. That is,
%           IRset('IRmethod')
% or
%           IRset(@IRmethod)
% returns an options structure containing all the parameter names and
% default values relevant to an IR method.
%
% Examples:
%   To create options with the default options for an IR nethod
%       options = IRset('IRmethod');
%   To change the maximum iterations to 150 in options
%       options = IRset(options,'MaxIter',150);
%
% See also: IRget, PRget, PRset

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD License. A separate license file should be provided as part
% of the package.

% Create a struct of all the fields. Note that if new options are added
% then IRoptions_fields.m should be edited to include the field name in
% the list created in IRoptions_fields.m.
%
allfields = IRoption_fields;

% Create cell array.
structinput = cell(2,length(allfields));
% Fields go in first row.
structinput(1,:) = allfields';
% []'s go in second row.
structinput(2,:) = {[]};
% Turn it into correctly ordered comma separated list and call struct
options = struct(structinput{:});

numberargs = nargin; % We might change this value, so assign it.

% If we pass in a function name then return the defaults.
if (numberargs==1) && (ischar(varargin{1}) || isa(varargin{1},'function_handle') )
  if ischar(varargin{1})
    funcname = lower(varargin{1});
  elseif isa(varargin{1},'function_handle')
    funcname = func2str(varargin{1});
  end
  if ~exist(funcname, 'file')
     error('Undefined function.  Please use IRset(''funcname'') or IRset(@funcname), where funcname is a IR method')
  end
  try
    optionsfcn = feval(varargin{1},'defaults');
  catch
    error('IRset ONLY works with an IR method.  Please use IRset(@IR).')
  end
  % The defaults from IR functions don't include all the fields typically,
  % so run the rest of IRset as if called with IRset(options,optionsfcn)
  % to get all the fields.
  varargin{1} = options;
  varargin{2} = optionsfcn;
  numberargs = 2;
end

Names = allfields;
m = size(Names,1);
names = lower(Names);

i = 1;
while i <= numberargs
  arg = varargin{i};
  if ischar(arg)                         % arg is an option name.
    break;
  end
  if ~isempty(arg)                      % [] is a valid options argument.
    if ~isa(arg,'struct')
      error(['Expected argument %d to be a string parameter name ' ...
             'or an options structure \n created with IRset.'], i);
    end
    for j = 1:m
      if any(strcmp(fieldnames(arg),Names{j,:}))
        val = arg.(Names{j,:});
      else
        val = [];
      end
      if ~isempty(val)
        if ischar(val)
          val = lower(deblank(val));
        end
        checkfield(Names{j,:},val)
        options.(Names{j,:}) = val;
      end
    end
  end
  i = i + 1;
end

% A finite state machine to parse name-value pairs.
if rem(numberargs-i+1,2) ~= 0
  error('Arguments must occur in name-value pairs.');
end

expectval = 0;                    % Start expecting a name, not a value.
while i <= numberargs
  arg = varargin{i};

  if ~expectval
    if ~ischar(arg)
      error('Expected argument %d to be a string parameter name.', i);
    end

    lowArg = lower(arg);
    j = find(strcmp(lowArg,names));
    if isempty(j)                       % If no matches.
      error('Invalid parameter name ''%s'' ', arg);
    elseif length(j) > 1                % If more than one match.
      % Check for any exact matches (in case any names are subsets of others).
      k = find(strcmp(lowArg,names));
      if length(k) == 1
        j = k;
      else
        error('Ambiguous parameter name ''%s'' ', arg);
      end
    end
    expectval = 1;                      % We expect a value next.

  else
    if ischar(arg)
      arg = lower(deblank(arg));
    end
    checkfield(Names{j,:},arg);
    options.(Names{j,:}) = arg;
    expectval = 0;
  end
  i = i + 1;
end

if expectval
  error( 'Expected value for parameter ''%s''.', arg);
end

%----SUBFUNCTIONS---------------------------------------------

function checkfield(field,value)
%CHECKFIELD Check validity of structure field contents.
%   CHECKFIELD('field',V) checks the contents of the specified
%   value V to be valid for the field 'field'.

% Empty matrix is always valid.
if isempty(value)
  return
end

% See if 'field' is a valid field.
%validfield = true;
switch field
  case {'MaxIter'} % real positive integer
    [validvalue, errmsg] = PosInteger(field,value);
  case {'MaxIterIn'} % real positive integer
    [validvalue, errmsg] = PosInteger(field,value);
  case {'MaxIterOut'} % real positive integer
    [validvalue, errmsg] = PosInteger(field,value);
  case {'InProjIt'} % real positive integer
    [validvalue, errmsg] = PosInteger(field,value);
  case{'TotIterMax'}
    [validvalue, errmsg] = PosInteger2(field, value);
  case{'ktotcount'}
    [validvalue, errmsg] = nonNegInteger2(field, value);
  case {'x0'} % numeric array, none
    [validvalue, errmsg] = x0type(field,value);
  case {'IterBar'} % off, on
    [validvalue, errmsg] = onOffType(field,value);
  case {'x_true'} % numeric array, off
    [validvalue, errmsg] = x_truetype(field,value);
  case {'AllX'} %  off, on
    [validvalue, errmsg] = onOffType(field,value);
  case {'NoStop'} %  off, on
    [validvalue, errmsg] = onOffType(field,value);
  case {'NoStopIn'} %  off, on
    [validvalue, errmsg] = onOffType(field,value);
  case {'NoStopOut'} %  off, on
    [validvalue, errmsg] = onOffType(field,value);
  case {'restart'} %  off, on
    [validvalue, errmsg] = onOffType(field,value);
  case {'MatUpdate'}
    [validvalue, errmsg] = onOffType(field,value);
  case {'stopGCV'}
    [validvalue, errmsg] = stopGCVtype(field,value);
  case {'stopCrit'}
    [validvalue, errmsg] = stopCritype(field,value);
  case {'stopOut'}
    [validvalue, errmsg] = stopOutype(field,value);
  case {'tollalpha'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'NE_Rtol'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'thr0'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'stabOut'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'RegMatrix'} % character string
    [validvalue, errmsg] = tikFunType(field,value);
  case {'enrichment'} % character string
    [validvalue, errmsg] = enrichType(field,value);
  case {'stdCGLS_out'}
    [validvalue, errmsg] = onOffType(field,value);
  case {'transf'} % character string
    [validvalue, errmsg] = spTransf(field,value);
  case {'Reorth'} % off, on
    [validvalue, errmsg] = onOffType(field,value);
  case {'stoprule'} % character string
    [validvalue, errmsg] = stopruleType(field,value);
  case {'xMin'} % real scalar
    [validvalue, errmsg] = realScalar(field,value);
  case {'xMax'} % real scalar
    [validvalue, errmsg] = realScalar(field,value);
  case {'xEnergy'} % real positive scalar or 'none'
    [validvalue, errmsg] = positiveScalar2(field,value);
  case {'taudelta'} % real non-negative scalar or 'none'
    [validvalue, errmsg] = nonNegScalar2(field,value);
  case {'nonneg'} % off, on
    [validvalue, errmsg] = onOffType(field,value);
  case {'verbosity'} % off, on
    [validvalue, errmsg] = onOffType(field,value);
  case {'Ubound'} % real positive scalar or 'none'
    [validvalue, errmsg] = positiveScalar2(field,value);
  case {'relaxParam'}% real non-negative scalar
    [validvalue, errmsg] = nonNegScalar2(field,value);
  case {'RegParam'} %
    [validvalue, errmsg] = RegPartype(field,value);
  case {'GCVweight'} % real non-negative scalar, adapt
    [validvalue, errmsg] = WeightType(field,value);
  case {'GCVflatTol'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'tolX'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'GCVminTol'}% real positive integer
    [validvalue, errmsg] = PosInteger(field,value);
  case {'resflatTol'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'jbdTol'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'jbdTolx'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'regPflatTol'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'LSQRtols'}% real non-negative vector of length 2
    [validvalue, errmsg] = nonNeg2vector(field,value);
  case {'RegParam0'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'RegParamVect0'}% real non-negative scalar
    [validvalue, errmsg] = nonNegVector(field,value);
  case {'NoiseLevel'}% real non-negative scalar
    [validvalue, errmsg] = NoiseLevelType(field,value);
  case {'eta'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'BestEnrm'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'trunc'}
    [validvalue, errmsg] = truncation(field,value);
  case {'noiseType'}
    [validvalue, errmsg] = noiseType(field,value);
  case {'noiseParam'}
    [validvalue, errmsg] = PosVector2(field,value);
  case {'inSolver'} % character string
    [validvalue, errmsg] = inSolverType(field,value);
  case {'solver'} % character string
    [validvalue, errmsg] = SolverType(field,value);
  case {'adaptConstr'} % character string
    [validvalue, errmsg] = adaptConstrType(field,value);
  case {'Ktot'} % vector with all real positive integers
    [validvalue, errmsg] = PosIntegerVector(field,value);
  case {'nonnegativity'} % off, on
    [validvalue, errmsg] = onOffType(field,value);
  case {'DecompOut'} % off, on
    [validvalue, errmsg] = onOffType(field,value);
  case{'sirt_method'}
    [validvalue, errmsg] = sirtType(field,value);
  case {'shrink'} % off, on
    [validvalue, errmsg] = onOffType(field,value);
  case {'stepsize'}% real non-negative scalar
    [validvalue, errmsg] = stepsizeoption(field,value);
  case {'backtracking'} % off, on
    [validvalue, errmsg] = onOffType(field,value);
  case {'backscalar'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'backit'} % real positive integer
    [validvalue, errmsg] = PosInteger(field,value);
  case{'SparsityTrans'}
    [validvalue, errmsg] = SparsityTransType(field,value);
  case{'hybridvariant'}
    [validvalue, errmsg] = hybridvariantType(field,value);
  case{'wname'}
    [validvalue, errmsg] = wnameType(field,value);
  case{'wlevels'}
    [validvalue, errmsg] = PosInteger(field,value);
  case{'warmrestart'} % off, on
    [validvalue, errmsg] = onOffType(field,value);
  case {'weight0'} % numeric array, none
    [validvalue, errmsg] = x0type(field,value);
  case {'svdbasis'} % numeric array, none
    [validvalue, errmsg] = PosIntegerBound(field,value);
  case {'qnorm'} % real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'dimension'} % real non-negative scalar
    [validvalue, errmsg] = dimensionsize(field,value);
  case {'reginskaExp'}
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'plotty'}
    [validvalue, errmsg] = onOffType(field,value);
  case {'RegParamRange'}
    [validvalue, errmsg] = nonNegBounds(field,value);
  case {'RegParRegRange'}
    [validvalue, errmsg] = nonNegBounds(field,value);
  case {'discrbilStopTol'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'discrflatTol'}
    [validvalue, errmsg] = nonNegscalar(field,value);
  case {'regbilStopTol'}% real non-negative scalar
    [validvalue, errmsg] = nonNegscalar(field,value);
  otherwise
    %validfield = false;
    validvalue = false;
    errmsg = sprintf('Unrecognized parameter name ''%s''.', field);
end

if validvalue
    return;
else
  error(errmsg)
end

%------------------------------------------------------------------------

function [valid, errmsg] = PosIntegerVector(field,value)
% Any vector with all postive real integers
valid =  (isreal(value) && (all(value > 0)) && all(value == floor(value))) | (ischar(value) && any(strcmp(value,{'none'}))) ;
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be positive real integers.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = PosVector2(field,value)
% Any vector with all postive real integers
valid =  (isreal(value) && (all(value > 0))) && (length(value) == 2) || (ischar(value) && any(strcmp(value,{'none'}))) ;
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be positive real vector of length 2.',field);
else
    if value(1)>value(2)
        warning('For the methods to be meaningful, the Poisson parameter (%s(2)) should be bigger than the Gaussian std deviation (%s(1)).', field, field);
        errmsg = '';
    else
        errmsg = '';
    end
end

%------------------------------------------------------------------------

function [valid, errmsg] = PosInteger(field,value)
% Any positive real integer
valid =  isreal(value) && isscalar(value) && (value > 0) && value == floor(value) ;
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a positive real integer.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = PosIntegerBound(field,value)
% Any positive real integer
valid =  isreal(value) && isscalar(value) && (value > 0) && (value < 7) && value == floor(value) ;
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a positive real integer between 1 and 6.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = PosInteger2(field,value)
% Any positive real integer
valid =  (isreal(value) && isscalar(value) && (value > 0) && value == floor(value)) || ...
    (ischar(value) && strcmp(value,{'none'}));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a positive real integer, or ''none''.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = nonNegInteger2(field,value)
% Any positive real integer
valid =  (isreal(value) && isscalar(value) && (value >= 0) && value == floor(value)) || ...
    (ischar(value) && strcmp(value,{'none'}));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a positive real integer, or ''none''.',field);
else
  errmsg = '';
end

%-----------------------------------------------------------------------

function [valid, errmsg] = onOffType(field,value)
% One of these strings: on, off
valid =  ischar(value) && any(strcmp(value,{'on';'off'}));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be ''off'' or ''on''.',field);
else
  errmsg = '';
end

%-----------------------------------------------------------------------

function [valid, errmsg] = SparsityTransType(field,value)
valid =  ischar(value) && any(strcmpi(value,{'none';'dwt';'dct';'svd';'tv1D';'tv2D'}));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be ''none'' or ''dwt'' or ''dct'' or ''svd'' or ''tv1D'' or ''tv2D''.',field);
else
  errmsg = '';
end

%-----------------------------------------------------------------------

function [valid, errmsg] = x_truetype(field,value)
% Either a numeric array or off
valid =  isnumeric(value) | (ischar(value) && any(strcmp(value,{'none'})));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a numeric array or ''none''.',field);
else
  errmsg = '';
end

%-----------------------------------------------------------------------

function [valid, errmsg] = x0type(field,value)
% Either a numeric array or none
valid =  isnumeric(value) | (ischar(value) && any(strcmp(value,{'none'})));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a numeric array or ''none''.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = nonNegscalar(field,value)
% Any real non-negative scalar
valid =  isreal(value) && isscalar(value) && (value >= 0);
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a real non-negative scalar.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = realScalar(field,value)
% Any real  scalar
valid =  (isreal(value) && isscalar(value)) || (ischar(value) && strcmpi(value, 'none'));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a real non-negative scalar.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = tikFunType(field,value)
% One of these strings: none, Identity, Laplacian
valid =  (ischar(value) && any(strcmpi(value,{'none';'Identity';'Laplacian1D';'Laplacian2D';'Gradient1D';'Gradient2D';'tv';'tv1D';'tv2D'}))) || ...
         (~ischar(value) && ismatrix(value));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be ''none'' or ''Identity'' or ''Laplacian1D'' or ''Laplacian2D'' or ''tv'' or a matrix.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = enrichType(field,value)
% One of these strings: none, Identity, Laplacian
valid =  (ischar(value) && (strcmpi(value,'none')||strcmpi(value,'ones')))|| ...
         (~ischar(value) && ismatrix(value));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be ''none'' or ''ones'' or a matrix of coherent size.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = spTransf(field,value)
% One of these strings: none, Identity, Laplacian
valid =  ischar(value) && any(strcmpi(value,{'none';'dct'}));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be ''none'' or ''dct''.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = stopruleType(field,value)
% One of these strings: none, Identity, Laplacian
valid =  ischar(value) && any(strcmpi(value,{'none';'NCP';'DP'}));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be ''none'' or ''NCP'' or ''DP''.',field);
else
  errmsg = '';
end

%-----------------------------------------------------------------------

function [valid, errmsg] = nonNegScalar2(field,value)
% Either a numeric array or none
valid =  (isreal(value) && isscalar(value) && value >= 0) || (ischar(value) && any(strcmp(value,{'none'})));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be non-negative scalar or ''none''.',field);
else
  errmsg = '';
end

%-----------------------------------------------------------------------

function [valid, errmsg] = positiveScalar2(field,value)
% Either a numeric array or none
valid =  (isreal(value) && isscalar(value) && value > 0) || (ischar(value) && any(strcmp(value,{'none'})));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be positive scalar or ''none''.',field);
else
  errmsg = '';
end%--------------------------------------------------------------------------

function [valid, errmsg] = RegPartype(field,value)
% One of these: real nonnegative scalar, GCV, WGCV, optimal
valid =  (isreal(value) && isscalar(value) && (value >= 0)) | (ischar(value) && any(strcmpi(value,{'gcv','wgcv','modgcv','optimal','discrep','discrepit','discrepbil','Lcurve','off','none','reginskait','reginskabil'})));

if ~valid
 errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a non-negative scalar or ''GCV'' or ''WGCV'' or ''discrep'' or ''discrepit'' or ''optimal'' or ''none''.',field);
else
  errmsg = '';
end

%--------------------------------------------------------------------------

function [valid, errmsg] = stopGCVtype(field,value)
% One of these: resflat, GCVflat
valid =  (ischar(value) && any(strcmpi(value,{'resflat', 'GCVvalues'})));
if ~valid
 errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be ''resflat'' or ''GCVvalues''.',field);
else
  errmsg = '';
end

%--------------------------------------------------------------------------

function [valid, errmsg] = stopCritype(field,value)
% One of these: resflat, GCVflat
valid =  (ischar(value) && any(strcmpi(value,{'none', 'resflat', 'discrep'}))); % 'outer'
if ~valid
 errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be ''resflat'' or ''discrep''.',field);
else
  errmsg = '';
end

%--------------------------------------------------------------------------

function [valid, errmsg] = stopOutype(field,value)
% One of these: xstab, Lxstab, regPstab
valid =  (ischar(value) && any(strcmpi(value,{'xstab', 'Lxstab', 'regPstab'})));
if ~valid
 errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be ''xstab'' or ''Lxstab'' or ''regPstab''.',field);
else
  errmsg = '';
end

%--------------------------------------------------------------------------

function [valid, errmsg] = WeightType(field,value)
% One of these: real non-negative scalar, adapt
valid =  (isreal(value) && isscalar(value) && (value >= 0) )| (ischar(value) && any(strcmp(value,{'adapt'})));
if ~valid
 errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a non-negative scalar or ''adapt''.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

%function [valid, errmsg] = InSolvetype(field,value)
% % One of these strings: tsvd, tikhonov, none
%valid =  ischar(value) && any(strcmp(value,{'tsvd','tikhonov','none'}));
%if ~valid
% errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be ''tsvd'', ''tikhonov'' or ''none''.',field);
%else
%  errmsg = '';
%end

%------------------------------------------------------------------------

function [valid, errmsg] = nonNeg2vector(field,value)
% Any real non-negative scalar
valid =  isreal(value) && isvector(value) && sum(value>=0)==length(value) && (length(value)==2);
if ~valid
 errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a real non-negative vector of length 2.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = nonNegBounds(field,value)
% Any real non-negative scalar
valid =  isreal(value) && isvector(value) && sum(value>=0)==length(value) && (length(value)==2) && (value(1)<value(2));
if ~valid
 errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a real non-negative vector of length 2, with the first entry smaller than the second entry.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = nonNegVector(field,value)
% Any real non-negative scalar
valid =  isreal(value) && isvector(value) && sum(value>=0)==length(value);
if ~valid
 errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a real non-negative vector of length 2.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = NoiseLevelType(field,value)
% One of these: real nonnegative scalar, GCV, WGCV, optimal
valid =  (isreal(value) && isscalar(value) && (value >= 0)) |...
         (ischar(value));

if ~valid
 errmsg = sprintf('Invalid value for OPTIONS parameter %s: a nonnegative scalar noise level should be provided.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = truncation(field,value)
% One of these: real nonnegative scalar, GCV, WGCV, optimal
valid =  (isreal(value) && isscalar(value) && (value > 0) && value == floor(value)) | (ischar(value) && any(strcmpi(value,{'full', 'none'})));
if ~valid
 errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a positive integer or ''none'' or ''full''.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = noiseType(field,value)
% One of these strings: none, Identity, Laplacian
valid =  ischar(value) && any(strcmpi(value,{'gauss';'poisson'}));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be ''gauss'' or ''poisson''.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = inSolverType(field,value)
% One of these strings: gmres, lsqr, fgmres
valid =  (ischar(value) && any(strcmpi(value,{'gmres';'fgmres';'lsqr';'cgls';'rrgmres'})));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS inSolver %s: must be ''gmres'' or ''fgmres'' or ''lsqr'' or ''cgls''.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = SolverType(field,value)
% One of these strings: gmres, lsqr, fgmres
valid =  (ischar(value) && any(strcmpi(value,{'fgmres';'flsqr'})));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS inSolver %s: must be ''fgmres'' or ''flsqr''.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = adaptConstrType(field,value)
% One of these strings: tv, nn
valid =  (ischar(value) && any(strcmpi(value,{'tv';'nn';'tvnn';'box';'energy';'project';'spnn';'sp';'none'})));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS inSolver %s: must be ''tv'' or ''nn'' or ''tvnn'' or ''box'' or ''spnn''.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = sirtType(field,value)
% One of these strings: cav, cimmino, drop, landweber, sart
valid =  (ischar(value) && any(strcmpi(value,{'cav';'cimmino';'drop';'landweber';'sart'})));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be ''cav'' or ''cimmino'' or ''drop'' or ''landweber'' or ''sart''.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = hybridvariantType(field,value)
% One of these strings: tv, nn
valid =  (ischar(value) && any(strcmpi(value,{'I';'R'})));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS inSolver %s: must be ''I'' or ''R''.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = wnameType(field,value)
% One of these strings: cav, cimmino, drop, landweber, sart
valid =  (ischar(value) && any(strcmpi(value,{'db1'})));
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be ''db1''.',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = dimensionsize(field,value)
% Any positive real integer
valid =  (value == 1 || value == 2);
if ~valid
  errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be either 1 or 2 (scalar).',field);
else
  errmsg = '';
end

%------------------------------------------------------------------------

function [valid, errmsg] = stepsizeoption(field,value)
% Any real non-negative scalar
valid =  (isreal(value) && isscalar(value) && (value >= 0)) || (ischar(value) && (strcmp(value, 'none') ));

if ~valid
 errmsg = sprintf('Invalid value for OPTIONS parameter %s: must be a real non-negative scalar.',field);
else
  errmsg = '';
end

repo: IRtools


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

%%
clear; close all; clc;

%% Syntax
% |cMap=gjet(varargin);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function varargout = struct2double(varargin)

% This function is a backward compatibility wrapper. It allows existing
% MATLAB scripts that do not use the new FieldTrip ft_xxx function naming
% scheme to work with recent versions of the FieldTrip toolbox.
%
% Please look in ft_xxx for the help of the function that you are looking
% for, where xxx is the name of the function that you were looking for.

% Copyright (C) 2009, Robert Oostenveld

prefix    = 'ft_';
funname   = mfilename;
warning([upper(mfilename) ' is only a compatibility wrapper, which will soon be removed. Please instead call ' upper(prefix) upper(funname) '.']);
funhandle = str2func([prefix funname]);
[varargout{1:nargout}] = funhandle(varargin{:});

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% GAMLPR = Gamma Distribution - Log Probability Density Ratio
% Copyright (c) 1998, Harvard University. Full copyright in the file Copyright
%
%   [ lpr ] = gamlpr( g1, g2, alph, gam )
%
% returns log ( p(g1)/p(g2) ) when g1 and g2 are
% distributed gamma(alph,gam)
%

function [ lpr ] = gamlpr(g1, g2, alph, gam)
lpr = (alph-1)*log(g1/g2)  - (g1-g2)/gam ;

author: Sable

repo: mcbench-benchmarks


text:
function ar=lpccc2ar(cc)
%LPCCC2AR Convert complex cepstrum to ar coefficients AR=(CC)
%
% MATLAB5 version

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: lpccc2ar.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p]=size(cc);
rp=-(1:p);
cc = cc .* rp(ones(nf,1),:);
if p<2
  ar = [ones(nf,1) cc(:,1)];
else
 ar=zeros(nf,p+1);
 ar(:,1:3) = [ones(nf,1) cc(:,1) (cc(:,2)+cc(:,1).^2)/2];
  for k=3:p
    ar(:,k+1) = (cc(:,k) + sum(cc(:,1:k-1).*ar(:,k:-1:2),2))/k;
  end
end

author: decouples


text:
% A LP example which shows all the potentials of GLPKMEX
clear;

disp('LP problem');
s=-1;
c=[10,6,4]';
a=[1,1,1;...
   10,4,5;...
   2,2,6];
b=[100,600,300]';
ctype=['U','U','U']';
lb=[0,0,0]';
ub=[]';
vartype=['C','C','C']';
% Output all GLPK messages on workspace
param.msglev=3;
% Set save options
param.save=1;
param.savefilename='SimpleLP';
param.savefiletype='fixedmps';
[xmin,fmin,status,extra]=glpk(c,a,b,lb,ub,ctype,vartype,s,param)

% OBSOLETE SYNTAX
%lpsolver = param.lpsolver;
%save_pb = param.save;
%[xmin,fmin,status,extra]=glpkmex(s,c,a,b,ctype,lb,ub,vartype,param,lpsolver,save_pb)

author: opencobra

repo: cobratoolbox


text:
% MEEG Object
%
% Functionalities of the MEEG object are described in the SPM manual,
% Chapter 13, Section 7.
% _________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: Contents.m 2696 2009-02-05 20:29:48Z guillaume $


text:
function pass = test_jacobian( pref )
% Test Jacobian

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 100*pref.cheb2Prefs.chebfun2eps;

% Check definition:
F = chebfun2v(@(x,y) cos(x), @(x,y) sin(y));
Fx = diffx(F);
Fy = diffy(F);
jacF = Fx(1).*Fy(2) - Fx(2).*Fy(1);

pass(1) = ( norm(jacF - jacobian(F) ) < tol );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function  write2hdf5
% write image data (and simple feature data) from one sequence to hdf5 file
% positive and negative example seperate
% IMPORTANT!!!
% All the data and labels saved into hdf5 should be sing
data_dir = '/scratch/BS/pool2/tang/multicut/pairwise_example/';
data_set = 'TownCenter-Sample';%'tud-crossing';%'ETH-Bahnhof';%'ETH-Pedcross2';%

chunksz = 1000;

positive_img_dir  = [data_dir data_set '/positive/'];
negative_img_dir  = [data_dir data_set '/negative/'];
save_dir = [data_dir data_set '/mixed_hdf5/'];
if ~exist(save_dir)
    mkdir(save_dir);
end

fprintf('Listing all the images ...\n')
positive_img_list = dir([positive_img_dir '*.jpg']);
negative_img_list = dir([negative_img_dir '*.jpg']);

num_img = length(positive_img_list) + length(negative_img_list);
num_round = (num_img/2)/chunksz;
parfor k = 1:num_round
    fprintf('Saving hdf5 file %d ...\n', k);
    last_end_img_idx = (k-1)*chunksz;

    batchdata = single(zeros(227,227,3,chunksz*2));
    batchlabs = single(zeros(1,chunksz*2));

    filename = [save_dir sprintf('%05d',last_end_img_idx + 1) '.h5'];
    count = 0;
    for i = 1:chunksz
        idx = last_end_img_idx + i;
        positive_img_name  = positive_img_list(idx);
        img = imread([ positive_img_dir positive_img_name.name]);
        % unit to single
        img = imresize(img,[227,227]);

        img =  single(img);
%         img(:,120:end,1) = 255;
        %
        img = img(:,:,[3 2 1]);
        %
        img = permute(img, [2 1 3]);
        %

        count = count +1;
        batchdata(:,:,:,count) = img;
        batchlabs(count) = 1;

        negative_img_name  = negative_img_list(idx);
        img = imread([ negative_img_dir negative_img_name.name]);
        % unit to single
        img = imresize(img,[227,227]);
        img =  single(img);
        %
        img = img(:,:,[3 2 1]);
        %
        img = permute(img, [2 1 3]);

        count = count +1;
        batchdata(:,:,:,count) = img;
        batchlabs(count) = 0;
    end
    curr_dat_sz=store2hdf5(filename, batchdata, batchlabs, 1);
end
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function prob_test064 ( )

%*****************************************************************************80
%
%% TEST064 tests EXTREME_VALUES_CDF, EXTREME_VALUES_CDF_INV, EXTREME_VALUES_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST064\n' );
  fprintf ( 1, '  For the Extreme Values CDF:\n' );
  fprintf ( 1, '  EXTREME_VALUES_CDF evaluates the CDF;\n' );
  fprintf ( 1, '  EXTREME_VALUES_CDF_INV inverts the CDF.\n' );
  fprintf ( 1, '  EXTREME_VALUES_PDF evaluates the PDF;\n' );

  a = 2.0;
  b = 3.0;

  check = extreme_values_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST064 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A = %14f\n', a );
  fprintf ( 1, '  PDF parameter B = %14f\n', b );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '       X            PDF           CDF            CDF_INV\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ x, seed ] = extreme_values_sample ( a, b, seed );

    pdf = extreme_values_pdf ( x, a, b );

    cdf = extreme_values_cdf ( x, a, b );

    x2 = extreme_values_cdf_inv ( cdf, a, b );

    fprintf ( 1, ' %14f  %14f  %14f  %14f\n', x, pdf, cdf, x2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function h = togglePlot(Y, E, X, names, colors, styles);
% Plot time curves with errorbars, and toggles for individual curves.
%
% h = togglePlot(Y, [E], [X], [names], [colors], [styles]);
%
% This is an updated version of my old function FANCYPLOT. It produces
% a plot of the data Y, at points X, with error bar size E. For X, Y, and
% E, different columns represent different errorbar lines while rows
% represent different points along X. If X is omitted, it initializes to
% [1:size(Y, 1)] for each column. A legend is appended, along with
% checkboxes for each column in Y, which toggles the visibility of that
% data line.
%
% names, colors, and styles are all optional cell arrays specifying the
% names of each column in Y for the legend, the colors associated with each
% line, and the line styles (as set per SETLINESTYLES).
%
% Returns an aray of handles to the plot elements, in the format
%   [figure, axes, legendPanel, (errorbar handles) (checkbox handles)].
%
% EXAMPLE:
%   X = linspace(0, 2*pi, 40);
%   Y = [cos(X); sin(X); cos(X).^2; tan(X)]';
%   E = .1 * rand(40, 4);
%   names = {'cos' 'sin' 'cos^2' 'tan'};
%   colors = {[1 0 0] [.9 .2 0] [.7 .3 .1] [.5 .5 .3]};
%   styles = {'-2' '--2' '.-2' ':2'};
%   h = togglePlot(Y, E, X, names, colors, styles);
%
% ras, 03/2007.
if notDefined('Y'),     error('Need to specify data Y.');       end
if notDefined('X'),     X = repmat([1:size(Y,1)]', [1 size(Y,2)]);  end
if notDefined('E'),     E = zeros(size(Y));                     end
if notDefined('styles'),    styles = {};                    end

nCols = size(Y, 2);

if notDefined('names'),
    for i = 1:nCols
        names{i} = sprintf('Col %i', i);
    end
end

% size check on X
if ~isequal(size(X), size(Y))
    if numel(X)==size(Y, 1)
        % single X vec for all columns
        X = repmat(X(:), [1 nCols]);

    else
        error('X and Y should be same size.')

    end
end

%% create figure, axes
h(1) = figure('Color', 'w');

h(2) = axes;

%% plot the data
hBars = errorbar(X, Y, E);

% get default colors if not specified
if notDefined('colors'),
    colorOrder = get(h(2), 'ColorOrder');
    while size(colorOrder, 1) < nCols
        colorOrder = [colorOrder; colorOrder];
    end

    for col = 1:nCols
        colors{col} = colorOrder(col,:);
    end
end

setLineColors(colors);

if ~isempty(styles)
    setLineStyles(styles);
end

%% add legend panel, checkboxes
hLeg = legendPanel(names, colors);

% callback for each checkbox
cb = ['if get(gcbo,''Value'')==1, TMP=''on''; else, TMP=''off''; end; ' ...
      'set(get(gcbo,''UserData''), ''Visible'', TMP); ' ...
      'clear TMP; '];

% add checkboxes
for n = 1:nCols
    % these are positions of the axes in LEGENDPANEL
    row = mod(n-1, 20) + 1;
    pos = [.4,  .96-row*.05,  .1,  .025];

    hCheck(n) = uicontrol('Parent', hLeg, 'Style', 'checkbox', ...
                    'Units', 'normalized', 'Position', pos, ...
                    'UserData', hBars(n), 'Callback', cb, ...
                    'String', '', 'BackgroundColor', 'w', 'Value', 1);
end

%% append all handles to h
h = [h hLeg hBars hCheck];

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function exact = p01_exact ( )

%*****************************************************************************80
%
%% P01_EXACT returns the exact integral for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  omega = 1.0;

  exact = sqrt ( pi ) * exp ( - omega * omega );

  return
end


text:
classdef math_model_opf_accs < mp.math_model_opf_acc
%MP.MATH_MODEL_OPF_ACCS  MATPOWER mathematical model for AC optimal power flow (OPF) problem.
%   ?
%
%   MP.MATH_MODEL_OPF_ACCS ... power flow ...
%
%   Properties
%       ? - ?
%
%   Methods
%       ?

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

%     properties
%     end

    methods
        function tag = form_tag(obj)
            tag = 'accs';
        end

        function name = form_name(obj)
            name = 'AC-cartesian-power';
        end

        function add_node_balance_constraints(obj, nm, dm, mpopt)
            %% power balance constraints
            nn = nm.node.N;             %% number of nodes
            fcn_mis = @(x)obj.nodal_power_balance_fcn(x, nm);
            hess_mis = @(x, lam)obj.nodal_power_balance_hess(x, lam, nm);
            obj.add_nln_constraint({'Pmis', 'Qmis'}, [nn;nn], 1, fcn_mis, hess_mis);
        end

        function [lam_p, lam_q] = node_power_balance_prices(obj, nm)
            %% shadow prices on node power balance
            nne = obj.get_idx('nle');
            lambda = obj.soln.lambda;
            lam_p = lambda.eqnonlin(nne.i1.Pmis:nne.iN.Pmis);
            lam_q = lambda.eqnonlin(nne.i1.Qmis:nne.iN.Qmis);
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
clc;clear all;close all;
addpath ../test
addpath ..
addpath CM_Curve
%netStruct = load('../data/res52_drop0.75_batch16_psudo_label/net-epoch-25.mat');
%netStruct.net.layers(end-3)=[];
%netStruct.net.vars(178)=[];
%net = dagnn.DagNN.loadobj(netStruct.net);
%net.addLayer('feature',dagnn.Concat('dim',3),{'pool5','pool5_local'},{'pool5_fine'});
%net.mode = 'test' ;
%net.move('gpu') ;
rank_size = 2000;
%net.conserveMemory = false;
%im_mean = net.meta.normalization.averageImage;
%im_mean = imresize(im_mean,[224,224]);

%% add necessary paths
query_dir = '/DATACENTER/1/qxl/PRID/Market-1501/query/';% query directory
test_dir = '/DATACENTER/1/qxl/PRID/Market-1501/bounding_box_test/';% database directory
feature_dir = '/home/qxl/Desktop/';

%% calculate query features
Hist_query = importdata([feature_dir 'Hist_query.mat']);
nQuery = size(Hist_query, 2);
%nQuery = 100;
%% calculate database features
Hist_test = importdata([feature_dir 'Hist_test.mat']);
nTest = size(Hist_test, 2);
%nTest = 19732;
%% calculate the ID and camera for database images
mkdir('./data')
test_files = dir([test_dir '*.jpg']);
testID = zeros(length(test_files), 1);
testCAM = zeros(length(test_files), 1);
if ~exist('data/testID.mat')
    for n = 1:length(test_files)
        img_name = test_files(n).name;
        if strcmp(img_name(1), '-') % junk images
            testID(n) = -1;
            testCAM(n) = str2num(img_name(5));
        else
            %img_name
            testID(n) = str2num(img_name(1:4));
            testCAM(n) = str2num(img_name(7));
        end
    end
    save('data/testID.mat', 'testID');
    save('data/testCAM.mat', 'testCAM');
else
    testID = importdata('data/testID.mat');
    testCAM = importdata('data/testCAM.mat');
end

%% calculate the ID and camera for query images
query_files = dir([query_dir '*.jpg']);
queryID = zeros(length(query_files), 1);
queryCAM = zeros(length(query_files), 1);
if ~exist('data/queryID.mat')
    for n = 1:length(query_files)
        img_name = query_files(n).name;
        if strcmp(img_name(1), '-') % junk images
            queryID(n) = -1;
            queryCAM(n) = str2num(img_name(5));
        else
            queryID(n) = str2num(img_name(1:4));
            queryCAM(n) = str2num(img_name(7));
        end
    end
    save('data/queryID.mat', 'queryID');
    save('data/queryCAM.mat', 'queryCAM');
else
    queryID = importdata('data/queryID.mat');
    queryCAM = importdata('data/queryCAM.mat');
end

%% search the database and calcuate re-id accuracy
ap = zeros(nQuery, 1); % average precision
ap_max_rerank  = zeros(nQuery, 1); % average precision with MultiQ_max + re-ranking
ap_pairwise = zeros(nQuery, 6); % pairwise average precision with single query (see Fig. 7 in the paper)

CMC = zeros(nQuery, rank_size);
CMC_max_rerank = zeros(nQuery, rank_size);

r1 = 0; % rank 1 precision with single query
r1_max_rerank = 0; % rank 1 precision with MultiQ_max + re-ranking
r1_pairwise = zeros(nQuery, 6);% pairwise rank 1 precision with single query (see Fig. 7 in the paper)

%dist = importdata('/home/qxl/Desktop/score_distmat.mat');
dist = sqdist(Hist_test, Hist_query); % distance calculate with single query. Note that Euclidean distance is equivalent to cosine distance if vectors are l2-normalized
%dist_cos_max = (2-dist_max)./2; % cosine distance with MultiQ_max, used for re-ranking

knn = 1; % number of expanded queries. knn = 1 yields best result
%queryCam = importdata('data/queryCAM_duke.mat'); % camera ID for each query
%testCam = importdata('data/testCAM_duke.mat'); % camera ID for each database image

for k = 1:nQuery
    % load ground truth for each query (good and junk)
    good_index = intersect(find(testID == queryID(k)), find(testCAM ~= queryCAM(k)))';% images with the same ID but different camera from the query
    junk_index1 = find(testID == -1);% images neither good nor bad in terms of bbox quality
    junk_index2 = intersect(find(testID == queryID(k)), find(testCAM == queryCAM(k))); % images with the same ID and the same camera as the query
    junk_index = [junk_index1; junk_index2]';
    score = dist(:, k);
    %score_avg = dist_avg(:, k);
    %score_max = dist_max(:, k);

    % sort database images according Euclidean distance
    [~, index] = sort(score, 'ascend');  % single query

    % re-rank  select rank_size=1000 index
    index = index(1:rank_size);
    %for l=1:length(index)
    %    fprintf('%d, %s\n',index(l),test_files(index(l)).name)
    %end

    [ap(k), CMC(k, :)] = compute_AP_rerank(good_index, junk_index, index);% compute AP for single query
    %ap_pairwise(k, :) = compute_AP_multiCam(good_index, junk_index1,junk_index2, index, queryCAM(k), testCAM); % compute pairwise AP for single query
    fprintf('%d,%s::%f\n',k,query_files(k).name,ap(k));
    %%%%%%%%%%% calculate pairwise r1 precision %%%%%%%%%%%%%%%%%%%%
    %r1_pairwise(k, :) = compute_r1_multiCam(good_index, junk_index1,junk_index2, index, queryCAM(k), testCAM); % pairwise rank 1 precision with single query
    %%%%%%%%%%%%%% calculate r1 precision %%%%%%%%%%%%%%%%%%%%
end
CMC = mean(CMC);
%% print result
fprintf('single query:    mAP = %f, r1 precision = %f\r\n', mean(ap), CMC(1));
%[ap_CM, r1_CM] = draw_confusion_matrix(ap_pairwise, r1_pairwise, queryCam);
%fprintf('average of confusion matrix with single query:  mAP = %f, r1 precision = %f\r\n', (sum(ap_CM(:))-sum(diag(ap_CM)))/30, (sum(r1_CM(:))-sum(diag(r1_CM)))/30);

%% plot CMC curves
figure;
s = 50;
CMC_curve = CMC ;
plot(1:s, CMC_curve(:, 1:s));

author: naiq


text:
% This file is part of TREEQSM.
%
% TREEQSM is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% TREEQSM is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with TREEQSM.  If not, see <http://www.gnu.org/licenses/>.

function [Curve,Ind] = boundary_curve(P,Curve0,rball,dmax)

% ---------------------------------------------------------------------
% BOUNDARY_CURVE.M      Determines the boundary curve based on the
%                           previously defined boundary curve.
%
% Version 1.1.0
% Latest update     3 May 2022
%
% Copyright (C) 2015-2022 Pasi Raumonen
% ---------------------------------------------------------------------
%
% Inputs:
% P         Point cloud of the cross section
% Curve0     Seed points from previous cross section curve
% rball     Radius of the balls centered at seed points
% dmax      Maximum distance between concecutive curve points, if larger,
%               then create a new one between the points
% ---------------------------------------------------------------------

% Changes from version 1.0.0 to 1.1.0, 3 May 2022:
% 1) Increased the cubical neighborhood in the generation of the segments

%% Partition the point cloud into cubes
Min = double(min([P(:,1:2); Curve0(:,1:2)]));
Max = double(max([P(:,1:2); Curve0(:,1:2)]));
N = double(ceil((Max-Min)/rball)+5);
% cube coordinates of the section points
CC = floor([P(:,1)-Min(1) P(:,2)-Min(2)]/rball)+3;
% Sorts the points according a lexicographical order
S = [CC(:,1) CC(:,2)-1]*[1 N(1)]';
[S,I] = sort(S);
% Define "partition"
np = size(P,1);
partition = cell(N(1),N(2));
p = 1;              % The index of the point under comparison
while p <= np
  t = 1;
  while (p+t <= np) && (S(p) == S(p+t))
    t = t+1;
  end
  q = I(p);
  partition{CC(q,1),CC(q,2)} = I(p:p+t-1);
  p = p+t;
end

%% Define segments using the previous points
% cube coordinates of the seed points:
CC = floor([Curve0(:,1)-Min(1) Curve0(:,2)-Min(2)]/rball)+3;
I = CC < 3;
CC(I) = 3;
nc = size(Curve0,1);  % number of sets
Dist = 1e8*ones(np,1);  % distance of point to the closest center
SoP = zeros(np,1);  % the segment the points belong to
Radius = rball^2;
for i = 1:nc
  points = partition(CC(i,1)-2:CC(i,1)+2,CC(i,2)-2:CC(i,2)+2);
  points = vertcat(points{:});
  V = [P(points,1)-Curve0(i,1) P(points,2)-Curve0(i,2)];
  dist = sum(V.*V,2);
  PointsInBall = dist < Radius;
  points = points(PointsInBall);
  dist = dist(PointsInBall);
  D = Dist(points);
  L = dist < D;
  I = points(L);
  Dist(I) = dist(L);
  SoP(I) = i;
end

%% Finalise the segments
% Number of points in each segment and index of each point in its segment
Num = zeros(nc,1);
IndPoints = zeros(np,1);
for i = 1:np
  if SoP(i) > 0
    Num(SoP(i)) = Num(SoP(i))+1;
    IndPoints(i) = Num(SoP(i));
  end
end
% Continue if enough non-emtpy segments
if nnz(Num) > 0.05*nc
  % Initialization of the "Seg"
  Seg = cell(nc,1);
  for i = 1:nc
    Seg{i} = zeros(Num(i),1);
  end
  % Define the "Seg"
  for i = 1:np
    if SoP(i) > 0
      Seg{SoP(i),1}(IndPoints(i),1) = i;
    end
  end

  %% Define the new curve points as the average of the segments
  Curve = zeros(nc,3);  % the new boundary curve
  Empty = false(nc,1);
  for i = 1:nc
    S = Seg{i};
    if ~isempty(S)
      Curve(i,:) = mean(P(S,:),1);
      if norm(Curve(i,:)-Curve0(i,:)) > 1.25*dmax
        Curve(i,:) = Curve0(i,:);
      end
    else
      Empty(i) = true;
    end
  end

  %% Interpolate for empty segments
  % For empty segments create points by interpolation from neighboring
  % non-empty segments
  if any(Empty)
    for i = 1:nc
      if Empty(i)
        if i > 1 && i < nc
          k = 0;
          while i+k <= nc && Empty(i+k)
            k = k+1;
          end
          if i+k <= nc
            LineEle = Curve(i+k,:)-Curve(i-1,:);
          else
            LineEle = Curve(1,:)-Curve(i-1,:);
          end
          if k < 5
            for j = 1:k
              Curve(i+j-1,:) = Curve(i-1,:)+j/(k+1)*LineEle;
            end
          else
            Curve(i:i+k-1,:) = Curve0(i:i+k-1,:);
          end
        elseif i == 1
          a = 0;
          while Empty(end-a)
            a = a+1;
          end
          b = 1;
          while Empty(b)
            b = b+1;
          end
          LineEle = Curve(b,:)-Curve(nc-a,:);
          n = a+b-1;
          if n < 5
            for j = 1:a-1
              Curve(nc-a+1+j,:) = Curve(nc-a,:)+j/n*LineEle;
            end
            for j = 1:b-1
              Curve(j,:) = Curve(nc-a,:)+(j+a-1)/n*LineEle;
            end
          else
            Curve(nc-a+2:nc,1:2) = Curve0(nc-a+2:nc,1:2);
            Curve(nc-a+2:nc,3) = Curve0(nc-a+2:nc,3);
            Curve(1:b-1,1:2) = Curve0(1:b-1,1:2);
            Curve(1:b-1,3) = Curve0(1:b-1,3);
          end
        elseif i == nc
          LineEle = Curve(1,:)-Curve(nc-1,:);
          Curve(i,:) = Curve(nc-1,:)+0.5*LineEle;
        end
      end
    end
  end

  % Correct the height
  Curve(:,3) = min(Curve(:,3));

  % Check self-intersection
  [Intersect,IntersectLines] = check_self_intersection(Curve(:,1:2));

  % If self-intersection, try to modify the curve
  j = 1;
  while Intersect && j <= 5
    n = size(Curve,1);
    InterLines = (1:1:n)';
    NumberOfIntersections = cellfun('length',IntersectLines(:,1));
    I = NumberOfIntersections > 0;
    InterLines = InterLines(I);
    CrossLen = vertcat(IntersectLines{I,2});
    if length(CrossLen) == length(InterLines)
      LineEle = Curve([2:end 1],:)-Curve(1:end,:);
      d = sqrt(sum(LineEle.*LineEle,2));
      m = length(InterLines);
      for i = 1:2:m
        if InterLines(i) ~= n
          Curve(InterLines(i)+1,:) = Curve(InterLines(i),:)+...
            0.9*CrossLen(i)/d(InterLines(i))*LineEle(InterLines(i),:);
        else
          Curve(1,:) = Curve(InterLines(i),:)+...
            0.9*CrossLen(i)/d(InterLines(i))*LineEle(InterLines(i),:);
        end
      end
      [Intersect,IntersectLines] = check_self_intersection(Curve(:,1:2));
      j = j+1;
    else
      j = 6;
    end
  end

  %% Add new points if too large distances
  LineEle = Curve([2:end 1],:)-Curve(1:end,:);
  d = sum(LineEle.*LineEle,2);
  Large = d > dmax^2;
  m = nnz(Large);
  if m > 0
    Curve0 = zeros(nc+m,3);
    Ind = zeros(nc+m,2);
    t = 0;
    for i = 1:nc
      if Large(i)
        t = t+1;
        Curve0(t,:) = Curve(i,:);
        if i < nc
          Ind(t,:) = [i i+1];
        else
          Ind(t,:) = [i 1];
        end
        t = t+1;
        Curve0(t,:) = Curve(i,:)+0.5*LineEle(i,:);
        if i < nc
          Ind(t,:) = [i+1 0];
        else
          Ind(t,:) = [1 0];
        end
      else
        t = t+1;
        Curve0(t,:) = Curve(i,:);
        if i < nc
          Ind(t,:) = [i i+1];
        else
          Ind(t,:) = [i 1];
        end
      end
    end
    Curve = Curve0;

  else
    Ind = [(1:1:nc)' [(2:1:nc)'; 1]];
  end

  %% Remove new points if too small distances
  nc = size(Curve,1);
  LineEle = Curve([2:end 1],:)-Curve(1:end,:);
  d = sum(LineEle.*LineEle,2);
  Small = d < (0.333*dmax)^2;
  m = nnz(Small);
  if m > 0
    for i = 1:nc-1
      if ~Small(i) && Small(i+1)
        Ind(i,2) = -1;
      elseif Small(i) && Small(i+1)
        Small(i+1) = false;
      end
    end
    if ~Small(nc) && Small(1)
      Ind(nc,2) = -1;
      Ind(1,2) = -1;
      Small(1) = false;
      Small(nc) = true;
      I = Ind(:,2) > 0;
      Ind(2:end,1) = Ind(2:end,1)+1;
      Ind(I,2) = Ind(I,2)+1;

    end
    Ind = Ind(~Small,:);
    Curve = Curve(~Small,:);
  end

else
  % If not enough new points, return the old curve
  Ind = [(1:1:nc)' [(2:1:nc)'; 1]];
  Curve = Curve0;
end

author: InverseTampere

repo: TreeQSM


text:
function h_cb= plotutil_colorbarAside(varargin)
%PLOT_COLORBARASIDE - Add a colorbar without shrinking the axis
%
%Synopsis:
% H_CB= plotutil_colorbarAside(<OPT>)
% H_CB= plotutil_colorbarAside(ORIENTATION, <OPT>)
%
%Input:
% ORIENTATION: see OPT.Orientation
% OPT: struct or property/value list of optional properties
%  .orientation - orientation resp. location of the colorbar relative
%                 to the axis: 'vert','horiz', 'NorthOutside',
%                 'EastOutside','SouthOutside', 'WestOutside'
%
%Output:
% HCB: handle of the colorbar
%
%Note:
% So far, plotutil_colorbarAside is compatible with Matlab 6, but it allows
% the new orientation modes of Matlab 7.

% blanker@cs.tu-berlin.de

props = {'Orientation'      'vert'      '!CHAR(vert horiz NorthOutside EastOutside SouthOutside WestOutside)';
         'Gap'              0.02        '!DOUBLE';
         };

if mod(nargin, 2)==1,
  varargin = {'Orientation'  varargin{:}};
end

opt= opt_proplistToStruct(varargin{:});
[opt,isdefault] = opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

ax= gca;
pos= get(ax, 'position');
% this has an influence on the behavious when resizing the figure:
epos= axis_getEffectivePosition(ax);
% this would require Matlab 7:
% h_cb= colorbar(opt.Orientation);
if ismember(lower(opt.Orientation), {'horiz','northoutside','southoutside'},'legacy'),
  h_cb= colorbar('horiz');
%  h_cb= colorbar('SouthOutside');
else
  h_cb= colorbar('EastOutside');
end
set(ax, 'position',pos);
drawnow;
cb_pos= get(h_cb, 'position');
%ii= strmatch(opt.Orientation, {'vert','horiz'}, 'exact');
%cb_pos(ii)= pos(ii)+pos(ii+2)+0.02;
%cb_pos(5-ii)= min(cb_pos(5-ii), epos(5-ii));
%cb_pos(3-ii)= pos(3-ii) + (pos(5-ii)-cb_pos(5-ii))/2;
switch(lower(opt.Orientation)),
 case {'vert','eastoutside'},
  cb_pos(1)= pos(1) + pos(3) + opt.Gap;
  cb_pos(4)= min(cb_pos(4), epos(4));
  cb_pos(2)= pos(2) + (pos(4)-cb_pos(4))/2;
 case 'westoutside',
  cb_pos(1)= pos(1) - cb_pos(3) - opt.Gap;
  cb_pos(4)= min(cb_pos(4), epos(4));
  cb_pos(2)= pos(2) + (pos(4)-cb_pos(4))/2;
  set(h_cb, 'YAxisLocation','left');
 case {'horiz','southoutside'},
  cb_pos(2)= pos(2) - cb_pos(4) - opt.Gap;
  cb_pos(3)= min(cb_pos(3), epos(3));
  cb_pos(1)= pos(1) + (pos(3)-cb_pos(3))/2;
 case 'northoutside',
  cb_pos(2)= pos(2) + pos(4) + opt.Gap;
  cb_pos(3)= min(cb_pos(3), epos(3));
  cb_pos(1)= pos(1) + (pos(3)-cb_pos(3))/2;
  set(h_cb, 'XAxisLocation','top');
 otherwise,
  error('unknown orientation/location');
end
set(h_cb, 'Position',cb_pos);
%delete(h_cb);
%cb_ax= axes('position',cb_pos);
%h_cb= colorbar(cb_ax, 'peer',ax);

if nargout==0,
  clear h_cb;
end

author: bbci


text:
 function [sig, kernel, kernel_ft] = nufft_best_gauss(J, K_N, sn_type)
%function [sig, kernel, kernel_ft] = nufft_best_gauss(J, K_N, sn_type)
%|
%| Return "sigma" of best (truncated) gaussian for NUFFT
%| with previously numerically-optimized width
%|
%| in
%|	K_N		K/N
%|	J		# of neighbors used per frequency location
%|	sn_type		'zn' or 'ft' (latter recommended)
%|
%| out
%|	sig		best sigma
%|	kernel		string for inline kernel function, args (k,J)
%|	kernel_ft	string for Fourier transform function, arg: (t)
%|
%| Copyright 2002-4-11, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end

if nargin == 1 && streq(J, 'test')
	[sig, kernel, kernel_ft] = nufft_best_gauss(6)
	clear sig kernel kernel_ft
return
end

if nargin < 2, K_N = 2; end
if nargin < 3, sn_type = 'ft'; end

if K_N ~= 2, error 'only K/N=2 done', end

%if ir_is_octave % avoid annoying load warning
%	warning('off', 'Octave:load-file-in-path', 'local')
%end

nufft_dir = path_find_dir('nufft');
s = load([nufft_dir filesep 'param-data/nufft_gauss2.mat']);

ii = find(J == s.Jgauss2);
if length(ii) ~= 1
	disp(s.Jgauss2(:)')
	error 'only above J values done'
end
if streq(sn_type, 'ft')
	sig = s.Sgauss2.ft(find(J == s.Jgauss2));
elseif streq(sn_type, 'zn')
	sig = s.Sgauss2.zn(find(J == s.Jgauss2));
else
	error 'bad sn_type'
end
[kernel, kernel_ft] = nufft_gauss('string', J, sig);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function this = plot(this, module, varargin)
% Plots different aspects of MrSeries using MrImage-plot functions
%
%   Y = MrSeries()
%   Y.plot(inputs)
%
% This is a method of class MrSeries.
%
% IN
%
% OUT
%
% EXAMPLE
%   plot
%
%   See also MrSeries

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-07-08
% Copyright (C) 2014 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 Licence (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/>.

if ~exist('module', 'var')
    module = 'data';
end

switch module
    case {'data', 'mean', 'sd', 'snr', 'coeffVar', 'diffLastFirst'}
        this.(module).plot(varargin{:});
    otherwise
        error('tapas:uniqc:MrSeries:InvalidPlotModule', ...
            'plotting module %s not implemented for MrSeries', ...
            module);
end

author: translationalneuromodeling

repo: tapas


function [X,Y] = ToyCarPairsToLabels(pairs, instance_matrix)
        imgs1 = [pairs.img1];
        imgs2 = [pairs.img2];

        keymat = [[imgs1.classId] [imgs2.classId]; [imgs1.id] [imgs2.id]];
        ids = unique(keymat','rows');

        X = instance_matrix(:,ids(:,2));
        Y = ids(:,1)';
end


text:
function [NV,NF,OV,OE,OV1,OE1,epsilon] = remesh_at_handles(V,F,C,P,BE,CE)
  % REMESH_AT_HANDLES  remesh a given mesh using its outline and constraining
  % interrior controls at given handles. Point handles are met exactly and
  % bones are sample at the given rate.
  %
  % [NV,NF] = remesh_at_handles(V,F,C,P,BE,CE)
  %
  % Inputs:
  %  V  list of vertex positions
  %  F  list of face indices
  %  C  list of control vertex positions
  %  P  list of indices into C for point controls, { 1:size(C,1) }
  %  BE  list of bones, pairs of indices into C, connecting control vertices,
  %    { [] }
  %  CE  list of "cage edges", pairs of indices into ***P***, connecting
  %    control ***points***. A "cage edge" just tells point boundary conditions
  %    to vary linearly along straight lines between the end points, and to be
  %    zero for all other handles. { [] }
  % Outputs
  %  NV  new list of vertex positions
  %  NF  new list of face indices
  %  OV  list of vertex positions sent to triangle
  %  OE  list of edges sent to triangle
  %  OV1 list of vertex positions before performing close point collapses and
  %    edge splits
  %  OE1 list of edges before performing close point collapses and edge splits
  %  epsilon  used to collapse points and split edges
  %
  % Example:
  %  [NV,NF,OV,OE,OV1,OE1] = remesh_at_handles(V,F,C,P,BE,CE);
  %  % show a plot of what's sent to triangle
  %  subplot(3,1,1);
  %  plot([OV1(OE1(:,1),1) OV1(OE1(:,2),1)]',[OV1(OE1(:,1),2) OV1(OE1(:,2),2)]','-','LineWidth',1);
  %  hold on;
  %  plot(OV1(:,1),OV1(:,2),'.');
  %  hold off;
  %  axis equal;
  %  title('Outline + handle samples');
  %  subplot(3,1,2);
  %  plot([OV(OE(:,1),1) OV(OE(:,2),1)]',[OV(OE(:,1),2) OV(OE(:,2),2)]','-','LineWidth',1);
  %  hold on;
  %  plot(OV(:,1),OV(:,2),'.');
  %  hold off;
  %  axis equal;
  %  title('Input to triangle (collapsed points + split edges)');
  %  % show a plot of the result
  %  subplot(3,1,3);
  %  tsurf(NF,NV);
  %  hold on;
  %  plot([OV(OE(:,1),1) OV(OE(:,2),1)]',[OV(OE(:,1),2) OV(OE(:,2),2)]','-','LineWidth',2);
  %  hold off;
  %  axis equal;
  %  title('Output of triangle');

  %% THIS DOES NOT WORK FOR MESHES THAT AREN"T TOPOLOGICAL DISKS!! HOLES WILL
  %% BE CLOSED

  dim = size(V,2);
  assert(dim == size(C,2));
  % only work with 2D
  assert(dim == 2);

  % Find all edges in mesh, note internal edges are repeated
  E = sort([F(:,1) F(:,2); F(:,2) F(:,3); F(:,3) F(:,1)]')';
  % determine uniqueness of edges
  [u,m,n] = unique(E,'rows');
  % determine counts for each unique edge
  counts = accumarray(n(:), 1);
  % extract edges that only occurred once
  O = u(counts==1,:);
  % extract unique vertex indices on outline
  [u,m,n] = unique(O(:));
  % original map O = IM(O)
  IM = 1:size(V,1);
  IM(O(:)) = n;
  OV = V(u,:);
  OE = IM(O);

  samples_per_edge = 10;
  % Bone samples
  [BES,BESE] = sample_edges(C,BE,samples_per_edge);
  % Cage edge samples
  [CES,CESE] = sample_edges(C(P,:),CE,samples_per_edge);
  % Append handle samples
  OV = [OV;C(P,:);BES;CES];
  OE = [ ...
    OE; ...
    (size(OV,1)-size(BES,1)-size(CES,1))+BESE; ...
    (size(OV,1)-size(CES,1))+CESE];

  % save old values
  OV1 = OV;
  OE1 = OE;

  % phony counts to enter while loop
  prev_v_count = size(OV,1)+1;
  prev_e_count = size(OE,1)+1;
  % this epsilon is a parameter that probably needs to be tweaked and should at
  % least be exposed.
  % Use fraction of average edge length
  epsilon = mean(sqrt(sum((OV(OE(:,1),:) - OV(OE(:,2),:)).^2,2)))/8;
  % continue to collapse points and split edges until nothing more to do
  while( prev_v_count ~= size(OV,1) || prev_e_count ~= size(OE,1))
    prev_v_count = size(OV,1);
    prev_e_count = size(OE,1);
    [OV,OE] = collapse_close_points(OV,OE,epsilon);
    [OV,OE] = snap_points_to_close_edges(OV,OE,epsilon);
  end

  % Again, max area term here is a heuristic
  %sum((OV(OE(:,1)) - OV(OE(:,2))).^2,2)
  % use a multiple of min edge length
  %max_area = 2*(sqrt(3)/4)*min(sum((OV(OE(:,1),:) - OV(OE(:,2),:)).^2,2));
  max_area = 8*(sqrt(3)/4)*min(sum((OV(OE(:,1),:) - OV(OE(:,2),:)).^2,2));
  [NV,NF] = triangle(OV,OE,[],'Quality',30,'MaxArea',max_area);
end

author: alecjacobson

repo: gptoolbox


text:
%GETGABORKERNEL  Returns Gabor filter coefficients
%
%     kernel = cv.getGaborKernel('OptionName',optionValue, ...)
%
% ## Output
% * __kernel__ output kernel of the specified size and type.
%
% ## Options
% * __KSize__ Size of the filter returned `[w,h]`. default [21,21]
% * __Sigma__ Standard deviation of the gaussian envelope. default 5.0
% * __Theta__ Orientation of the normal to the parallel stripes of a Gabor
%   function. default pi/4
% * __Lambda__ Wavelength of the sinusoidal factor. default 10.0
% * __Gamma__ Spatial aspect ratio. default 0.75
% * __Psi__ Phase offset. default pi/2
% * __KType__ Type of filter coefficients. It can be `single` or `double`.
%   default `double`.
%
% For more details about gabor filter equations and parameters, see
% [Gabor filter](https://en.wikipedia.org/wiki/Gabor_filter).
%
% ## Example
% Try the `gabor_filter_gui.m` sample file to interact with the various
% parameters of the Gabor filter and visualize the result.
%
% See also: gabor, imgaborfilt
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function T = prep_trpca_data( impath, dataname, scale )

% impath = the root dir of the images

if ~exist( 'scale', 'var' ) || isempty( scale )
    scale = 1.0;
end

files = dir( [impath, '\*.bmp'] );
N = length( files );
% N = N - 1;
img = imresize( imread( [impath, '\', files(1).name], 'bmp' ), scale );
img = rgb2gray( img );
sz = size( img );   lsz = length(sz);
sz(lsz+1) = N;
X = zeros( sz );

for i = 1:N
    img = imresize( imread( [impath, '\', files(i).name], 'bmp' ), scale );
    img = rgb2gray( img );
    if lsz == 2
        X( :, :, i ) = double(img) / 255;
    else
        X( :, :, :, i ) = double(img) / 255;
    end
end

% img = imresize( imread( [impath, '\', files(N+1).name], 'bmp' ), scale );
% img = rgb2gray( img );
% data.truth = double(img) / 255;

T = tensor( X );

data.X = T;

save( ['..\..\data\',dataname], 'data' );
end

author: andrewssobral

repo: lrslibrary


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

%%
clear; close all; clc;

%% Syntax
% |[F1,V1,ind1]=triSurfSelfTriangulateBoundary(F1,V1,ind1,angleThreshold,isClosedLoop);|

%% Description
% The |triSurfSelfTriangulateBoundary| function
%

%% Examples

%% Example 1: Self-triangulate the boundary of a surface by creating triangles at sharp boundary segments
% Create test data set
w=1;
[X,Y]=ndgrid(linspace(0,w,15));
Z=ones(size(X));
C=tril(Z);
[F,V,C]=surf2patch(X,Y,Z,C); %Quads
C=vertexToFaceMeasure(F,C)>0;

logicKeep=C>0;
F=F(logicKeep,:);
C=C(logicKeep,:);
[F,V]=patchCleanUnused(F,V);

F=[F(:,[1 2 3]);F(:,[3 4 1])]; %Triangles

%%
% Get boundary curve

Eb=patchBoundary(F);
indBoundaryCurve=edgeListToCurve(Eb);
indBoundaryCurve=indBoundaryCurve(1:end-1)'; %Start=End for closed curve so remove double entry

%%
% Self triangulate
angleThreshold=(100/180)*pi;
isClosedPath=1;
[F_new,V_new,indBoundaryCurve_new]=triSurfSelfTriangulateBoundary(F,V,indBoundaryCurve,angleThreshold,isClosedPath);

%%

cFigure;
subplot(1,2,1); hold on;
title('Original');
gpatch(F,V,'kw');
plotV(V(indBoundaryCurve,:),'k-','LineWidth',3);
axisGeom; view(2);

subplot(1,2,2); hold on;
title('Self-triangulated');
gpatch(F_new,V_new,'kw');
plotV(V_new(indBoundaryCurve_new,:),'k-','LineWidth',3);
axisGeom; view(2);

drawnow;

%% Example 2: Altered shape

V(:,1)=V(:,1)-V(:,2);

%%
% Calculate mesh path angles

[F_new,V_new,indBoundaryCurve_new]=triSurfSelfTriangulateBoundary(F,V,indBoundaryCurve,angleThreshold,isClosedPath);

%%

cFigure;
subplot(1,2,1); hold on;
title('Original');
gpatch(F,V,'kw');
plotV(V(indBoundaryCurve,:),'k-','LineWidth',3);
axisGeom; view(2);

subplot(1,2,2); hold on;
title('Self-triangulated');
gpatch(F_new,V_new,'kw');
plotV(V_new(indBoundaryCurve_new,:),'k-','LineWidth',3);
axisGeom; view(2);

drawnow;

%% Example 3: Self-triangulate a part of boundary surface

%%
% Create path segment
indBoundaryCurve=indBoundaryCurve(1:15);

%%
% Calculate mesh path angles
isClosedPath=0;
[F_new,V_new,indBoundaryCurve_new]=triSurfSelfTriangulateBoundary(F,V,indBoundaryCurve,angleThreshold,isClosedPath);

%%

cFigure;
subplot(1,2,1); hold on;
title('Original');
gpatch(F,V,'kw');
plotV(V(indBoundaryCurve,:),'k-','LineWidth',3);
axisGeom; view(2);

subplot(1,2,2); hold on;
title('Self-triangulated');
gpatch(F_new,V_new,'kw');
plotV(V_new(indBoundaryCurve_new,:),'k-','LineWidth',3);
axisGeom; view(2);

drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode


text:
function [varargout]= fileutil_concatBV(file_list, varargin)
% FILEUTIL_CONCATBV - concatenate files which are stored in BrainVision format
%
% Synopsis:
%   [DAT, MRK, MNT]= fileutil_concatBV(FILE_LIST, 'Property, 'Value', ...)
%
% Arguments:
%   FILE_LIST: list of file names (no extension)
%
% Returns:
%   DAT: structure of continuous or epoched signals
%   MRK: marker structure
%   MNT: electrode montage structure
%
% Properties:
%   are passed to file_loadBV
%
% Description:
%   This function is called by file_loadBV in case the file name argument
%   is a cell array of file names. Typically there is no need to call this
%   function directly.
%
%

if ~iscell(file_list),
  file_list= {file_list};
end

T= zeros(1, length(file_list));
dataOffset = 0;
for ii= 1:length(file_list),
  [cnt, mrk, hdr]= file_loadBV(file_list{ii}, varargin{:});
  T(ii)= size(cnt.x,1);
  if ii==1,
    ccnt= cnt;
    curmrk= mrk;
  else
    if ~isequal(cnt.clab, ccnt.clab),
      warning(['inconsistent clab structure will be repaired ' ...
               'by using the intersection']);
      commonclab= intersect(cnt.clab, ccnt.clab);
      cnt= proc_selectChannels(cnt, commonclab{:});
      ccnt= proc_selectChannels(ccnt, commonclab{:});
    end
    if ~isequal(cnt.fs, ccnt.fs)
        error('inconsistent sampling rate');
    end
    ccnt.x= cat(1, ccnt.x, cnt.x);

    mrk.time= mrk.time + dataOffset*1000/cnt.fs;

%     % find markers in the loaded interval
%     inival= find(curmrk.time > skip*1000/cnt.fs & ...
%         curmrk.time <= (skip+maxlen)*1000/cnt.fs);
%     curmrk= mrk_selectEvents(curmrk, inival);
%     %let the markers start at zero
%     curmrk.time= curmrk.time - skip*1000/cnt.fs;

    curmrk = mrk_mergeMarkers(curmrk, mrk);

  end
  dataOffset = dataOffset + T(ii);

end

ccnt.T= T;
if length(file_list)>1,
  ccnt.title= [ccnt.title ' et al.'];
  ccnt.file= strcat(fileparts(ccnt.file), file_list);
end

varargout= cell(1, nargout);
varargout{1}= ccnt;
if nargout>1,
  varargout{2}= curmrk;
  if nargout>2,
    varargout{3}= hdr;
  end
end


text:
function p = addDisjointBilinearSDPcut(p)

if any(p.K.s) && any(p.variabletype == 1)
    top = 1+p.K.f+p.K.l+sum(p.K.q)+p.K.e+sum(p.K.p);
    newData = [];
    newKs = [];
    for i = 1:length(p.K.s)
        data = p.F_struc(top:top+p.K.s(i)^2-1,:);
        used = find(any(data(:,2:end)));
        if ~any(p.variabletype(used))
            % This is a linear SDP F(x) >= 0
            possible_t = setdiff(p.linears,used);
            weCanDoIt = 1;
            for j = 1:length(possible_t)
                t = possible_t(j);
                generateMonomials = t_times_x(t,used,p);
                if ~any(isnan(generateMonomials))
                    if ~isinf(p.lb(t))
                        % F0*(t-L) + Fi*xi*(t-L) >= 0
                        temp = data;
                        temp(:,1 + t) = temp(:,1);
                        temp(:,1) = temp(:,1)*-p.lb(t);
                        temp(:,1 + generateMonomials) = temp(:,1+used);
                        temp(:,1+used) = temp(:,1+used)*-p.lb(t);
                        newData = [newData;temp];
                        newKs = [newKs p.K.s(i)];
                    end
                    if ~isinf(p.ub(t))
                        % F0*(U-t) + Fi*xi*(U-t) >= 0
                        temp = data;
                        temp(:,1 + t) = -temp(:,1);
                        temp(:,1) = temp(:,1)*p.ub(t);
                        temp(:,1 + generateMonomials) = -temp(:,1+used);
                        temp(:,1+used) = temp(:,1+used)*p.ub(t);
                        newData = [newData;temp];
                        newKs = [newKs p.K.s(i)];
                    end
                end
            end
        end
        top = top + p.K.s(i)^2;
    end
    if ~isempty(newData);
        p.F_struc = [p.F_struc;newData];
        p.K.s = [p.K.s newKs];
    end
end

function index = t_times_x(t,x,p)
index = [];
for i = 1:length(x)
    xi = x(i);
    p_t = p.monomtable(t,:);
    p_x = p.monomtable(xi,:);
    j = findrows(p.monomtable,p_t+p_x);
    if isempty(j)
        index(i) = nan;
    else
        index(i) = j;
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [suspicious_index abof] = fastABOD(A,n_k)
%
% Angle Based Outlier Detection
% Authors: Hans-Peter, Kriegel Matthias, Schubert Arthur Zimek
% Original paper :
% Angle-Based Outlier Detection in High-dimensional Data In KDD2008
% Website : http://www.dbs.ifi.lmu.de/
% e-mail : {kriegel,schubert,zimek}@dbs.ifi.lmu.de
% Programmer: Yi-Ren Yeh
% modified by: Wei-Chih Lai
%
% Time Complexity : O(n^2 + n * n_k^2)
%
% Inputs
%   A: Represent NxM data
%      N is number of instance
%      M is number of feature
%   n_k: 1x1 integer (0, N]
%      iterator time for approximate abof value
%
% Outputs
%   abof: Nx1 vector
%      suspicious abof value for each instance
%   suspicious_index: Nx1 vector
%      [~,suspicious_index]=sort(abof,'ascend')
% Note
%   1. each value in output abof was normalized to [0, 1]
%   2. if input n_k == N - 1, then the effect is equal to ABOD
%
    [A, ia, ic] = unique(A,'rows');
    instance_number = size(ia, 1);
    origin_instance_number = size(ic, 1);
    var_array = zeros(instance_number, 1);
    n_k = min(n_k, instance_number);

    for i=1:instance_number
        var_front = 0;
        var_back = 0;
        denominator = 0;
        Temp = repmat(A(i, :), instance_number, 1) - A;
        Temp = sum(Temp .^ 2, 2);
        [~, index] = sort(Temp);
        index = index(2:n_k);
        index = index';
        count = 0;
        for j=index
            count = count + 1;
            for k=index(count+1:end)
                vector1 = A(j, :) - A(i, :);
                vector2 = A(k, :) - A(i, :);
                norm_vector1Xnorm_vector2 = norm(vector1) * norm(vector2);
                vector1Xvector2T = vector1 * vector2';
                var_front = var_front + (1 / norm_vector1Xnorm_vector2) * (vector1Xvector2T / (norm_vector1Xnorm_vector2 ^ 2)) ^ 2;
                var_back = var_back + (vector1Xvector2T / norm_vector1Xnorm_vector2 ^ 3);
                denominator = denominator + (1 / norm_vector1Xnorm_vector2);
            end
        end
        var_array(i) = var_front / denominator - (var_back / denominator) ^ 2;
    end

    min_var_array = min(var_array);
    abof = (var_array - min_var_array) / (max(var_array) - min_var_array);
    origin_abof = zeros(origin_instance_number, 1);
    for i=1:origin_instance_number
        origin_abof(i, 1) = abof(ic(i, 1), 1);
    end
    abof = origin_abof;
    [~, suspicious_index] = sort(abof);

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
function this = rmfield(this, vals);

    eval( [ 'this.EEG.' vals '=[];' ] );

author: PatternRecognition

repo: OpenBMI


text:
function iprime = primer ( n )

%*****************************************************************************80
%
%% PRIMER computes the prime numbers up to a given limit.
%
%  Discussion:
%
%    PRIMER returns the results of its computations in the vector
%    IPRIME.  IPRIME(I) is -1 if the number I is not prime, and
%    1 if I is prime.
%
%    The algorithm is a simple-minded sieve of Eristothenes, with
%    no attempt at efficiency.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of IPRIME, and the maximum
%    value that will be considered.
%
%    Output, integer IPRIME(N), records the results for each integer.
%    IPRIME(I) = -1 if I is not prime, and IPRIME(I) = 1 if I is
%    prime.  By convention, IPRIME(1) will be set to -1.
%

%
%  IPRIME(I) = 0 means we don't know if I is prime.
%
  iprime(1:n) = 0;
%
%  By convention, 1 is not prime.
%
  iprime(1) = -1;
  next = 1;
%
%  Examine the integers in order.
%
  for next = 2 : n

    if ( iprime(next) == 0 )
      iprime(next) = 1;
      for i = 2 * next : next : n
        iprime(i) = -1;
      end
    end

  end

  return
end


text:
function test_bug2370

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_filetype ft_read_headshape ft_write_headshape read_ply write_ply

filename = [tempname '.ply'];

[pos, tri] = mesh_sphere(162);

% combine them in a structure
mesh1.pos = pos;
mesh1.tri = tri;

ft_write_headshape(filename, mesh1, 'format', 'ply');
mesh2 = ft_read_headshape(filename);
assert(isequal(mesh1.tri, mesh2.tri));
ratio = mesh1.pos(:)\mesh2.pos(:);
assert(ratio>0.99 && ratio<1.01);
delete(filename);

ft_write_headshape(filename, mesh1, 'format', 'ply_ascii');
mesh2 = ft_read_headshape(filename);
assert(isequal(mesh1.tri, mesh2.tri));
ratio = mesh1.pos(:)\mesh2.pos(:);
assert(ratio>0.99 && ratio<1.01);
delete(filename);

ft_write_headshape(filename, mesh1, 'format', 'ply_binary');
mesh2 = ft_read_headshape(filename);
assert(isequal(mesh1.tri, mesh2.tri));
ratio = mesh1.pos(:)\mesh2.pos(:);
assert(ratio>0.99 && ratio<1.01);
delete(filename);

author: fieldtrip

repo: fieldtrip


text:
%% Example demonstrating MOL-WENO3-LF and MOL-WENO5-LF schemes
%% from J. Comp. Phys. 126, pp. 202-228 (1996) by Jiang and Shu
%% "Efficient Implementation of Weighted ENO Schemes"
function WENO_Buckley_Leverett
close all; clc;
global dx

% Spatial variable interval (-1, 1) from (8.1) example 2
x = linspace(-1.0, 1.0, 80);
N = length(x);
dx = x(2)-x(1);

% Initial conditions for 1D Buckey-Leverett problem
u0(1:N) = 0.0;
u0((N/4):(N/2)) = 1.0;

% Time interval (as in Example 2)
t = linspace(0.0, 0.4, 100);

% Solve MOL ODEs system for WENO3-LF and WENO5-LF schemes with RKF45 solver
tic;
[ZZ, u1] = ode45(@MOL_WENO3, t, u0);
[ZZ, u2] = ode45(@MOL_WENO5, t, u0);
toc

% Plot each solution
for I=1:100
    plot(x, u0, 'b-', x, u1(I,:), 'ro', x, u2(I,:), 'k-', 'LineWidth', 2);
    legend('Initial', 'MOL + WENO3-LF', 'MOL + WENO5-LF', 'Location', 'NorthEast');
    axis([x(1) x(end) 0 1.1]);
    grid on;
    pause(0.01)
    drawnow;
end

%% Buckey-Leverett flux F(u)
function res = Flux(u)
res = 4*u.^2./(4*u.^2+(1-u).^2);

%% Exact expression for partial derivative from F(u) by u
function res = Jacobian(u)
res = 8*u.*(1-u)./(5*u.^2-2*u+1).^2;

%% MOL ODEs for semi-discrete Buckley-Leverett equation with WENO3 flux
%% reconstruction
function [res] = MOL_WENO3(t, u)
global dx

N = length(u);   I = 2:(N+1);

% Preallocate memory and add ghost nodes at the ends
U(1) = u(1); U(I) = u(I-1); U(N+2) = u(end); U(N+3) = u(end);
alpha1 = zeros(size(U)); alpha2 = zeros(size(U));
R_minus = zeros(size(U)); R_plus = zeros(size(U));

% Lax-Friedrichs (LF) flux splitting
a = max(abs(Jacobian(U)));
Fp = 0.5*( Flux(U) + a*U ); Fm = 0.5*( Flux(U) - a*U );

% WENO3 "right" flux reconstruction
alpha1(I) = (1/3)./(eps + (Fp(I)-Fp(I-1)).^2).^2;   alpha2(I) = (2/3)./(eps + (Fp(I+1)-Fp(I)).^2).^2;
omega1 = alpha1./(alpha1 + alpha2);  omega2 = alpha2./(alpha1 + alpha2);
R_plus(I) = omega1(I).*(3/2*Fp(I) - 1/2*Fp(I-1)) + omega2(I).*(1/2*Fp(I) + 1/2*Fp(I+1));

% WENO3 "left" flux reconstruction
alpha1(I) = (1/3)./(eps + (Fm(I+2)-Fm(I+1)).^2).^2;   alpha2(I) = (2/3)./(eps + (Fm(I+1)-Fm(I)).^2).^2;
omega1 = alpha1./(alpha1 + alpha2);  omega2 = alpha2./(alpha1 + alpha2);
R_minus(I) = omega1(I).*(3/2*Fm(I+1) - 1/2*Fm(I+2)) + omega2(I).*(1/2*Fm(I) + 1/2*Fm(I+1));

% Combine fluxes and find finite volume spatial derivative
res(I-1) = -(R_plus(I)+R_minus(I)-R_plus(I-1)-R_minus(I-1))/dx;
res = res';

%% MOL ODEs for semi-discrete Buckley-Leverett equation with WENO5 flux
%% reconstruction
function [res] = MOL_WENO5(t, u)
global dx

N = length(u);   I = 3:(N+2);

% Preallocate memory and add ghost nodes at the ends
U(1) = u(1); U(2) = u(1); U(I) = u(I-2); U(N+3) = u(end); U(N+4) = u(end); U(N+5) = u(end);
alpha1 = zeros(size(U)); alpha2 = zeros(size(U)); alpha3 = zeros(size(U));
beta1 = zeros(size(U)); beta2 = zeros(size(U)); beta3 = zeros(size(U));
R_minus = zeros(size(U)); R_plus = zeros(size(U));

% Lax-Friedrichs (LF) flux splitting
a = max(abs(Jacobian(U)));
Fp = 0.5*( Flux(U) + a*U ); Fm = 0.5*( Flux(U) - a*U );

% WENO5 "right" flux reconstruction
beta1(I) = (13.0/12.0)*(Fp(I) - 2.0*Fp(I+1) + Fp(I+2)).^2 ...
         + (1.0/4.0)*(3.0*Fp(I) - 4.0*Fp(I+1) + Fp(I+2)).^2;
beta2(I) = (13.0/12.0)*(Fp(I-1) - 2.0*Fp(I) + Fp(I+1)).^2 ...
         + (1.0/4.0)*(Fp(I-1) - Fp(I+1)).^2;
beta3(I) = (13.0/12.0)*(Fp(I-2) - 2.0*Fp(I-1) + Fp(I)).^2 ...
         + (1.0/4.0)*(Fp(I-2) - 4.0*Fp(I-1) + 3.0*Fp(I)).^2;

alpha1(I) = (3.0/10.0)./(eps + beta1(I)).^2;
alpha2(I) = (3.0/5.0)./(eps + beta2(I)).^2;
alpha3(I) = (1.0/10.0)./(eps + beta3(I)).^2;

omega1 = alpha1./(alpha1 + alpha2 + alpha3);
omega2 = alpha2./(alpha1 + alpha2 + alpha3);
omega3 = alpha3./(alpha1 + alpha2 + alpha3);

R_plus(I) = omega1(I).*(1.0/3.0*Fp(I) + 5.0/6.0*Fp(I+1) - 1.0/6.0*Fp(I+2)) ...
          + omega2(I).*(-1.0/6.0*Fp(I-1) + 5.0/6.0*Fp(I) + 1.0/3.0*Fp(I+1)) ...
          + omega3(I).*(1.0/3.0*Fp(I-2) - 7.0/6.0*Fp(I-1) + 11.0/6.0*Fp(I));

% WENO5 "left" flux reconstruction
beta1(I) = (13.0/12.0)*(Fm(I+1) - 2.0*Fm(I+2) + Fm(I+3)).^2 ...
         + (1.0/4.0)*(3.0*Fm(I+1) - 4.0*Fm(I+2) + Fm(I+3)).^2;
beta2(I) = (13.0/12.0)*(Fm(I) - 2.0*Fm(I+1) + Fm(I+2)).^2 ...
         + (1.0/4.0)*(Fm(I) - Fm(I+2)).^2;
beta3(I) = (13.0/12.0)*(Fm(I-1) - 2.0*Fm(I) + Fm(I+1)).^2 ...
         + (1.0/4.0)*(Fm(I-1) - 4.0*Fm(I) + 3.0*Fm(I+1)).^2;

alpha1 = (1.0/10.0)./(eps + beta1).^2;
alpha2 = (3.0/5.0)./(eps + beta2).^2;
alpha3 = (3.0/10.0)./(eps + beta3).^2;

omega1 = alpha1./(alpha1 + alpha2 + alpha3);
omega2 = alpha2./(alpha1 + alpha2 + alpha3);
omega3 = alpha3./(alpha1 + alpha2 + alpha3);

R_minus(I) = omega1(I).*(1.0/3.0*Fm(I+3) - 7.0/6.0*Fm(I+2) + 11.0/6.0*Fm(I+1)) ...
           + omega2(I).*(-1.0/6.0*Fm(I+2) + 5.0/6.0*Fm(I+1) + 1.0/3.0*Fm(I)) ...
           + omega3(I).*(1.0/3.0*Fm(I+1) + 5.0/6.0*Fm(I) - 1.0/6.0*Fm(I-1));

% Combine fluxes and find finite volume spatial derivative
res(I-2) = -(R_plus(I)+R_minus(I)-R_plus(I-1)-R_minus(I-1))/dx;
res = res';

author: Sable

repo: mcbench-benchmarks


text:
function indx = whichInterval(dom, x, direction)
%WHICHINTERVAL   Determine which interval a point lies in.
%   INDX = WHICHINTERVAL(DOM, X) returns a matrix of size(X) whos j,k entry is a
%   positive integer denoting which subinterval of the domain DOM (which should
%   be a sorted real-valued vector) the real part of the point X(j,k) is
%   positioned. INDX(j,k) = -/+INF if X(j,k) < DOM(1) or X(j,k) > DOM(end),
%   respectively.
%
%   WHICHINTERVAL(DOM, X, DIRECTION) determines which interval those points in X
%   which lie on the breakpoints of DOM should be counted in. DIRECTION = -1
%   (the default) places such points in the interval to the left, whilst
%   DIRECTION = 1 places them to the right.

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

numInts = numel(dom) - 1;
xReal = real(x);
indx = NaN(size(x));

% Points to the left of the domain:
indx(xReal < dom(1)) = -inf;

% Deal with the points that lie on breakpoints:
if ( nargin > 2 && direction > 0 )
    mygt = @(x, y) x >= y;
    mylt = @(x, y) x < y;
else
    mygt = @(x, y) x > y;
    mylt = @(x, y) x <= y;
end

indx(xReal == dom(1)) = 1;
% Points within the domain:
for j = 1:numInts
    indx( mygt(xReal, dom(j)) & mylt(xReal, dom(j+1)) ) = j;
end
indx(xReal == dom(end)) = numInts;

% Points to the right of the domain:
indx(xReal > dom(end)) = inf;

end

author: chebfun


text:
function [X_norm, mu, sigma] = featureNormalize(X)
%FEATURENORMALIZE Normalizes the features in X
%   FEATURENORMALIZE(X) returns a normalized version of X where
%   the mean value of each feature is 0 and the standard deviation
%   is 1. This is often a good preprocessing step to do when
%   working with learning algorithms.

% You need to set these values correctly
X_norm = X;
mu = zeros(1, size(X, 2));
sigma = zeros(1, size(X, 2));

% ====================== YOUR CODE HERE ======================
% Instructions: First, for each feature dimension, compute the mean
%               of the feature and subtract it from the dataset,
%               storing the mean value in mu. Next, compute the
%               standard deviation of each feature and divide
%               each feature by it's standard deviation, storing
%               the standard deviation in sigma.
%
%               Note that X is a matrix where each column is a
%               feature and each row is an example. You need
%               to perform the normalization separately for
%               each feature.
%
% Hint: You might find the 'mean' and 'std' functions useful.
%

mu = mean(X_norm);
sigma = std(X_norm);

tf_mu = X_norm - repmat(mu,length(X_norm),1);
tf_std = repmat(sigma,length(X_norm),1);

X_norm = tf_mu ./ tf_std;

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

end

author: Borye

path: github-repos/MATLAB/Borye-machine-learning-coursera-1/machine-learning-coursera-1-033fdc2e6da393eeb1179a09aafe92362021effb/Week 2 Assignments/Linear Regression with Multiple Variables/mlclass-ex1/featureNormalize.m