text:
function [grad,err,finaldelta] = gradest(fun,x0)
% gradest: estimate of the gradient vector of an analytical function of n variables
% usage: [grad,err,finaldelta] = gradest(fun,x0)
%
% Uses derivest to provide both derivative estimates
% and error estimates. fun needs not be vectorized.
%
% arguments: (input)
%  fun - analytical function to differentiate. fun must
%        be a function of the vector or array x0.
%
%  x0  - vector location at which to differentiate fun
%        If x0 is an nxm array, then fun is assumed to be
%        a function of n*m variables.
%
% arguments: (output)
%  grad - vector of first partial derivatives of fun.
%        grad will be a row vector of length numel(x0).
%
%  err - vector of error estimates corresponding to
%        each partial derivative in grad.
%
%  finaldelta - vector of final step sizes chosen for
%        each partial derivative.
%
%
% Example:
%  [grad,err] = gradest(@(x) sum(x.^2),[1 2 3])
%  grad =
%      2     4     6
%  err =
%      5.8899e-15    1.178e-14            0
%
%
% Example:
%  At [x,y] = [1,1], compute the numerical gradient
%  of the function sin(x-y) + y*exp(x)
%
%  z = @(xy) sin(diff(xy)) + xy(2)*exp(xy(1))
%
%  [grad,err ] = gradest(z,[1 1])
%  grad =
%       1.7183       3.7183
%  err =
%    7.537e-14   1.1846e-13
%
%
% Example:
%  At the global minimizer (1,1) of the Rosenbrock function,
%  compute the gradient. It should be essentially zero.
%
%  rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2;
%  [g,err] = gradest(rosen,[1 1])
%  g =
%    1.0843e-20            0
%  err =
%    1.9075e-18            0
%
%
% See also: derivest, gradient
%
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 2/9/2007

% get the size of x0 so we can reshape
% later.
sx = size(x0);

% total number of derivatives we will need to take
nx = numel(x0);

grad = zeros(1,nx);
err = grad;
finaldelta = grad;
for ind = 1:nx
  [grad(ind),err(ind),finaldelta(ind)] = lf_derivest( ...
    @(xi) fun(swapelement(x0,ind,xi)), ...
    x0(ind),'deriv',1,'vectorized','no', ...
    'methodorder',2);
end

end % mainline function end

% =======================================
%      sub-functions
% =======================================
function vec = swapelement(vec,ind,val)
% swaps val as element ind, into the vector vec
vec(ind) = val;

end % sub-function end

function [der,errest,finaldelta] = lf_derivest(fun,x0,varargin)
% DERIVEST: estimate the n'th derivative of fun at x0, provide an error estimate
% usage: [der,errest] = DERIVEST(fun,x0)  % first derivative
% usage: [der,errest] = DERIVEST(fun,x0,prop1,val1,prop2,val2,...)
%
% Derivest will perform numerical differentiation of an
% analytical function provided in fun. It will not
% differentiate a function provided as data. Use gradient
% for that purpose, or differentiate a spline model.
%
% The methods used by DERIVEST are finite difference
% approximations of various orders, coupled with a generalized
% (multiple term) Romberg extrapolation. This also yields
% the error estimate provided. DERIVEST uses a semi-adaptive
% scheme to provide the best estimate that it can by its
% automatic choice of a differencing interval.
%
% Finally, While I have not written this function for the
% absolute maximum speed, speed was a major consideration
% in the algorithmic design. Maximum accuracy was my main goal.
%
%
% Arguments (input)
%  fun - function to differentiate. May be an inline function,
%        anonymous, or an m-file. fun will be sampled at a set
%        of distinct points for each element of x0. If there are
%        additional parameters to be passed into fun, then use of
%        an anonymous function is recommended.
%
%        fun should be vectorized to allow evaluation at multiple
%        locations at once. This will provide the best possible
%        speed. IF fun is not so vectorized, then you MUST set
%        'vectorized' property to 'no', so that derivest will
%        then call your function sequentially instead.
%
%        Fun is assumed to return a result of the same
%        shape as its input x0.
%
%  x0  - scalar, vector, or array of points at which to
%        differentiate fun.
%
% Additional inputs must be in the form of property/value pairs.
%  Properties are character strings. They may be shortened
%  to the extent that they are unambiguous. Properties are
%  not case sensitive. Valid property names are:
%
%  'DerivativeOrder', 'MethodOrder', 'Style', 'RombergTerms'
%  'FixedStep', 'MaxStep'
%
%  All properties have default values, chosen as intelligently
%  as I could manage. Values that are character strings may
%  also be unambiguously shortened. The legal values for each
%  property are:
%
%  'DerivativeOrder' - specifies the derivative order estimated.
%        Must be a positive integer from the set [1,2,3,4].
%
%        DEFAULT: 1 (first derivative of fun)
%
%  'MethodOrder' - specifies the order of the basic method
%        used for the estimation.
%
%        For 'central' methods, must be a positive integer
%        from the set [2,4].
%
%        For 'forward' or 'backward' difference methods,
%        must be a positive integer from the set [1,2,3,4].
%
%        DEFAULT: 4 (a second order method)
%
%        Note: higher order methods will generally be more
%        accurate, but may also suffere more from numerical
%        problems.
%
%        Note: First order methods would usually not be
%        recommended.
%
%  'Style' - specifies the style of the basic method
%        used for the estimation. 'central', 'forward',
%        or 'backwards' difference methods are used.
%
%        Must be one of 'Central', 'forward', 'backward'.
%
%        DEFAULT: 'Central'
%
%        Note: Central difference methods are usually the
%        most accurate, but sometiems one must not allow
%        evaluation in one direction or the other.
%
%  'RombergTerms' - Allows the user to specify the generalized
%        Romberg extrapolation method used, or turn it off
%        completely.
%
%        Must be a positive integer from the set [0,1,2,3].
%
%        DEFAULT: 2 (Two Romberg terms)
%
%        Note: 0 disables the Romberg step completely.
%
%  'FixedStep' - Allows the specification of a fixed step
%        size, preventing the adaptive logic from working.
%        This will be considerably faster, but not necessarily
%        as accurate as allowing the adaptive logic to run.
%
%        DEFAULT: []
%
%        Note: If specified, 'FixedStep' will define the
%        maximum excursion from x0 that will be used.
%
%  'Vectorized' - Derivest will normally assume that your
%        function can be safely evaluated at multiple locations
%        in a single call. This would minimize the overhead of
%        a loop and additional function call overhead. Some
%        functions are not easily vectorizable, but you may
%        (if your matlab release is new enough) be able to use
%        arrayfun to accomplish the vectorization.
%
%        When all else fails, set the 'vectorized' property
%        to 'no'. This will cause derivest to loop over the
%        successive function calls.
%
%        DEFAULT: 'yes'
%
%
%  'MaxStep' - Specifies the maximum excursion from x0 that
%        will be allowed, as a multiple of x0.
%
%        DEFAULT: 100
%
%  'StepRatio' - Derivest uses a proportionally cascaded
%        series of function evaluations, moving away from your
%        point of evaluation. The StepRatio is the ratio used
%        between sequential steps.
%
%        DEFAULT: 2.0000001
%
%        Note: use of a non-integer stepratio is intentional,
%        to avoid integer multiples of the period of a periodic
%        function under some circumstances.
%
%
% See the document DERIVEST.pdf for more explanation of the
% algorithms behind the parameters of DERIVEST. In most cases,
% I have chosen good values for these parameters, so the user
% should never need to specify anything other than possibly
% the DerivativeOrder. I've also tried to make my code robust
% enough that it will not need much. But complete flexibility
% is in there for your use.
%
%
% Arguments: (output)
%  der - derivative estimate for each element of x0
%        der will have the same shape as x0.
%
%  errest - 95% uncertainty estimate of the derivative, such that
%
%        abs(der(j) - f'(x0(j))) < erest(j)
%
%  finaldelta - The final overall stepsize chosen by DERIVEST
%
%
% Example usage:
%  First derivative of exp(x), at x == 1
%   [d,e]=derivest(@(x) exp(x),1)
%   d =
%       2.71828182845904
%
%   e =
%       1.02015503167879e-14
%
%  True derivative
%   exp(1)
%   ans =
%       2.71828182845905
%
% Example usage:
%  Third derivative of x.^3+x.^4, at x = [0,1]
%   derivest(@(x) x.^3 + x.^4,[0 1],'deriv',3)
%   ans =
%       6       30
%
%  True derivatives: [6,30]
%
%
% See also: gradient
%
%
% Author: John D'Errico
% e-mail: woodchips@rochester.rr.com
% Release: 1.0
% Release date: 12/27/2006

par.DerivativeOrder = 1;
par.MethodOrder = 4;
par.Style = 'central';
par.RombergTerms = 2;
par.FixedStep = [];
par.MaxStep = 100;
% setting a default stepratio as a non-integer prevents
% integer multiples of the initial point from being used.
% In turn that avoids some problems for periodic functions.
par.StepRatio = 2.0000001;
par.NominalStep = [];
par.Vectorized = 'yes';

na = length(varargin);
if (rem(na,2)==1)
  error 'Property/value pairs must come as PAIRS of arguments.'
elseif na>0
  par = parse_pv_pairs(par,varargin);
end
par = check_params(par);

% Was fun a string, or an inline/anonymous function?
if (nargin<1)
  help derivest
  return
elseif isempty(fun)
  error 'fun was not supplied.'
elseif ischar(fun)
  % a character function name
  fun = str2func(fun);
end

% no default for x0
if (nargin<2) || isempty(x0)
  error 'x0 was not supplied'
end
par.NominalStep = max(x0,0.02);

% was a single point supplied?
nx0 = size(x0);
n = prod(nx0);

% Set the steps to use.
if isempty(par.FixedStep)
  % Basic sequence of steps, relative to a stepsize of 1.
  delta = par.MaxStep*par.StepRatio .^(0:-1:-25)';
  ndel = length(delta);
else
  % Fixed, user supplied absolute sequence of steps.
  ndel = 3 + ceil(par.DerivativeOrder/2) + ...
     par.MethodOrder + par.RombergTerms;
  if par.Style(1) == 'c'
    ndel = ndel - 2;
  end
  delta = par.FixedStep*par.StepRatio .^(-(0:(ndel-1)))';
end

% generate finite differencing rule in advance.
% The rule is for a nominal unit step size, and will
% be scaled later to reflect the local step size.
fdarule = 1;
switch par.Style
  case 'central'
    % for central rules, we will reduce the load by an
    % even or odd transformation as appropriate.
    if par.MethodOrder==2
      switch par.DerivativeOrder
        case 1
          % the odd transformation did all the work
          fdarule = 1;
        case 2
          % the even transformation did all the work
          fdarule = 2;
        case 3
          % the odd transformation did most of the work, but
          % we need to kill off the linear term
          fdarule = [0 1]/fdamat(par.StepRatio,1,2);
        case 4
          % the even transformation did most of the work, but
          % we need to kill off the quadratic term
          fdarule = [0 1]/fdamat(par.StepRatio,2,2);
      end
    else
      % a 4th order method. We've already ruled out the 1st
      % order methods since these are central rules.
      switch par.DerivativeOrder
        case 1
          % the odd transformation did most of the work, but
          % we need to kill off the cubic term
          fdarule = [1 0]/fdamat(par.StepRatio,1,2);
        case 2
          % the even transformation did most of the work, but
          % we need to kill off the quartic term
          fdarule = [1 0]/fdamat(par.StepRatio,2,2);
        case 3
          % the odd transformation did much of the work, but
          % we need to kill off the linear & quintic terms
          fdarule = [0 1 0]/fdamat(par.StepRatio,1,3);
        case 4
          % the even transformation did much of the work, but
          % we need to kill off the quadratic and 6th order terms
          fdarule = [0 1 0]/fdamat(par.StepRatio,2,3);
      end
    end
  case {'forward' 'backward'}
    % These two cases are identical, except at the very end,
    % where a sign will be introduced.

    % No odd/even trans, but we already dropped
    % off the constant term
    if par.MethodOrder==1
      if par.DerivativeOrder==1
        % an easy one
        fdarule = 1;
      else
        % 2:4
        v = zeros(1,par.DerivativeOrder);
        v(par.DerivativeOrder) = 1;
        fdarule = v/fdamat(par.StepRatio,0,par.DerivativeOrder);
      end
    else
      % par.MethodOrder methods drop off the lower order terms,
      % plus terms directly above DerivativeOrder
      v = zeros(1,par.DerivativeOrder + par.MethodOrder - 1);
      v(par.DerivativeOrder) = 1;
      fdarule = v/fdamat(par.StepRatio,0,par.DerivativeOrder+par.MethodOrder-1);
    end

    % correct sign for the 'backward' rule
    if par.Style(1) == 'b'
      fdarule = -fdarule;
    end

end % switch on par.style (generating fdarule)
nfda = length(fdarule);

% will we need fun(x0)?
if (rem(par.DerivativeOrder,2) == 0) || ~strncmpi(par.Style,'central',7)
  if strcmpi(par.Vectorized,'yes')
    f_x0 = fun(x0);
  else
    % not vectorized, so loop
    f_x0 = zeros(size(x0));
    for j = 1:numel(x0)
      f_x0(j) = fun(x0(j));
    end
  end
else
  f_x0 = [];
end

% Loop over the elements of x0, reducing it to
% a scalar problem. Sorry, vectorization is not
% complete here, but this IS only a single loop.
der = zeros(nx0);
errest = der;
finaldelta = der;
for i = 1:n
  x0i = x0(i);
  h = par.NominalStep(i);

  % a central, forward or backwards differencing rule?
  % f_del is the set of all the function evaluations we
  % will generate. For a central rule, it will have the
  % even or odd transformation built in.
  if par.Style(1) == 'c'
    % A central rule, so we will need to evaluate
    % symmetrically around x0i.
    if strcmpi(par.Vectorized,'yes')
      f_plusdel = fun(x0i+h*delta);
      f_minusdel = fun(x0i-h*delta);
    else
      % not vectorized, so loop
      f_minusdel = zeros(size(delta));
      f_plusdel = zeros(size(delta));
      for j = 1:numel(delta)
        f_plusdel(j) = fun(x0i+h*delta(j));
        f_minusdel(j) = fun(x0i-h*delta(j));
      end
    end

    if ismember(par.DerivativeOrder,[1 3])
      % odd transformation
      f_del = (f_plusdel - f_minusdel)/2;
    else
      f_del = (f_plusdel + f_minusdel)/2 - f_x0(i);
    end
  elseif par.Style(1) == 'f'
    % forward rule
    % drop off the constant only
    if strcmpi(par.Vectorized,'yes')
      f_del = fun(x0i+h*delta) - f_x0(i);
    else
      % not vectorized, so loop
      f_del = zeros(size(delta));
      for j = 1:numel(delta)
        f_del(j) = fun(x0i+h*delta(j)) - f_x0(i);
      end
    end
  else
    % backward rule
    % drop off the constant only
    if strcmpi(par.Vectorized,'yes')
      f_del = fun(x0i-h*delta) - f_x0(i);
    else
      % not vectorized, so loop
      f_del = zeros(size(delta));
      for j = 1:numel(delta)
        f_del(j) = fun(x0i-h*delta(j)) - f_x0(i);
      end
    end
  end

  % check the size of f_del to ensure it was properly vectorized.
  f_del = f_del(:);
  if length(f_del)~=ndel
    error 'fun did not return the correct size result (fun must be vectorized)'
  end

  % Apply the finite difference rule at each delta, scaling
  % as appropriate for delta and the requested DerivativeOrder.
  % First, decide how many of these estimates we will end up with.
  ne = ndel + 1 - nfda - par.RombergTerms;

  % Form the initial derivative estimates from the chosen
  % finite difference method.
  der_init = vec2mat(f_del,ne,nfda)*fdarule.';

  % scale to reflect the local delta
  der_init = der_init(:)./(h*delta(1:ne)).^par.DerivativeOrder;

  % Each approximation that results is an approximation
  % of order par.DerivativeOrder to the desired derivative.
  % Additional (higher order, even or odd) terms in the
  % Taylor series also remain. Use a generalized (multi-term)
  % Romberg extrapolation to improve these estimates.
  switch par.Style
    case 'central'
      rombexpon = 2*(1:par.RombergTerms) + par.MethodOrder - 2;
    otherwise
      rombexpon = (1:par.RombergTerms) + par.MethodOrder - 1;
  end
  [der_romb,errors] = rombextrap(par.StepRatio,der_init,rombexpon);

  % Choose which result to return

  % first, trim off the
  if isempty(par.FixedStep)
    % trim off the estimates at each end of the scale
    nest = length(der_romb);
    switch par.DerivativeOrder
      case {1 2}
        trim = [1 2 nest-1 nest];
      case 3
        trim = [1:4 nest+(-3:0)];
      case 4
        trim = [1:6 nest+(-5:0)];
    end

    [der_romb,tags] = sort(der_romb);

    der_romb(trim) = [];
    tags(trim) = [];
    errors = errors(tags);
    trimdelta = delta(tags);

    [errest(i),ind] = min(errors);

    finaldelta(i) = h*trimdelta(ind);
    der(i) = der_romb(ind);
  else
    [errest(i),ind] = min(errors);
    finaldelta(i) = h*delta(ind);
    der(i) = der_romb(ind);
  end
end

end % mainline end

% ============================================
% subfunction - romberg extrapolation
% ============================================
function [der_romb,errest] = rombextrap(StepRatio,der_init,rombexpon)
% do romberg extrapolation for each estimate
%
%  StepRatio - Ratio decrease in step
%  der_init - initial derivative estimates
%  rombexpon - higher order terms to cancel using the romberg step
%
%  der_romb - derivative estimates returned
%  errest - error estimates
%  amp - noise amplification factor due to the romberg step

srinv = 1/StepRatio;

% do nothing if no romberg terms
nexpon = length(rombexpon);
rmat = ones(nexpon+2,nexpon+1);
switch nexpon
  case 0
    % rmat is simple: ones(2,1)
  case 1
    % only one romberg term
    rmat(2,2) = srinv^rombexpon;
    rmat(3,2) = srinv^(2*rombexpon);
  case 2
    % two romberg terms
    rmat(2,2:3) = srinv.^rombexpon;
    rmat(3,2:3) = srinv.^(2*rombexpon);
    rmat(4,2:3) = srinv.^(3*rombexpon);
  case 3
    % three romberg terms
    rmat(2,2:4) = srinv.^rombexpon;
    rmat(3,2:4) = srinv.^(2*rombexpon);
    rmat(4,2:4) = srinv.^(3*rombexpon);
    rmat(5,2:4) = srinv.^(4*rombexpon);
end

% qr factorization used for the extrapolation as well
% as the uncertainty estimates
[qromb,rromb] = qr(rmat,0);

% the noise amplification is further amplified by the Romberg step.
% amp = cond(rromb);

% this does the extrapolation to a zero step size.
ne = length(der_init);
rhs = vec2mat(der_init,nexpon+2,max(1,ne - (nexpon+2)));
rombcoefs = rromb\(qromb.'*rhs);
der_romb = rombcoefs(1,:).';

% uncertainty estimate of derivative prediction
s = sqrt(sum((rhs - rmat*rombcoefs).^2,1));
rinv = rromb\eye(nexpon+1);
cov1 = sum(rinv.^2,2); % 1 spare dof
errest = s.'*12.7062047361747*sqrt(cov1(1));

end % rombextrap

% ============================================
% subfunction - vec2mat
% ============================================
function mat = vec2mat(vec,n,m)
% forms the matrix M, such that M(i,j) = vec(i+j-1)
[i,j] = ndgrid(1:n,0:m-1);
ind = i+j;
mat = vec(ind);
if n==1
  mat = mat.';
end

end % vec2mat

% ============================================
% subfunction - fdamat
% ============================================
function mat = fdamat(sr,parity,nterms)
% Compute matrix for fda derivation.
% parity can be
%   0 (one sided, all terms included but zeroth order)
%   1 (only odd terms included)
%   2 (only even terms included)
% nterms - number of terms

% sr is the ratio between successive steps
srinv = 1./sr;

switch parity
  case 0
    % single sided rule
    [i,j] = ndgrid(1:nterms);
    c = 1./factorial(1:nterms);
    mat = c(j).*srinv.^((i-1).*j);
  case 1
    % odd order derivative
    [i,j] = ndgrid(1:nterms);
    c = 1./factorial(1:2:(2*nterms));
    mat = c(j).*srinv.^((i-1).*(2*j-1));
  case 2
    % even order derivative
    [i,j] = ndgrid(1:nterms);
    c = 1./factorial(2:2:(2*nterms));
    mat = c(j).*srinv.^((i-1).*(2*j));
end

end % fdamat

% ============================================
% subfunction - check_params
% ============================================
function par = check_params(par)
% check the parameters for acceptability
%
% Defaults
% par.DerivativeOrder = 1;
% par.MethodOrder = 2;
% par.Style = 'central';
% par.RombergTerms = 2;
% par.FixedStep = [];

% DerivativeOrder == 1 by default
if isempty(par.DerivativeOrder)
  par.DerivativeOrder = 1;
else
  if (length(par.DerivativeOrder)>1) || ~ismember(par.DerivativeOrder,1:4)
    error 'DerivativeOrder must be scalar, one of [1 2 3 4].'
  end
end

% MethodOrder == 2 by default
if isempty(par.MethodOrder)
  par.MethodOrder = 2;
else
  if (length(par.MethodOrder)>1) || ~ismember(par.MethodOrder,[1 2 3 4])
    error 'MethodOrder must be scalar, one of [1 2 3 4].'
  elseif ismember(par.MethodOrder,[1 3]) && (par.Style(1)=='c')
    error 'MethodOrder==1 or 3 is not possible with central difference methods'
  end
end

% style is char
valid = {'central', 'forward', 'backward'};
if isempty(par.Style)
  par.Style = 'central';
elseif ~ischar(par.Style)
  error 'Invalid Style: Must be character'
end
ind = find(strncmpi(par.Style,valid,length(par.Style)));
if (length(ind)==1)
  par.Style = valid{ind};
else
  error(['Invalid Style: ',par.Style])
end

% vectorized is char
valid = {'yes', 'no'};
if isempty(par.Vectorized)
  par.Vectorized = 'yes';
elseif ~ischar(par.Vectorized)
  error 'Invalid Vectorized: Must be character'
end
ind = find(strncmpi(par.Vectorized,valid,length(par.Vectorized)));
if (length(ind)==1)
  par.Vectorized = valid{ind};
else
  error(['Invalid Vectorized: ',par.Vectorized])
end

% RombergTerms == 2 by default
if isempty(par.RombergTerms)
  par.RombergTerms = 2;
else
  if (length(par.RombergTerms)>1) || ~ismember(par.RombergTerms,0:3)
    error 'RombergTerms must be scalar, one of [0 1 2 3].'
  end
end

% FixedStep == [] by default
if (length(par.FixedStep)>1) || (~isempty(par.FixedStep) && (par.FixedStep<=0))
  error 'FixedStep must be empty or a scalar, >0.'
end

% MaxStep == 10 by default
if isempty(par.MaxStep)
  par.MaxStep = 10;
elseif (length(par.MaxStep)>1) || (par.MaxStep<=0)
  error 'MaxStep must be empty or a scalar, >0.'
end

end % check_params

% ============================================
% Included subfunction - parse_pv_pairs
% ============================================
function params=parse_pv_pairs(params,pv_pairs)
% parse_pv_pairs: parses sets of property value pairs, allows defaults
% usage: params=parse_pv_pairs(default_params,pv_pairs)
%
% arguments: (input)
%  default_params - structure, with one field for every potential
%             property/value pair. Each field will contain the default
%             value for that property. If no default is supplied for a
%             given property, then that field must be empty.
%
%  pv_array - cell array of property/value pairs.
%             Case is ignored when comparing properties to the list
%             of field names. Also, any unambiguous shortening of a
%             field/property name is allowed.
%
% arguments: (output)
%  params   - parameter struct that reflects any updated property/value
%             pairs in the pv_array.
%
% Example usage:
% First, set default values for the parameters. Assume we
% have four parameters that we wish to use optionally in
% the function examplefun.
%
%  - 'viscosity', which will have a default value of 1
%  - 'volume', which will default to 1
%  - 'pie' - which will have default value 3.141592653589793
%  - 'description' - a text field, left empty by default
%
% The first argument to examplefun is one which will always be
% supplied.
%
%   function examplefun(dummyarg1,varargin)
%   params.Viscosity = 1;
%   params.Volume = 1;
%   params.Pie = 3.141592653589793
%
%   params.Description = '';
%   params=parse_pv_pairs(params,varargin);
%   params
%
% Use examplefun, overriding the defaults for 'pie', 'viscosity'
% and 'description'. The 'volume' parameter is left at its default.
%
%   examplefun(rand(10),'vis',10,'pie',3,'Description','Hello world')
%
% params =
%     Viscosity: 10
%        Volume: 1
%           Pie: 3
%   Description: 'Hello world'
%
% Note that capitalization was ignored, and the property 'viscosity'
% was truncated as supplied. Also note that the order the pairs were
% supplied was arbitrary.

npv = length(pv_pairs);
n = npv/2;

if n~=floor(n)
  error 'Property/value pairs must come in PAIRS.'
end
if n<=0
  % just return the defaults
  return
end

if ~isstruct(params)
  error 'No structure for defaults was supplied'
end

% there was at least one pv pair. process any supplied
propnames = fieldnames(params);
lpropnames = lower(propnames);
for i=1:n
  p_i = lower(pv_pairs{2*i-1});
  v_i = pv_pairs{2*i};

  ind = strmatch(p_i,lpropnames,'exact');
  if isempty(ind)
    ind = find(strncmp(p_i,lpropnames,length(p_i)));
    if isempty(ind)
      error(['No matching property found for: ',pv_pairs{2*i-1}])
    elseif length(ind)>1
      error(['Ambiguous property name: ',pv_pairs{2*i-1}])
    end
  end
  p_i = propnames{ind};

  % override the corresponding default in params
  params = setfield(params,p_i,v_i); %#ok

end

end % parse_pv_pairs


text:
clear; close all; clc;

%% Excel data load
WP = xlsread('H:\NeuroImage\WM_acc.xlsx', 3, 'B3:D19');
PM = xlsread('H:\NeuroImage\WM_acc.xlsx', 3, 'J3:L19');
LM = xlsread('H:\NeuroImage\WM_acc.xlsx', 3, 'R3:T19');

col=[253, 174, 97;
    171, 221, 164;
    43, 131, 186];
col=col/255;

i = [2, 0, -2];
%% boxplot [WP]
figure()
boxplot(WP, 'Notch', 'on')
set(gca, 'xtick', []);

set(findobj(gca,'type','line'),'linew',3);
h = findobj(gca,'Tag','Box');

for j=1:length(h)
    x = patch(get(h(j+i(j)),'XData'),get(h(j+i(j)),'YData'), col(j,:),'FaceAlpha',.5);
    % MODIFICATION TO SET THE LINE COLOR TO PATCH COLOR:
    h(length(h)-j+1).Color = x.FaceColor; % reordered to match
    hold on;

    x=repmat(1:3,length(WP),1);
    scatter(x(:,j),WP(:,j),'filled','MarkerEdgeColor',col(:,j), 'MarkerFaceColor',col(:,j),'jitter','on','jitterAmount',0.15);
end
% set(findobj(gca,'type','line')','linew',2);
% hLegend = legend(findall(gca,'Tag','Box'), {'BN','AN','24H'});

%% boxplot [PM]
figure()
boxplot(PM, 'Notch', 'on')
set(gca, 'xtick', []);

set(findobj(gca,'type','line'),'linew',3);
h = findobj(gca,'Tag','Box');

for j=1:length(h)
    x = patch(get(h(j+i(j)),'XData'),get(h(j+i(j)),'YData'), col(j,:),'FaceAlpha',.5);
    % MODIFICATION TO SET THE LINE COLOR TO PATCH COLOR:
    h(length(h)-j+1).Color = x.FaceColor; % reordered to match
    hold on;

    x=repmat(1:3,length(PM),1);
    scatter(x(:,j),PM(:,j),'filled','MarkerEdgeColor',col(:,j), 'MarkerFaceColor',col(:,j),'jitter','on','jitterAmount',0.15);
end
% set(findobj(gca,'type','line')','linew',2);
% hLegend = legend(findall(gca,'Tag','Box'), {'BN','AN','24H'});

%% boxplot [LM]
figure()
boxplot(LM, 'Notch', 'on')
set(gca, 'xtick', []);

set(findobj(gca,'type','line'),'linew',3);
h = findobj(gca,'Tag','Box');

for j=1:length(h)
    x = patch(get(h(j+i(j)),'XData'),get(h(j+i(j)),'YData'), col(j,:),'FaceAlpha',.5);
    % MODIFICATION TO SET THE LINE COLOR TO PATCH COLOR:
    h(length(h)-j+1).Color = x.FaceColor; % reordered to match
    hold on;

    x=repmat(1:3,length(LM),1);
    scatter(x(:,j),LM(:,j),'filled','MarkerEdgeColor',col(:,j), 'MarkerFaceColor',col(:,j),'jitter','on','jitterAmount',0.15);
end
% set(findobj(gca,'type','line')','linew',2);
% hLegend = legend(findall(gca,'Tag','Box'), {'BN','AN','24H'});

author: PatternRecognition

repo: OpenBMI


text:
function iszscored = BF_iszscored(x)
% BF_iszscored  Crude check for whether a data vector is z-scored.
%
% (~eps-close to being) z-scored.
% Used for displaying warning messages for functions that require z-scored inputs.
%
%---INPUT:
% x, the input time series (or any vector)
%
%---OUTPUT:
% iszscored, a logical with the verdict.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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/>.
% ------------------------------------------------------------------------------

% Give it a bit of numerical lee-way... Down in the 2e-14 region:
numericThreshold = 100*eps;

iszscored = ((abs(mean(x)) < numericThreshold) && (abs(std(x)-1) < numericThreshold));

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function tagS = general_series_module_tags_subset
%"general_series_module_tags_subset"
%   Return a subset of the tags specified in C.7.3.1 in PS3.3 of the 2006
%   DICOM specification required to identify different series'.
%
%   Tags are returned in a struct array with 3 fields:
%
%   Tag: String containing hex DICOM tag of a field.
%  Type: String describing type of field, with 5 options:
%         '1' Field must exist, data must exist and be valid.
%         '2' Field must exist, data can exist and be valid, or be NULL.
%         '3' Field is optional, if the field exists data can exist and be
%             valid or be NULL.
%         '1C' Field must exist under certain conditions and contain valid
%             data.
%         '2C' Field must exist under certain conditions and can contain
%             valid data or be NULL.
%Children: For sequences, a tagS with the same format as this struct array.
%
%AI 11/20/2020
%
%Usage:
%   tagS = general_series_module_tags_subset
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

tagC = {'00080060','0020000E'};

tagDecC = {524384,2097166};

fieldnameC = {'Modality','SeriesInstanceUID'};

typeC = {'1','1'};

childC = cell(1,length(tagC));

tagS = struct('tag',tagC,'fieldname',fieldnameC,'tagdec',tagDecC,'type',typeC,'children',childC);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function varargout=lin2pcma(varargin)
%
% For calling details please see v_lin2pcma.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lin2pcma.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lin2pcma.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lin2pcma(varargin{:});
else
    v_lin2pcma(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%--- help for check_property ---
%
%  INTERNAL FUNCTION
%
%
%    Other functions named check_property
%
%       dsge/check_property    generic/check_property
%

author: jmaih


text:
function LDLplot(X1, Y1, S1, C1)
%CREATEFIGURE(X1,Y1,S1,C1)
%  X1:  scatter x
%  Y1:  scatter y
%  S1:  scatter s
%  C1:  scatter c

%  Auto-generated by MATLAB on 21-Sep-2010 11:16:44
%  Copyright 2010 - 2011 MathWorks, Inc.

% Create figure
figure1 = figure;

% Create axes
axes1 = axes('Parent',figure1);
% Uncomment the following line to preserve the X-limits of the axes
% xlim(axes1,[3.6 5]);
% Uncomment the following line to preserve the Y-limits of the axes
% ylim(axes1,[1 5]);
box(axes1,'on');
hold(axes1,'all');

% Create scatter
scatter1 = scatter(X1,Y1,S1,C1,'MarkerEdgeColor',[0 0.498039215803146 0],...
    'Marker','.',...
    'Parent',axes1,...
    'DisplayName','data 1');

% Create xlabel
xlabel('LDL C at BL');

% Create ylabel
ylabel('LDL C at 12 weeks');

% Get xdata from plot
xdata1 = get(scatter1, 'xdata');
% Get ydata from plot
ydata1 = get(scatter1, 'ydata');
% Make sure data are column vectors
xdata1 = xdata1(:);
ydata1 = ydata1(:);

% Get axes ylim
axYLim1 = get(axes1, 'ylim');
% Get axes xlim
axXLim1 = get(axes1, 'xlim');

% Find the mean
xmean1 = mean(xdata1);
% Get coordinates for the mean line
meanValue1 = [xmean1 xmean1];
% Plot the mean
statLine1 = plot(meanValue1,axYLim1,'DisplayName','   x mean',...
    'Parent',axes1,...
    'Tag','mean x',...
    'LineWidth',3,...
    'LineStyle','--',...
    'Color',[1 0 0]);

% Set new line in proper position
setLineOrder(axes1, statLine1, scatter1);

% Find the mean
ymean1 = mean(ydata1);
% Get coordinates for the mean line
meanValue2 = [ymean1 ymean1];
% Plot the mean
statLine2 = plot(axXLim1,meanValue2,'DisplayName','   y mean',...
    'Parent',axes1,...
    'Tag','mean y',...
    'LineWidth',3,...
    'LineStyle','-.',...
    'Color',[1 0 0]);

% Set new line in proper position
setLineOrder(axes1, statLine2, scatter1);

% Create legend
legend(axes1,'show');

%-------------------------------------------------------------------------%
function setLineOrder(axesh1, newLine1, associatedLine1)
%SETLINEORDER(AXESH1,NEWLINE1,ASSOCIATEDLINE1)
%  Set line order
%  AXESH1:  axes
%  NEWLINE1:  new line
%  ASSOCIATEDLINE1:  associated line

% Get the axes children
hChildren = get(axesh1,'Children');
% Remove the new line
hChildren(hChildren==newLine1) = [];
% Get the index to the associatedLine
lineIndex = find(hChildren==associatedLine1);
% Reorder lines so the new line appears with associated data
hNewChildren = [hChildren(1:lineIndex-1);newLine1;hChildren(lineIndex:end)];
% Set the children:
set(axesh1,'Children',hNewChildren);

author: Sable

repo: mcbench-benchmarks


text:
function [ MINSCORE, MINTRANS] = compRoomMatchScore( HYPS, VALID_TRANS_GNDS )
%COMPROOMMATCHSCORE Summary of this function goes here
%   Detailed explanation goes here
HYP_COMP = getHypInfo(HYPS);
MINSCORE = zeros(length(HYPS), length(VALID_TRANS_GNDS));
MINTRANS = zeros(length(HYPS), length(VALID_TRANS_GNDS));
for j = 1:length(VALID_TRANS_GNDS)
%     fprintf('GND: %d/%d\n', j, length(VALID_TRANS_GNDS));
    fprintf('>>%d', j);
    tCOST = roomAlignmentMex(HYP_COMP, VALID_TRANS_GNDS(j).SAMPLE_COMP, 3, 0);
    [B,I] = min( tCOST, [], 2);
    MINSCORE(:,j) = B;
    MINTRANS(:,j) = I;
end

end

author: zouchuhang

repo: LayoutNet


text:
close all;
% function F = discriminative_texture_feature(I_TEXT,theta,verbose,colored,include_intensity,tau_diff,steps_diff,sigma_diff)

Z = imread('zebra.bmp');
figure; imagesc(Z); set(gcf,'Name','original image');
% Z0 = discriminative_texture_feature(double(Z),6,[1 2],0,0,1,500,0); %

Z1 = discriminative_texture_feature(double(Z),0,[1 2],0,1,1,500,0); % this takes a lot of time!, compare to p 64 of Thomas Brox's Phd Thesis
Z2 = discriminative_texture_feature(double(Z));% this is a rasonable approximation(depending on the application) to the previous one while being much faster

F = imread('frog.bmp');
figure; imagesc(F); set(gcf,'Name','original image');
% F0 = discriminative_texture_feature(double(F),6,[1 2],0,0,1,500,0); %
F1 = discriminative_texture_feature(double(F),0,[1 2],0,1,1,1000,0); % this takes a lot of time!, compare to p 65 of Thomas Brox's Phd Thesis

F2 = discriminative_texture_feature(double(F),2,[2],0,1,10,100,0.5); % this is a rasonable approximation(depending on the application to the previous one while being much faster

[sy sx d] = size(Z);
figure;colormap gray;
subplot(2,3,1); imagesc(Z); title('original image');
for i = 1 :5, subplot(2,3,i+1); imagesc(reshape(Z1(i,:),[sy sx])); title(sprintf('F%d',i)); end


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: Examples for grid generation
%
%==============================================================================

omega = [0,6,0,4,0,8]
m     = [3,2,2]
xc    = getCellCenteredGrid(omega(1:2),m(1));        xc = reshape(xc,1,[])
xc    = getCellCenteredGrid(omega(1:4),m(1:2));      xc = reshape(xc,[m(1:2),2])
xc    = getCellCenteredGrid(omega(1:6),m(1:3));      xc = reshape(xc,[m(1:3),3])
%==============================================================================


text:
% Launch AMS example selector
% Copywrite 2001-2013 The MathWorks, Inc.
cd \ams
example_sel2html('AMS_sel.txt')


text:
function sizes = getVarSizes(obj, inputSizes)
%GETVARSIZES  Get the size of the variables
%   SIZES = GETVARSIZES(OBJ, INPUTSIZES) computes the SIZES of the
%   DagNN variables given the size of the inputs. `inputSizes` is
%   a cell array of the type `{'inputName', inputSize, ...}`
%   Returns a cell array with sizes of all network variables.
%
%   Example, compute the storage needed for a batch size of 256 for an
%   imagenet-like network:
%   ```
%   batch_size = 256; single_num_bytes = 4;
%   input_size = [net.meta.normalization.imageSize, batch_size];
%   var_sizes = net.getVarSizes({'data', input_size});
%   fprintf('Network activations will take %.2fMiB in single.\n', ...
%     sum(prod(cell2mat(var_sizes, 1))) * single_num_bytes ./ 1024^3);
%   ```

% Copyright (C) 2015 Andrea Vedaldi, Karel Lenc.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

nv = numel(obj.vars) ;
sizes = num2cell(NaN(nv, 4),2)' ;

for i = 1:2:numel(inputSizes)
  v = obj.getVarIndex(inputSizes{i}) ;
  if isnan(v)
    error('Variable `%s` not found in the network.', inputSizes{i});
  end;
  sizes{v} = [inputSizes{i+1}(:)' ones(1, 4 - numel(inputSizes{i+1}))] ;
end

for layer = obj.layers(obj.executionOrder)
  in = layer.inputIndexes ;
  out = layer.outputIndexes ;
  sizes(out) = layer.block.getOutputSizes(sizes(in)) ;
end

repo: CREST-Release


text:
function keypointsAll = keypoints14to16(keypointsAllcell)

keypointsAll = repmat(struct('det',nan(16,3),'imgname',''),size(keypointsAllcell,1),1);

for imgidx = 1:length(keypointsAll)
    det = keypointsAllcell{imgidx};
    keypointsAll(imgidx).det(1:6,1:2) = det(1:6,1:2);
    %keypointsAll(imgidx).det(9,1:2) = det(13,1:2);
    %keypointsAll(imgidx).det(10,1:2) = det(14,1:2);
    keypointsAll(imgidx).det(9:16,1:2) = det(7:14,1:2);
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [event] = read_yokogawa_event(filename, varargin)

% READ_YOKOGAWA_EVENT reads event information from continuous,
% epoched or averaged MEG data that has been generated by the Yokogawa
% MEG system and software and allows those events to be used in
% combination with FieldTrip.
%
% Use as
%   [event] = read_yokogawa_event(filename)
%
% See also READ_YOKOGAWA_HEADER, READ_YOKOGAWA_DATA

% Copyright (C) 2005, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

event   = [];
handles = definehandles;

% get the options, the default is set below
chanindx    = ft_getopt(varargin, 'chanindx');
threshold   = ft_getopt(varargin, 'threshold');
detectflank = ft_getopt(varargin, 'detectflank');
combinebinary = ft_getopt(varargin, 'combinebinary', false);
trigshift   = ft_getopt(varargin, 'trigshift');
if isempty(combinebinary)
  combinebinary = false;
end

% ensure that the required toolbox is on the path
if ft_hastoolbox('yokogawa_meg_reader')
  % read the dataset header
  hdr = read_yokogawa_header_new(filename);
  ch_info = hdr.orig.channel_info.channel;
  type = [ch_info.type];

  % determine the trigger channels (if not specified by the user)
  if isempty(chanindx)
    chanindx = find(type==handles.TriggerChannel);
  end

  % Use the MEG Reader documentation if more detailed support is required.
  if hdr.orig.acq_type==handles.AcqTypeEvokedRaw
    % read the trigger id from all trials
    event = getYkgwHdrEvent(filename);
    % use the standard FieldTrip header for trial events
    % make an event for each trial as defined in the header
    for i=1:hdr.nTrials
      event(end+1).type     = 'trial';
      event(end  ).sample   = (i-1)*hdr.nSamples + 1;
      event(end  ).offset   = -hdr.nSamplesPre;
      event(end  ).duration =  hdr.nSamples;
      if ~isempty(value)
        event(end  ).value    =  event(i).code;
      end
    end

  % Use the MEG Reader documentation if more detailed support is required.
  elseif hdr.orig.acq_type==handles.AcqTypeEvokedAve
    % make an event for the average
    event(1).type     = 'average';
    event(1).sample   = 1;
    event(1).offset   = -hdr.nSamplesPre;
    event(1).duration =  hdr.nSamples;

  elseif hdr.orig.acq_type==handles.AcqTypeContinuousRaw
    % Events annotated by users during measurements
    bookmark_tmp = getYkgwHdrBookmark(filename);
    if ~isempty(bookmark_tmp)
      for i = 1:length(bookmark_tmp)
        event(end+1).sample = bookmark_tmp(i).sample_no + 1;
        event(end  ).type   =  'bookmarks';
        if ~isempty(bookmark_tmp(i).label)
          event(end  ).value = bookmark_tmp(i).label;
          %% 0 in default
        else
          event(end  ).value = 99;
        end
      end
    end
    clear bookmark_tmp;
  end

elseif ft_hastoolbox('yokogawa')

  % read the dataset header
  hdr = read_yokogawa_header(filename);

  % determine the trigger channels (if not specified by the user)
  if isempty(chanindx)
    chanindx = find(hdr.orig.channel_info(:,2)==handles.TriggerChannel);
  end

  if hdr.orig.acq_type==handles.AcqTypeEvokedRaw
    % read the trigger id from all trials
    fid   = fopen_or_error(filename, 'r');
    value = GetMeg160TriggerEventM(fid);
    fclose(fid);
    % use the standard FieldTrip header for trial events
    % make an event for each trial as defined in the header
    for i=1:hdr.nTrials
      event(end+1).type     = 'trial';
      event(end  ).sample   = (i-1)*hdr.nSamples + 1;
      event(end  ).offset   = -hdr.nSamplesPre;
      event(end  ).duration =  hdr.nSamples;
      if ~isempty(value)
        event(end  ).value    =  value(i);
      end
    end

  elseif hdr.orig.acq_type==handles.AcqTypeEvokedAve
    % make an event for the average
    event(1).type     = 'average';
    event(1).sample   = 1;
    event(1).offset   = -hdr.nSamplesPre;
    event(1).duration =  hdr.nSamples;

  elseif hdr.orig.acq_type==handles.AcqTypeContinuousRaw
    % the data structure does not contain events, but flank detection on the trigger channel might reveal them
    % this is done below
  end

else
  ft_error('cannot determine whether the required Yokogawa toolbox is present');
end

% read the trigger channels and detect the flanks
if ~isempty(chanindx)
  trigger = read_trigger(filename, 'header', hdr, 'denoise', false, 'chanindx', chanindx, 'detectflank', detectflank, 'threshold', threshold, 'combinebinary', combinebinary, 'trigshift', trigshift);
  % combine the triggers and the other events
  event = appendstruct(event, trigger);
end

if isempty(event)
  ft_warning('no triggers were detected, please specify the "trigindx" option');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this defines some usefull constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function handles = definehandles
handles.output = [];
handles.sqd_load_flag = false;
handles.mri_load_flag = false;
handles.NullChannel         = 0;
handles.MagnetoMeter        = 1;
handles.AxialGradioMeter    = 2;
handles.PlannerGradioMeter  = 3;
handles.RefferenceChannelMark = hex2dec('0100');
handles.RefferenceMagnetoMeter       = bitor( handles.RefferenceChannelMark, handles.MagnetoMeter );
handles.RefferenceAxialGradioMeter   = bitor( handles.RefferenceChannelMark, handles.AxialGradioMeter );
handles.RefferencePlannerGradioMeter = bitor( handles.RefferenceChannelMark, handles.PlannerGradioMeter );
handles.TriggerChannel      = -1;
handles.EegChannel          = -2;
handles.EcgChannel          = -3;
handles.EtcChannel          = -4;
handles.NonMegChannelNameLength = 32;
handles.DefaultMagnetometerSize       = (4.0/1000.0);       % Square of 4.0mm in length
handles.DefaultAxialGradioMeterSize   = (15.5/1000.0);      % Circle of 15.5mm in diameter
handles.DefaultPlannerGradioMeterSize = (12.0/1000.0);      % Square of 12.0mm in length
handles.AcqTypeContinuousRaw = 1;
handles.AcqTypeEvokedAve     = 2;
handles.AcqTypeEvokedRaw     = 3;
handles.sqd = [];
handles.sqd.selected_start  = [];
handles.sqd.selected_end    = [];
handles.sqd.axialgradiometer_ch_no      = [];
handles.sqd.axialgradiometer_ch_info    = [];
handles.sqd.axialgradiometer_data       = [];
handles.sqd.plannergradiometer_ch_no    = [];
handles.sqd.plannergradiometer_ch_info  = [];
handles.sqd.plannergradiometer_data     = [];
handles.sqd.eegchannel_ch_no   = [];
handles.sqd.eegchannel_data    = [];
handles.sqd.nullchannel_ch_no   = [];
handles.sqd.nullchannel_data    = [];
handles.sqd.selected_time       = [];
handles.sqd.sample_rate         = [];
handles.sqd.sample_count        = [];
handles.sqd.pretrigger_length   = [];
handles.sqd.matching_info   = [];
handles.sqd.source_info     = [];
handles.sqd.mri_info        = [];
handles.mri                 = [];

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ q_min, q_max, q_ave, q_area ] = q_measure ( n, z, triangle_order, ...
  triangle_num, triangle_node )

%*****************************************************************************80
%
%% Q_MEASURE determines the triangulated pointset quality measure Q.
%
%  Discussion:
%
%    The Q measure evaluates the uniformity of the shapes of the triangles
%    defined by a triangulated pointset.
%
%    For a single triangle T, the value of Q(T) is defined as follows:
%
%      TAU_IN = radius of the inscribed circle,
%      TAU_OUT = radius of the circumscribed circle,
%
%      Q(T) = 2 * TAU_IN / TAU_OUT
%        = ( B + C - A ) * ( C + A - B ) * ( A + B - C ) / ( A * B * C )
%
%    where A, B and C are the lengths of the sides of the triangle T.
%
%    The Q measure computes the value of Q(T) for every triangle T in the
%    triangulation, and then computes the standard deviation of this
%    set of values:
%
%      Q_MEASURE = min ( all T in triangulation ) Q(T)
%
%    In an ideally regular mesh, all triangles would have the same
%    equilateral shape, for which Q = 1.  In a good mesh, 0.5 < Q.
%
%    Given the 2D coordinates of a set of N nodes, stored as Z(1:2,1:N),
%    a triangulation is a list of NT triples of node indices that form
%    triangles.  Generally, a maximal triangulation is expected, namely,
%    a triangulation whose image is a planar graph, but for which the
%    addition of any new triangle would mean the graph was no longer planar.
%    A Delaunay triangulation is a maximal triangulation which maximizes
%    the minimum angle that occurs in any triangle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, pages 329-345, June 2004.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real Z(DIM_NUM,N), the points.
%
%    Input, integer TRIANGLE_ORDER, the order of the triangles.
%
%    Input, integer TRIANGLE_NUM, the number of triangles.
%
%    Input, integer TRIANGLE_NODE(TRIANGLE_ORDER,TRIANGLE_NUM), the triangulation.
%
%    Output, real ( kind = 8 ) Q_MIN, Q_MAX, the minimum and maximum values
%    of Q over all triangles.
%
%    Output, real ( kind = 8 ) Q_AVE, the average value of Q.
%
%    Output, real ( kind = 8 ) Q_AREA, the average value of Q, weighted by
%    the area of each triangle.
%
  q_min = r8_huge ( );
  q_max = - r8_huge ( );
  q_ave = 0.0;
  q_area = 0.0;
  area_total = 0.0;

  for triangle = 1 : triangle_num

    a_index = triangle_node(1,triangle);
    b_index = triangle_node(2,triangle);
    c_index = triangle_node(3,triangle);

    ab_length = sqrt ( ...
        ( z(1,a_index) - z(1,b_index) )^2 ...
      + ( z(2,a_index) - z(2,b_index) )^2 );

    bc_length = sqrt ( ...
        ( z(1,b_index) - z(1,c_index) )^2 ...
      + ( z(2,b_index) - z(2,c_index) )^2 );

    ca_length = sqrt ( ...
        ( z(1,c_index) - z(1,a_index) )^2 ...
      + ( z(2,c_index) - z(2,a_index) )^2 );

    q = ( bc_length + ca_length - ab_length ) ...
      * ( ca_length + ab_length - bc_length ) ...
      * ( ab_length + bc_length - ca_length ) ...
      / ( ab_length * bc_length * ca_length );

    x1 = z(1,triangle_node(1,triangle));
    y1 = z(2,triangle_node(1,triangle));
    x2 = z(1,triangle_node(2,triangle));
    y2 = z(2,triangle_node(2,triangle));
    x3 = z(1,triangle_node(3,triangle));
    y3 = z(2,triangle_node(3,triangle));

    area = 0.5 * abs ( x1 * ( y2 - y3 ) ...
                     + x2 * ( y3 - y1 ) ...
                     + x3 * ( y1 - y2 ) );

    q_min = min ( q_min, q );
    q_max = max ( q_max, q );
    q_ave = q_ave + q;
    q_area = q_area + q * area;

    area_total = area_total + area;

  end

  q_ave = q_ave / triangle_num;

  if ( 0.0 < area_total )
    q_area = q_area / area_total;
  else
    q_area = 0.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


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

clear all;
%close all;
clf

num=1;
den=[1 0 0];
w=logspace(-2,2,100);
[m,p]=bode(num,den,w);
loglog(w,m);
xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.47 Spacecraft frequency-response magnitude');
bodegrid;


text:
% This package contains code for the "CRF-RNN" semantic image segmentation method, published in the
% ICCV 2015 paper Conditional Random Fields as Recurrent Neural Networks. Our software is built on
% top of the Caffe deep learning library.
%
% Contact:
% Shuai Zheng (szheng@robots.ox.ac.uk), Sadeep Jayasumana (sadeep@robots.ox.ac.uk), Bernardino Romera-Paredes (bernard@robots.ox.ac.uk)
%
% Supervisor:
% Philip Torr (philip.torr@eng.ox.ac.uk)
%
% For more information about CRF-RNN please vist the project website http://crfasrnn.torr.vision.
%

caffe_path = '../caffe/';

model_def_file = 'TVG_CRFRNN_new_deploy.prototxt';
model_file = 'TVG_CRFRNN_COCO_VOC.caffemodel';

if exist(model_file, 'file') ~= 2
    error('You need a network model file. Please download our default model by running ./download_trained_model.sh');
end

use_gpu = 1; % Set this to 0 if you don't have a GPU.
gpu_id = 0;% which gpu device id you are using?

addpath(fullfile(caffe_path, 'matlab'));

caffe.reset_all();

net = tvg_matcaffe_init(use_gpu, gpu_id, model_def_file, model_file);

im = imread('input.jpg');

[h, w, d] = size(im);

if (d ~= 3)
   error('Error! Wrong depth.\n');
end

if (h > 500 || w > 500)
   error('Error! Wrong image size.\n');
end

prepared_im = tvg_prepare_image_fixed(im);

inputData = {prepared_im};
scores = net.forward(inputData);

Q = scores{1};

[dumb, pred] = max(Q, [], 3);
pred = pred';
pred = pred(1:h, 1:w);

load map.mat
imwrite(pred, map, 'output.png', 'png');

author: torrvision

repo: crfasrnn

save_path: github-repos/MATLAB/torrvision-crfasrnn


text:
function [Texp,Lexp]=lyapunov(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,ystart,ioutp);
%
%    Lyapunov exponent calcullation for ODE-system.
%
%    The alogrithm employed in this m-file for determining Lyapunov
%    exponents was proposed in
%
%         A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano,
%        "Determining Lyapunov Exponents from a Time Series," Physica D,
%        Vol. 16, pp. 285-317, 1985.
%
%    For integrating ODE system can be used any MATLAB ODE-suite methods.
% This function is a part of MATDS program - toolbox for dynamical system investigation
%    See:    http://www.math.rsu.ru/mexmat/kvm/matds/
%
%    Input parameters:
%      n - number of equation
%      rhs_ext_fcn - handle of function with right hand side of extended ODE-system.
%              This function must include RHS of ODE-system coupled with
%              variational equation (n items of linearized systems, see Example).
%      fcn_integrator - handle of ODE integrator function, for example: @ode45
%      tstart - start values of independent value (time t)
%      stept - step on t-variable for Gram-Schmidt renormalization procedure.
%      tend - finish value of time
%      ystart - start point of trajectory of ODE system.
%      ioutp - step of print to MATLAB main window. ioutp==0 - no print,
%              if ioutp>0 then each ioutp-th point will be print.
%
%    Output parameters:
%      Texp - time values
%      Lexp - Lyapunov exponents to each time value.
%
%    Users have to write their own ODE functions for their specified
%    systems and use handle of this function as rhs_ext_fcn - parameter.
%
%    Example. Lorenz system:
%               dx/dt = sigma*(y - x)     = f1
%               dy/dt = r*x - y - x*z = f2
%               dz/dt = x*y - b*z     = f3
%
%    The Jacobian of system:
%        | -sigma  sigma  0 |
%    J = |   r-z    -1   -x |
%        |    y      x   -b |
%
%    Then, the variational equation has a form:
%
%    F = J*Y
%    where Y is a square matrix with the same dimension as J.
%    Corresponding m-file:
%        function f=lorenz_ext(t,X)
%         SIGMA = 10; R = 28; BETA = 8/3;
%         x=X(1); y=X(2); z=X(3);
%
%         Y= [X(4), X(7), X(10);
%             X(5), X(8), X(11);
%             X(6), X(9), X(12)];
%         f=zeros(9,1);
%         f(1)=SIGMA*(y-x); f(2)=-x*z+R*x-y; f(3)=x*y-BETA*z;
%
%         Jac=[-SIGMA,SIGMA,0; R-z,-1,-x; y, x,-BETA];
%
%         f(4:12)=Jac*Y;
%
%    Run Lyapunov exponent calculation:
%
%    [T,Res]=lyapunov(3,@lorenz_ext,@ode45,0,0.5,200,[0 1 0],10);
%
%    See files: lorenz_ext, run_lyap.
%
% --------------------------------------------------------------------
% Copyright (C) 2004, Govorukhin V.N.
% This file is intended for use with MATLAB and was produced for MATDS-program
% http://www.math.rsu.ru/mexmat/kvm/matds/
% lyapunov.m is free software. lyapunov.m is distributed in the hope that it
% will be useful, but WITHOUT ANY WARRANTY.
%

%
%       n=number of nonlinear odes
%       n2=n*(n+1)=total number of odes
%

n1=n; n2=n1*(n1+1);

%  Number of steps

nit = round((tend-tstart)/stept);

% Memory allocation

y=zeros(n2,1); cum=zeros(n1,1); y0=y;
gsc=cum; znorm=cum;

% Initial values

y(1:n)=ystart(:);

for i=1:n1 y((n1+1)*i)=1.0; end;

t=tstart;

% Main loop

for ITERLYAP=1:nit

% Solutuion of extended ODE system

  [T,Y] = feval(fcn_integrator,rhs_ext_fcn,[t t+stept],y);

  t=t+stept;
  y=Y(size(Y,1),:);

  for i=1:n1
      for j=1:n1 y0(n1*i+j)=y(n1*j+i); end;
  end;

%
%       construct new orthonormal basis by gram-schmidt
%

  znorm(1)=0.0;
  for j=1:n1 znorm(1)=znorm(1)+y0(n1*j+1)^2; end;

  znorm(1)=sqrt(znorm(1));

  for j=1:n1 y0(n1*j+1)=y0(n1*j+1)/znorm(1); end;

  for j=2:n1
      for k=1:(j-1)
          gsc(k)=0.0;
          for l=1:n1 gsc(k)=gsc(k)+y0(n1*l+j)*y0(n1*l+k); end;
      end;

      for k=1:n1
          for l=1:(j-1)
              y0(n1*k+j)=y0(n1*k+j)-gsc(l)*y0(n1*k+l);
          end;
      end;

      znorm(j)=0.0;
      for k=1:n1 znorm(j)=znorm(j)+y0(n1*k+j)^2; end;
      znorm(j)=sqrt(znorm(j));

      for k=1:n1 y0(n1*k+j)=y0(n1*k+j)/znorm(j); end;
  end;

%
%       update running vector magnitudes
%

  for k=1:n1 cum(k)=cum(k)+log(znorm(k)); end;

%
%       normalize exponent
%

  for k=1:n1
      lp(k)=cum(k)/(t-tstart);
  end;

% Output modification

  if ITERLYAP==1
     Lexp=lp;
     Texp=t;
  else
     Lexp=[Lexp; lp];
     Texp=[Texp; t];
  end;

  if (mod(ITERLYAP,ioutp)==0)
     fprintf('t=%6.4f',t);
     for k=1:n1 fprintf(' %10.6f',lp(k)); end;
     fprintf('\n');
  end;

  for i=1:n1
      for j=1:n1
          y(n1*j+i)=y0(n1*i+j);
      end;
  end;

end;


text:
% EX_PLANE_STRAIN_RING: solve the plane-strain problem on one quarter of a cylinder.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_ring.txt';

% Type of boundary conditions
problem_data.nmnn_sides   = [2];
problem_data.drchlt_sides = [];
problem_data.press_sides  = [1];
problem_data.symm_sides   = [3 4];

% Physical parameters
E  =  1; nu = 0;
problem_data.lambda_lame = @(x, y) ((nu*E)/((1+nu)*(1-2*nu)) * ones (size (x)));
problem_data.mu_lame = @(x, y) (E/(2*(1+nu)) * ones (size (x)));

% Source and boundary terms
P = 1;
problem_data.f = @(x, y) zeros (2, size (x, 1), size (x, 2));
problem_data.g = @(x, y, ind) test_plane_strain_ring_g_nmnn (x, y, P, nu, ind);
problem_data.h = @(x, y, ind) test_plane_strain_ring_uex (x, y, E, nu, P);
problem_data.p = @(x, y, ind) P * ones (size (x));

% Exact solution (optional)
problem_data.uex = @(x, y) test_plane_strain_ring_uex (x, y, E, nu, P);

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree     = [3 3];     % Degree of the basis functions
method_data.regularity = [2 2];     % Regularity of the basis functions
method_data.nsub       = [9 9];     % Number of subdivisions
method_data.nquad      = [4 4];     % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space, u] = solve_linear_elasticity (problem_data, method_data);

% 4) POST-PROCESSING.
% 4.1) Export to Paraview
output_file = 'plane_strain_ring_Deg3_Reg2_Sub9';

vtk_pts = {linspace(0, 1, 21), linspace(0, 1, 21)};
fprintf ('results being saved in: %s \n \n', output_file)
sp_to_vtk (u, space, geometry, vtk_pts, output_file, {'displacement', 'stress'}, {'value', 'stress'}, ...
    problem_data.lambda_lame, problem_data.mu_lame)

% 4.2) Plot in Matlab. Comparison with the exact solution.
[eu, F] = sp_eval (u, space, geometry, vtk_pts);
[X, Y]  = deal (squeeze(F(1,:,:)), squeeze(F(2,:,:)));

subplot (1,2,1)
quiver (X, Y, squeeze(eu(1,:,:)), squeeze(eu(2,:,:)))
title ('Numerical solution'), axis equal tight
subplot (1,2,2)
eu2 = problem_data.uex (X, Y);
quiver (X, Y, squeeze(eu2(1,:,:)), squeeze(eu2(2,:,:)))
title ('Exact solution'), axis equal tight

error_l2 = sp_l2_error (space, msh, u, problem_data.uex)

%!demo
%! ex_plane_strain_ring

%!test
%! problem_data.geo_name = 'geo_ring.txt';
%! problem_data.nmnn_sides   = [2];
%! problem_data.drchlt_sides = [];
%! problem_data.press_sides  = [1];
%! problem_data.symm_sides   = [3 4];
%! E  =  1; nu = 0;
%! problem_data.lambda_lame = @(x, y) ((nu*E)/((1+nu)*(1-2*nu)) * ones (size (x)));
%! problem_data.mu_lame = @(x, y) (E/(2*(1+nu)) * ones (size (x)));
%! P = 1;
%! problem_data.f = @(x, y) zeros (2, size (x, 1), size (x, 2));
%! problem_data.g = @(x, y, ind) test_plane_strain_ring_g_nmnn (x, y, P, nu, ind);
%! problem_data.h = @(x, y, ind) test_plane_strain_ring_uex (x, y, E, nu, P);
%! problem_data.p = @(x, y, ind) P * ones (size (x));
%! problem_data.uex = @(x, y) test_plane_strain_ring_uex (x, y, E, nu, P);
%! method_data.degree     = [3 3];     % Degree of the basis functions
%! method_data.regularity = [2 2];     % Regularity of the basis functions
%! method_data.nsub       = [9 9];     % Number of subdivisions
%! method_data.nquad      = [4 4];     % Points for the Gaussian quadrature rule
%! [geometry, msh, space, u] = solve_linear_elasticity (problem_data, method_data);
%! error_l2 = sp_l2_error (space, msh, u, problem_data.uex);
%! assert (msh.nel, 81)
%! assert (space.ndof, 288)
%! assert (error_l2, 1.53640513551293e-06, 1e-16)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function info = rtt_buildinfodata_uno()
%RTT_BUILDINFODATA_UNO Defines build info data

% Copyright 2011 The MathWorks, Inc.

info.SourceFiles            = i_getSorceFiles;
info.SourceFilesToSkip      = {'main.cpp'};
info.IncludePaths           = i_getIncludePaths;
info.CompileFlags           = '-mmcu=atmega328p -std=gnu99 -Wall -Wstrict-prototypes -gstabs -Os';
info.Defines                = 'F_CPU=16000000';
info.LinkFlags              = '-lm -mmcu=atmega328p -gstabs';
end

% -------------------------------------------------------------------------
function files = i_getSorceFiles()
files = {};

% Get ARDUINO 0022 installation folder
hSetup = hwconnectinstaller.SetupInfo;
spPkg  = hSetup.getSpPkgInfo('Arduino');
arduinoRootDir = hSetup.getTpPkgRootDir('ARDUINO 0022', spPkg);

% Set paths
filesPath = fullfile(arduinoRootDir, 'hardware', 'arduino', 'cores', 'arduino');
cfiles = dir(fullfile(filesPath, '*.c'));
for i=1:length(cfiles)
    files{end+1}.Name = cfiles(i).name; %#ok<*AGROW>
    files{end}.Path = filesPath;
end
cfiles = dir(fullfile(filesPath, '*.cpp'));
for i=1:length(cfiles)
    files{end+1}.Name = cfiles(i).name;
    files{end}.Path = filesPath;
end
files{end+1}.Name = 'io_wrappers.cpp';
files{end}.Path = fullfile(spPkg.RootDir, 'src');
end

% -------------------------------------------------------------------------
function paths = i_getIncludePaths()
paths = {};

% Get ARDUINO 0022 installation folder
hSetup = hwconnectinstaller.SetupInfo;
spPkg  = hSetup.getSpPkgInfo('Arduino');
arduinoRootDir = hSetup.getTpPkgRootDir('ARDUINO 0022', spPkg);

% Set paths
paths{end+1} = fullfile(arduinoRootDir, 'hardware', 'arduino', 'cores', 'arduino');
paths{end+1} = fullfile(spPkg.RootDir, 'include');
end


text:
function p = prior_t(varargin)
%PRIOR_T  Student-t prior structure
%
%  Description
%    P = PRIOR_T('PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    creates Student's t-distribution prior structure in which the
%    named parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%    P = PRIOR_T(P, 'PARAM1', VALUE1, 'PARAM2', VALUE2, ...)
%    modify a prior structure with the named parameters altered
%    with the specified values.
%
%    The parameterization is as in Gelman, Carlin, Stern, Dunson, Vehtari,
%    and Rubin (2013). Bayesian Data Analysis, third edition.
%
%    Parameters for Student-t prior [default]
%      mu       - location [0]
%      s2       - scale [1]
%      nu       - degrees of freedom [4]
%      mu_prior - prior for mu [prior_fixed]
%      s2_prior - prior for s2 [prior_fixed]
%      nu_prior - prior for nu [prior_fixed]
%
%  See also
%    PRIOR_*
%
% Copyright (c) 2000-2001,2010 Aki Vehtari
% Copyright (c) 2009 Jarno Vanhatalo
% Copyright (c) 2010 Jaakko Riihimäki

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'PRIOR_T';
  ip.addOptional('p', [], @isstruct);
  ip.addParamValue('mu',0, @(x) isscalar(x));
  ip.addParamValue('mu_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('s2',1, @(x) isscalar(x) && x>0);
  ip.addParamValue('s2_prior',[], @(x) isstruct(x) || isempty(x));
  ip.addParamValue('nu',4, @(x) isscalar(x) && x>0);
  ip.addParamValue('nu_prior',[], @(x) isstruct(x) || isempty(x));
  ip.parse(varargin{:});
  p=ip.Results.p;

  if isempty(p)
    init=true;
    p.type = 't';
  else
    if ~isfield(p,'type') && ~isequal(p.type,'t')
      error('First argument does not seem to be a valid prior structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('mu',ip.UsingDefaults)
    p.mu = ip.Results.mu;
  end
  if init || ~ismember('s2',ip.UsingDefaults)
    p.s2 = ip.Results.s2;
  end
  if init || ~ismember('nu',ip.UsingDefaults)
    p.nu = ip.Results.nu;
  end
  % Initialize prior structure
  if init
    p.p=[];
  end
  if init || ~ismember('mu_prior',ip.UsingDefaults)
    p.p.mu=ip.Results.mu_prior;
  end
  if init || ~ismember('s2_prior',ip.UsingDefaults)
    p.p.s2=ip.Results.s2_prior;
  end
  if init || ~ismember('nu_prior',ip.UsingDefaults)
    p.p.nu=ip.Results.nu_prior;
  end

  if init
    % set functions
    p.fh.pak = @prior_t_pak;
    p.fh.unpak = @prior_t_unpak;
    p.fh.lp = @prior_t_lp;
    p.fh.lpg = @prior_t_lpg;
    p.fh.recappend = @prior_t_recappend;
  end

end

function [w, s, h] = prior_t_pak(p)
% This is a mandatory subfunction used for example
% in energy and gradient computations.

  w=[];
  s={};
  h=[];
  if ~isempty(p.p.mu)
    w = p.mu;
    s=[s; 't.mu'];
    h = 1;
  end
  if ~isempty(p.p.s2)
    w = [w log(p.s2)];
    s=[s; 'log(t.s2)'];
    h = [h 1];
  end
  if ~isempty(p.p.nu)
    w = [w log(p.nu)];
    s=[s; 'log(t.nu)'];
    h = [h 1];
  end
end

function [p, w] = prior_t_unpak(p, w)
% This is a mandatory subfunction used for example
% in energy and gradient computations.

  if ~isempty(p.p.mu)
    i1=1;
    p.mu = w(i1);
    w = w(i1+1:end);
  end
  if ~isempty(p.p.s2)
    i1=1;
    p.s2 = exp(w(i1));
    w = w(i1+1:end);
  end
  if ~isempty(p.p.nu)
    i1=1;
    p.nu = exp(w(i1));
    w = w(i1+1:end);
  end
end

function lp = prior_t_lp(x, p)
% This is a mandatory subfunction used for example
% in energy computations.

  lp=sum(gammaln((p.nu+1)./2) -gammaln(p.nu./2) -0.5*log(p.nu.*pi.*p.s2) -(p.nu+1)./2.*log(1+(x-p.mu).^2./p.nu./p.s2));

  if ~isempty(p.p.mu)
    lp = lp + p.p.mu.fh.lp(p.mu, p.p.mu);
  end
  if ~isempty(p.p.s2)
    lp = lp + p.p.s2.fh.lp(p.s2, p.p.s2) +log(p.s2);
  end
  if ~isempty(p.p.nu)
    lp = lp + p.p.nu.fh.lp(p.nu, p.p.nu) +log(p.nu);
  end
end

function lpg = prior_t_lpg(x, p)
% This is a mandatory subfunction used for example
% in gradient computations.

 %lpg=(p.nu+1)./p.nu .* (x-p.mu)./p.s2 ./ (1 + (x-p.mu).^2./p.nu./p.s2);
  lpg=-(p.nu+1).* (x-p.mu) ./ (p.nu.*p.s2 + (x-p.mu).^2);

  if ~isempty(p.p.mu)
    lpgmu = sum( (p.nu+1).* (x-p.mu) ./ (p.nu.*p.s2 + (x-p.mu).^2) ) + p.p.mu.fh.lpg(p.mu, p.p.mu);
    lpg = [lpg lpgmu];
  end
  if ~isempty(p.p.s2)
    lpgs2 = (sum( -1./(2.*p.s2) +((p.nu + 1).*(p.mu - x).^2)./(2.*p.s2.*((p.mu-x).^2 + p.nu.*p.s2))) + p.p.s2.fh.lpg(p.s2, p.p.s2)).*p.s2 + 1;
    lpg = [lpg lpgs2];
  end
  if ~isempty(p.p.nu)
    lpgnu = (0.5*sum( digamma1((p.nu+1)./2)-digamma1(p.nu./2)-1./p.nu-log(1+(x-p.mu).^2./p.nu./p.s2)+(p.nu+1)./(1+(x-p.mu).^2./p.nu./p.s2).*(x-p.mu).^2./p.s2./p.nu.^2) + p.p.nu.fh.lpg(p.nu, p.p.nu)).*p.nu + 1;
    lpg = [lpg lpgnu];
  end
end

function rec = prior_t_recappend(rec, ri, p)
% This subfunction is needed when using MCMC sampling (gp_mc).

% The parameters are not sampled in any case.
  rec = rec;
  if ~isempty(p.p.mu)
    rec.mu(ri,:) = p.mu;
  end
  if ~isempty(p.p.s2)
    rec.s2(ri,:) = p.s2;
  end
  if ~isempty(p.p.nu)
    rec.nu(ri,:) = p.nu;
  end
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function test_boundary()

grid = logical(randi(2,100,100)-1);

for j = 1:20
  ndx = randi(100^2);

  filled = imfill(~grid, ndx, 8);
  component = ~(xor(filled,grid));
  t0 = tic();
  boundary = imdilate(component, strel('square', 3)) - component;
  t_dil = toc(t0);
  fprintf('imdilate: %f,   ', t_dil);
  t0 = tic();
  fast_boundary = iris.terrain_grid.find_boundary(component);
  t_bound = toc(t0);
  fprintf('fast: %f,   ratio: %f\n', t_bound, t_dil/t_bound);

  figure(23)
  subplot 211
  imshow(boundary, 'InitialMagnification', 'fit');
  subplot 212
  imshow(fast_boundary, 'InitialMagnification', 'fit');
  assert(all(all(boundary == fast_boundary)));
end

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
% VL_TEST_INIT  Intialize test file
%   This is a VLFeat internal function used to initialize
%   a file containing tests.
%
%   See also: VL_TEST(), VL_HELP().

% Author:: Andrea Vedaldi

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% obtain path to caller m-file
stack = dbstack(1) ;
mfilePath = stack(1).file ;
mfileName = stack(1).name ;
fid = fopen(mfilePath, 'r') ;

% Parse source code to get subfunction names. It could use
% which('-subfun', stack(1).name) instead, but it is undocumented.

setupFunc = [] ;
teardownFunc = [] ;

tests = {} ;
l = 0 ;
t = 0 ;
while 1
  tline = fgetl(fid) ;
  l = l + 1 ;
  if ~ischar(tline), break ; end
  mo = regexp(tline, '^\s*function\s+(\w\s*=\s*)?test_(?<name>[\w_]+)\s*\(.*\).*$', ...
              'names') ;
  if length(mo) > 0
    t = t + 1 ;
    tests{t}.name = mo(1).name ;
    tests{t}.func = str2func(['test_' mo(1).name]) ;
    tests{t}.lineNumber = l ;
  end

  mo = regexp(tline, '^\s*function\s+(\w+\s*=\s*)?(?<name>setup|teardown)\s*\(.*\).*$', ...
              'names') ;
  if length(mo) > 0
    switch mo(1).name
      case 'setup'
        if isempty(setupFunc)
          setupFunc = @setup ;
        else
          warning('%s:%d: ignoring redundant definition of the setup function', ...
                  mfilePath) ;
        end
      case 'teardown'
        if isempty(teardownFunc)
          teardownFunc = @teardown ;
        else
          warning('%s:%d: ignoring redundant definition of the setup function', ...
                  mfilePath) ;
        end
      otherwise
        assert(0) ;
    end
  end
end

if ~isempty(setupFunc), state = setupFunc() ; end

if length(varargin) > 0
  testFunc = str2func(sprintf('test_%s', varargin{1})) ;
  if exist('state','var')
    testFunc(state) ;
  else
    testFunc() ;
  end
  if ~isempty(teardownFunc), teardownFunc(state) ; end
  return ;
end

% Run all tests collecting the results.
if length(tests) == 0
  error('No tests found in ''%s''', mfileName)  ;
end

for t = 1:length(tests)
  subTestName = sprintf('test_%s', tests{t}.name) ;
  results{t}.testName = sprintf('%s:%s', mfileName(9:end), tests{t}.name) ;
  results{t}.succeded = true ;
  results{t}.exception = [] ;
  try
    if exist('state', 'var')
      tests{t}.func(state) ;
    else
      tests{t}.func() ;
    end
  catch exception
    results{t}.succeded = false ;
    results{t}.exception = exception ;
  end
end

if ~isempty(teardownFunc), teardownFunc(state) ; end

results = cat(2, results{:}) ;

author: Cloud-CV

repo: object-proposals


text:
function [ wav_files ] = find_wav( path )
%FIND_WAV, find all wav file recursively
wav_files = [];
if(isdir(path) == 0)

    return;
end
path_files = dir(path);
fileNum = length(path_files);
for k= 3:fileNum
    file = [path,'\', path_files(k).name];
    if (path_files(k).isdir == 1)
        ret = find_wav(file);
        if(isempty(ret) ~= 1)
            if(isempty(wav_files))
                wav_files = char(ret);
            else
                wav_files = char(wav_files, ret);
            end
        end
    elseif strfind(path_files(k).name, '.wav')
        if(isempty(wav_files))
            wav_files = char(file);
        else
            wav_files = char(wav_files, file);
        end
    end
end
end

author: vipchengrui

repo: traditional-speech-enhancement


text:
%% This function will determine the noise level
% Code By: K
% Last Updated: 2019/07/11
%%
function noise_level=Determine_Noise_Level(percent_iter)
% Set the noise level
if percent_iter==1
    noise_level=0;
elseif percent_iter==2
    noise_level=1e-7;
elseif percent_iter==3
    noise_level=5e-7;
elseif percent_iter==4
    noise_level=1e-6;
elseif percent_iter==5
    noise_level=5e-6;
elseif percent_iter==6
    noise_level=1e-5;
elseif percent_iter==7
    noise_level=5e-5;
elseif percent_iter==8
    noise_level=1e-4;
elseif percent_iter==9
    noise_level=5e-4;
elseif percent_iter==10
    noise_level=1e-3;
elseif percent_iter<=14
    noise_level=2e-3*(percent_iter-10);
elseif percent_iter==15
    noise_level=1e-2;
elseif percent_iter<=19
    noise_level=2e-2*(percent_iter-15);
elseif percent_iter==20
    noise_level=1e-1;
elseif percent_iter<=24
    noise_level=1e-1+1e-1*(percent_iter-20);
end

author: dynamicslab

repo: SINDy-PI


text:
function [model, globalOpt] = svargplvmReOptimiseModel(model, initVardistIters, itNo, newExperimentNo, varargin)

% SVARGPLVMREOPTIMISEMODEL Optimise an already optimised svargplvm model for more iterations
%
% SEEALSO : svargplvmOptimiseModel
%
% COPYRIGHT : Andreas C. Damianou, 2011

% VARGPLVM

globalOptOrig = model.globalOpt;
globalOpt = globalOptOrig;

if nargin < 4
    newExperimentNo = globalOpt.experimentNo;
end
globalOpt.experimentNo = newExperimentNo;

globalOpt.initVardistIters = initVardistIters;
globalOpt.itNo = itNo;
model.globalOpt = globalOpt;

model = svargplvmOptimiseModel(model, varargin{:});

globalOptOrig.initVardistIters = [globalOptOrig.initVardistIters initVardistIters];
globalOptOrig.itNo = [globalOptOrig.itNo itNo];

globalOpt = globalOptOrig;
model.globalOpt = globalOpt;

author: SheffieldML

repo: vargplvm


text:
function [ y ] = back_trans( x, TransformSpecifics, reconDIM, n1, n2, n3 )
% used to set up transform handle

% M. Fischer, April 2016

%%

if ~iscell(x) && isreal(x) || iscell(x) && isreal(x{1}) % is this error safe?
    if strcmp(reconDIM,'2D')
        switch TransformSpecifics.TransformName
            case 'dwt'
                y = waverec2(x,TransformSpecifics.waveS,TransformSpecifics.FilterName);
            case 'cplxdt_matlab' % matlab variant:
                x_helper = TransformSpecifics.wt_struct;
%                 for k = 1:TransformSpecifics.Stages+1 % change to array to apply thresholding
%                     x_helper.cfs{k} = x{k}; % cfs{i} contains a 5D array.
%                 end;
                x_helper.cfs = x;
                y = idddtree2(x_helper);
            case 'cplxdt' % NOTE! n1 & n2 must be divisible by 2^J.
                for l = 1:TransformSpecifics.Stages % wavletStages (j)
                    for k = 1:2 % real & imag (i)
                        for j = 1:2 % orientation#1 (d1)
                            for h = 1:3 % orientation#2 (d2)
                                x_helper{1,l}{1,k}{1,j}{1,h} = x{l}(:,:,h,j,k); % contains a 5D array.
                            end;
                        end;
                    end;
                end;
                for j = 1:2 % orientation#1 (d1)
                    for h = 1:2 % orientation#2 (d2)
                        x_helper{1,TransformSpecifics.Stages+1}{1,j}{1,h} = x{TransformSpecifics.Stages+1}(:,:,h,j,1);
                    end;
                end;
                y = icplxdual2D(x_helper, TransformSpecifics.Stages, TransformSpecifics.Fsf, TransformSpecifics.sf);
                y = y(1:n1,1:n2);
            case 'shearlet'
                y = nsst_rec2(x.',TransformSpecifics.shear_f,TransformSpecifics.lpfilt);
                y = y(1:n1,1:n2);
            case 'shearlet_nontight'
                y = SLshearrec2D(x,TransformSpecifics.shearletSystem);
            otherwise
                warning('No correct Transform chosen');
        end;
    elseif strcmp(reconDIM,'3D')
        switch TransformSpecifics.TransformName
            case 'dwt'
                x_helper = TransformSpecifics.WDEC;
                x_helper.dec = x;
                y = waverec3(x_helper);
            case 'cplxdt'
                for l = 1:TransformSpecifics.Stages % wavletStages (j)
                    for k = 1:2 % (m)
                        for j = 1:2 % (n)
                            for h = 1:2 % (p)
                                for g = 1:7 % (d)
                                    x_helper{1,l}{1,k}{1,j}{1,h}{1,g} = x{l}(:,:,:,g,h,j,k); % contains a 7D array.
                                end;
                            end;
                        end;
                    end;
                end;
                for k = 1:2 % (m)
                    for j = 1:2 % (n)
                        for g = 1:2 % (d) not 1:7
                        x_helper{1,TransformSpecifics.Stages+1}{1,k}{1,j}{1,g} = x{TransformSpecifics.Stages+1}(:,:,:,g,1,j,k);
                        end;
                    end;
                end;
                y = icplxdual3D(x_helper, TransformSpecifics.Stages, TransformSpecifics.Fsf, TransformSpecifics.sf);
                y = y(1:n1,1:n2,1:n3);
            case 'shearlet'
                decomp_cells = 3*(4*4+2*2); % (fixed in get_transformSpecifics) dims:3, level:2 ,directions: 4x4 & 2x2
                % hard coded solution: (feel free to replace it)
                y = cell(1+decomp_cells,1);
                i = 1;
                for l = 1:3 % dims
                    k = 1; % (level)
                    for j = 1:4 % (direction)
                        for h = 1:4 % (direction)
                            x_helper.D{l,k}{j,h} = x{i};
                            i = i+1;
                        end;
                    end;
                    k =2;
                    for j = 1:2 % (direction)
                        for h = 1:2 % (direction)
                            x_helper.D{l,k}{j,h} = x{i};
                            i = i+1;
                        end;
                    end;
                end;
                x_helper.A = x{i+1};
                y = InvShearTransform3D(x_helper);
            case 'shearlet_nontight'
                if TransformSpecifics.flagExtend
                    x = img_extend(x,0,0,TransformSpecifics.extend_z);
                end;
                y = SLshearrec3D(x,TransformSpecifics.shearletSystem);
            otherwise
                warning('No correct Transform chosen');
        end;
    else
        error('ReconDim is not supported');
    end;
else % should always be real if operator_frame.m is used.
    error('use operator_frame.m');
end;

end

author: thomaskuestner


text:
function basis_brick8_test ( )

%*****************************************************************************80
%
%% BASIS_BRICK8_TEST verifies BASIS_BRICK8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    None
%
  node_num = 8;
  test_num = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'BASIS_BRICK8_TEST:\n' );
  fprintf ( 1, '  Verify basis functions for element BRICK8.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of nodes = %d\n', node_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The basis function values at basis nodes\n' );
  fprintf ( 1, '  should form the identity matrix.\n' );
  fprintf ( 1, '\n' );

  n = node_num;

  p = nodes_brick8 ( );

  phi = basis_brick8 ( n, p );

  for j = 1 : n
    fprintf ( 1, '  ' );
    for i = 1 : node_num
      fprintf ( 1, '%7.3f', phi(i,j) );
    end
    fprintf ( 1, '\n' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The basis function values at ANY point P\n' );
  fprintf ( 1, '  should sum to 1:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    ------------P-------------     PHI_SUM\n' );
  fprintf ( 1, '\n' );

  n = test_num;

  seed = 123456789;

  [ p, seed ] = r8mat_uniform_01 ( 3, n, seed );

  phi = basis_brick8 ( n, p );

  for j = 1 : n
    phi_sum = sum ( phi(1:node_num,j) );
    fprintf ( 1, '  %8.4f  %8.4f  %8.4f  %8.4f\n',  p(1:3,j), phi_sum );
  end

  return
end


text:
function y = acoshplus1(x)
%ACOSHPLUS1 Inverse hyperbolic cosine of argument plus one.
%
%   ACOSHPLUS1(X) is ACOSH(X+1) calculated in a way that is numerically better
%   when X is close to one.
%
%   This illustrates the difference
%
%      x = 2e-15*(0:0.01:1);
%      plot(x, acosh(x+1), 'b-', x, acoshplus1(x), 'r-');
%
%   See also COSHMINUS1.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-10-13 14:54:17 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   % check number of input arguments
   error(nargchk(1, 1, nargin));

   y = 2 * asinh(sqrt(x / 2));

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function mu = imCSMoment(img, p, q, varargin)
% Compute centered and scaled moments of an image.
%
%   MU = imCSMoment(IMG, P, Q)
%   compute (p+q)-th centered moment of image IMG.
%
%   MU = imCSMoment(IMG, P, Q, CENTER, MASS)
%   where CENTER = [CX CY], provides pre-computed center of mass.
%
%   Example
%   imCMoment
%
%   See also
%     imMoment, imCMoment, imEquivalentEllipse

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2008-10-08,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% get the centroid, either by argument, or by computation
if ~isempty(varargin)
    var = varargin{1};
    if length(var)>1
        cx = var(1);
        cy = var(2);
        varargin(1) = [];
    else
        cx = var;
        cy = varargin{2};
        varargin(1:2) = [];
    end
else
    s = sum(img(:));
    cx = imMoment(img, 1, 0)/s;
    cy = imMoment(img, 0, 1)/s;
end

if ~isempty(varargin)
    m00 = varargin{1};
else
    m00 = sum(img(:));
end

% image dimension
dim = size(img);
Dy = dim(1);
Dx = dim(2);

% compute x and y for each pixel
Ix = repmat((1:Dx), [Dy 1])-cx;
Iy = repmat((1:Dy)', [1 Dx])-cy;

% compute moment
mu = zeros(size(p));
for i = 1:length(p(:))
    d = (p(i)+q(i)) / 2 + 1;
    mu(i) = sum(Ix(:).^p(i) .* Iy(:).^q(i) .* img(:)) / m00^d;
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
%This directory contains files helping functions and utilities used for the
%visualization functions.
%
%CMAP_* - Colormaps for the visualization of scalp maps or score matrices
%DEFOPT_* - Default property lists for visualizing specific types of data
%
%AXIS_* - Helping functions for handling axes
%FIG_* - Helping functions for handling figures
%GRIDUTIL_* - Helping functions for grid_plot
%MNTUTIL_* - Helping functions for handling montages
%PLOTUTIL_* - Helping functions for plot_channel
%VISUTIL_* - Various helping functions

author: bbci


text:
function x = p03_x ( n )

%*****************************************************************************80
%
%% P03_X returns the least squares solution X for problem 3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Output, real X(N,1), the least squares solution.
%
  x = [ -7.5555556; 0.1111111; 7.7777778 ];

  return
end


text:
% test_cgal_fixed_alpha_shape3.m

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2013 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% tetrahedron

xyz = [
    0 0 0
    1 0 0
    0 1 0
    0 0 1
    ];

% compute alpha shape
tri = cgal_fixed_alpha_shape3(xyz, 1);

% plot meshes
subplot(1, 1, 1)
hold off
trisurf(tri{1}, xyz(:, 1), xyz(:, 2), xyz(:, 3))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% pyramid with an indentation in the base

xyz = [
    0 0 0
    1 0 0
    0 1 0
    .25 .25 0
    0 0 1
    ];

% compute alpha shape
tri = cgal_fixed_alpha_shape3(xyz, [0 0.5625    1.5469]);

% plot meshes
subplot(1, 2, 1)
hold off
trisurf(tri{2}, xyz(:, 1), xyz(:, 2), xyz(:, 3))
view(78, 40)
subplot(1, 2, 2)
hold off
trisurf(tri{3}, xyz(:, 1), xyz(:, 2), xyz(:, 3))
view(78, 40)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% crescent shape

xyz = [
    1.5 3 0
    3.5 3 0
    0 2 0
    1 2 0
    4 2 0
    5 2 0
    0 1 0
    1 1 0
    4 1 0
    5 1 0
    1 0 0
    4 0 0
    ...
    1.5 3 1
    3.5 3 1
    0 2 1
    1 2 1
    4 2 1
    5 2 1
    0 1 1
    1 1 1
    4 1 1
    5 1 1
    1 0 1
    4 0 1
    ];

% compute minimal alpha shape that creates one connected object
tri = cgal_fixed_alpha_shape3(xyz, 2.5157);
tri = tri{1};

% plot mesh
subplot(1, 1, 1)
trisurf(tri, xyz(:, 1), xyz(:, 2), xyz(:, 3))

% compute convex hull
tri = cgal_fixed_alpha_shape3(xyz, Inf);
tri = tri{1};

% plot mesh
subplot(1, 1, 1)
trisurf(tri, xyz(:, 1), xyz(:, 2), xyz(:, 3))

% two connected components
tri = cgal_fixed_alpha_shape3(xyz, 1.2656);
tri = tri{1};

% plot mesh
subplot(1, 1, 1)
trisurf(tri, xyz(:, 1), xyz(:, 2), xyz(:, 3))

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [x,aspc,spec] = invmelfcc(cep, sr, varargin)
% [x,aspc,spec] = invmelfcc(cep, sr[, opts ...])
%    Attempt to invert plp cepstra back to a full spectrum
%    and even a waveform.  Takes all the same options as melfcc.
%    x is (noise-excited) time domain waveform; aspc is the
%    auditory spectrogram, spec is the |STFT| spectrogram.
% 2005-05-15 dpwe@ee.columbia.edu

% Parse out the optional arguments
[wintime, hoptime, numcep, lifterexp, sumpower, preemph, dither, ...
 minfreq, maxfreq, nbands, bwidth, dcttype, fbtype, usecmp, modelorder, broaden, excitation] = ...
    process_options(varargin, 'wintime', 0.025, 'hoptime', 0.010, ...
          'numcep', 13, 'lifterexp', 0.6, 'sumpower', 1, 'preemph', 0.97, ...
	  'dither', 0, 'minfreq', 0, 'maxfreq', 4000, ...
	  'nbands', 40, 'bwidth', 1.0, 'dcttype', 2, ...
	  'fbtype', 'mel', 'usecmp', 0, 'modelorder', 0, 'broaden', ...
                    0, 'excitation', []);

winpts = round(wintime*sr);
nfft = 2^(ceil(log(winpts)/log(2)));

cep = lifter(cep, lifterexp, 1);   % 3rd arg nonzero means undo liftering

% Need to reconstruct the two extra flanking bands for invpostaud to delete
% (if we're doing usecmp)
pspc = cep2spec(cep, nbands+2*broaden, dcttype);

if (usecmp)
  aspc = invpostaud(pspc, maxfreq, fbtype, broaden);
else
  aspc = pspc;
end

% Undo the auditory spectrum
spec = invaudspec(aspc, sr, nfft, fbtype, minfreq, maxfreq, sumpower, bwidth);

% Back to waveform (modulate white noise, or specified excitation)
x = invpowspec(spec, sr, wintime, hoptime, excitation);

if preemph ~= 0
  % Undo the original preemphasis
  x = filter(1, [1 -preemph], x);
end

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
%% Mean-Shift Video Tracking
% by Sylvain Bernhardt
% July 2008
%% Description
% Estimate the density of data samples
% (here colour histogram) in a patch T
% with a kernel profile k. Lmap is the
% colormap length and H,W the patch size.

function q = Density_estim(T,Lmap,k,H,W,graph)
q = zeros(Lmap,1);
colour = linspace(1,Lmap,Lmap);
for x=1:W
    for y=1:H
        q(T(y,x)+1) = q(T(y,x)+1)+k(y,x);
    end
end

% Normalizing
C = 1/sum(sum(k));
q = C.*q;

% Plotting the estimated densities
if graph==1
    figure
    plot(colour,q);
end

author: Sable

repo: mcbench-benchmarks


text:
function [G11, G12, G21, G22] = d2Sbus_dV2(Ybus, V, lam, vcart)
%D2SBUS_DV2   Computes 2nd derivatives of power injection w.r.t. voltage.
%
%   The derivatives can be take with respect to polar or cartesian coordinates
%   of voltage, depending on the 4th argument.
%
%   [GAA, GAV, GVA, GVV] = D2SBUS_DV2(YBUS, V, LAM)
%   [GAA, GAV, GVA, GVV] = D2SBUS_DV2(YBUS, V, LAM, 0)
%
%   Returns 4 matrices containing the partial derivatives w.r.t. voltage angle
%   and magnitude of the product of a vector LAM with the 1st partial
%   derivatives of the complex bus power injections.
%
%   [GRR, GIR, GIR, GII] = D2SBUS_DV2(YBUS, V, LAM, 1)
%
%   Returns 4 matrices containing the partial derivatives w.r.t. real and
%   imaginary parts of voltage of the product of a vector LAM with the 1st
%   partial derivatives of the complex bus power injections.
%
%   Takes sparse bus admittance matrix YBUS, voltage vector V and nb x 1 vector
%   of multipliers LAM. Output matrices are sparse.
%
%   Examples:
%       [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
%       [Gaa, Gav, Gva, Gvv] = d2Sbus_dV2(Ybus, V, lam);
%
%       Here the output matrices correspond to:
%           Gaa = d/dVa (dSbus_dVa.' * lam)
%           Gav = d/dVm (dSbus_dVa.' * lam)
%           Gva = d/dVa (dSbus_dVm.' * lam)
%           Gvv = d/dVm (dSbus_dVm.' * lam)
%
%       [Grr, Gri, Gir, Gii] = d2Sbus_dV2(Ybus, V, lam, 1);
%
%       Here the output matrices correspond to:
%           Grr = d/dVr (dSbus_dVr.' * lam)
%           Gri = d/dVi (dSbus_dVr.' * lam)
%           Gir = d/dVr (dSbus_dVi.' * lam)
%           Gii = d/dVi (dSbus_dVi.' * lam)
%
%   For more details on the derivations behind the derivative code used
%   in MATPOWER, see:
%
%   [TN2]  R. D. Zimmerman, "AC Power Flows, Generalized OPF Costs and
%          their Derivatives using Complex Matrix Notation", MATPOWER
%          Technical Note 2, February 2010. [Online]. Available:
%          https://matpower.org/docs/TN2-OPF-Derivatives.pdf
%          doi: 10.5281/zenodo.3237866
%   [TN4]  B. Sereeter and R. D. Zimmerman, "AC Power Flows and their
%          Derivatives using Complex Matrix Notation and Cartesian
%          Coordinate Voltages," MATPOWER Technical Note 4, April 2018.
%          [Online]. Available: https://matpower.org/docs/TN4-OPF-Derivatives-Cartesian.pdf
%          doi: 10.5281/zenodo.3237909

%   MATPOWER
%   Copyright (c) 2008-2019, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%   and Baljinnyam Sereeter, Delft University of Technology
%
%   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.

%% default input args
if nargin < 4
    vcart = 0;      %% default to polar coordinates
end

n = length(V);

diaglam = sparse(1:n, 1:n, lam, n, n);
if vcart
    E = diaglam * conj(Ybus);
    F = E + E.';
    G = 1j * (E - E.');

    G11 = F;        %% Grr
    G21 = G;        %% Gir
    G12 = G21.';    %% Gri
    G22 = G11;      %% Gii
else
    Ibus    = Ybus * V;
    diagV   = sparse(1:n, 1:n, V, n, n);

    A = sparse(1:n, 1:n, lam .* V, n, n);
    B = Ybus * diagV;
    C = A * conj(B);
    D = Ybus' * diagV;
    E = conj(diagV) * (D * diaglam - sparse(1:n, 1:n, D*lam, n, n));
    F = C - A * sparse(1:n, 1:n, conj(Ibus), n, n);
    G = sparse(1:n, 1:n, ones(n, 1)./abs(V), n, n);

    G11 = E + F;
    G21 = 1j * G * (E - F);
    G12 = G21.';
    G22 = G * (C + C.') * G;
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function r=getRange(xCart,useHalfRange,zTx,zRx)
%%GETRANGE Obtain the bistatic (or monostatic) range measurements of
%          targets in the absence of refraction under non-relativistic
%          mechanics, ignoring atmospheric effects. The transmitter and the
%          target can be collocated.
%
%INPUTS: xCart The numDimXN set of N target positions. numDim is typically
%              2 or 3.
% useHalfRange A boolean value specifying whether the bistatic range value
%             should be divided by two. This normally comes up when
%             operating in monostatic mode, so that the range reported is
%             a one-way range. The default if this parameter is not
%             provided is false.
%         zTx A numDimXN matrix of the positions of the transmitters . If
%             this parameter is omitted or an empty matrix is passed, the
%             transmitters are assumed to be at the origin. If only a
%             single vector is passed, then the transmitter position is
%             assumed the same for all of the target states being converted.
%         zRx A numDimXN matrix of the positions of the receivers. If this
%             parameter is omitted or an empty matrix is passed, the
%             receivers are assumed to be at the origin. If only a single
%             vector is passed, then the receiver position is assumed the
%             same for all of the target states being converted.
%
%OUTPUTS: r The 1XN bistatic ranges of the targets.
%
%Bistatic range is discussed in [1].
%
%REFERENCES:
%[1] D. F. Crouse, "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%February 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

N=size(xCart,2);
numDim=size(xCart,1);

if(nargin<4||isempty(zRx))
    zRx=zeros(numDim,N);
elseif(size(zRx,2)==1)
    zRx=repmat(zRx,[1,N]);
end

if(nargin<3||isempty(zTx))
    zTx=zeros(numDim,N);
elseif(size(zTx,2)==1)
    zTx=repmat(zTx,[1,N]);
end

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

r=sqrt(sum((xCart-zTx(1:numDim,:)).^2,1))+sqrt(sum((xCart-zRx(1:numDim,:)).^2,1));

if(useHalfRange)
   r=r/2;
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function obj = compile(obj, export_path, varargin)
    % export the symbolic expressions of the system dynamics matrices and
    % vectors and compile as MEX files.
    %
    % Parameters:
    %  export_path: the path to export the file @type char
    %  varargin: variable input parameters @type varargin
    %   StackVariable: whether to stack variables into one @type logical
    %   File: the (full) file name of exported file @type char
    %   ForceExport: force the export @type logical
    %   BuildMex: flag whether to MEX the exported file @type logical
    %   Namespace: the namespace of the function @type char

    % export the impact constraints
    i_constrs = fieldnames(obj.ImpactConstraints);
    if ~isempty(i_constrs)
        for i=1:length(i_constrs)
            constr = i_constrs{i};
            export(obj.ImpactConstraints.(constr),export_path,varargin{:});

        end

    end

    % The mass matrix will be exported by the continuoud dynamics object,
    % but call it here again will not re-export if it has been exported.
    if ~isempty(obj.Mmat)
        cellfun(@(x)export(x,export_path,varargin{:}),obj.Mmat,'UniformOutput',false);
        %         export(obj.Mmat,export_path,varargin{:});
    end

    % call superclass method
    compile@DiscreteDynamics(obj, export_path, varargin{:});
end

author: ayonga

repo: frost-dev

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


text:
clc;
clear all;
for num=1:19
  I_ir=im2double(imread(strcat('',num2str(num-1),'.png')));
  [I,v1,v2]=rgb2ihs(I_ir);
  imwrite(I, strcat('',num2str(num),'.bmp'));
%   imwrite(v1,strcat('',num2str(num),'.bmp'));
%   imwrite(v2,strcat('',num2str(num),'.bmp'));

end

author: thfylsty


text:
% Demonstrates Picket Fence signal

% Initialization
clear all; close all; clc;

N  = 256;

x = SPX_SimpleSignals.picket_fence(N);

fx = fft(x);

spx.graphics.figure.full_screen;
subplot(211);
stem(x, '.');
title('Time domain');
subplot(212);
stem(fx, '.');
title('Frequency domain');

author: indigits

repo: sparse-plex