text:
% demonstrate usage of prior distributions
%
% See also priorDistributions.m.
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2015-07-13.
%                                      File automatically generated using noweb.
clear all, close all

% 1) specify some priors
% a) univariate priors
mu = 1.0; s2 = 0.01^2; nu = 3;
pg = {@priorGauss,mu,s2};                          % Gaussian prior
pl = {'priorLaplace',mu,s2};                        % Laplace prior
pt = {@priorT,mu,s2,nu};                        % Student's t prior
p1 = {@priorSmoothBox1,0,3,15};  % smooth box constraints lin decay
p2 = {@priorSmoothBox2,0,2,15};  % smooth box constraints qua decay
pd = {'priorDelta'}; % fix value of prior exclude from optimisation
pc = {@priorClamped};                         % equivalent to above
lam = 1.05; k = 2.5;
pw = {@priorWeibull,lam,k};                         % Weibull prior

% b) meta priors
pmx = {@priorMix,[0.5,0.5],{pg,pl}};        % mixture of two priors
g = @exp; dg = @exp; ig = @log;
ptr = {@priorTransform,g,dg,ig,pg};    % Gaussian in the exp domain

% c) multivariate priors
m = [1;2]; V = [2,1;1,2];
pG = {@priorGaussMulti,m,V};                    % 2d Gaussian prior
pD = {'priorDeltaMulti'};   % fix value of prior exclude from optim
pC = {@priorClampedMulti};                    % equivalent to above

% 2) evaluation
% pri = pt;   hp = randn(1,3);
% pri = pmx;  hp = randn(1,3);
% pri = ptr;  hp = randn(1,3);
pri = pG;   hp = randn(2,3);

% a) draw a sample from the prior
feval(pri{:})

% b) evaluate prior and derivative if requires
[lp,dlp] = feval(pri{:},hp)

% 3) comprehensive example
x = (0:0.1:10)'; y = 2*x+randn(size(x));   % generate training data
mean = {@meanSum,{@meanConst,@meanLinear}}; % specify mean function
cov = {@covSEiso}; lik = {@likGauss};  % specify covariance and lik
hyp.cov = [log(1);log(1.2)]; hyp.lik = log(0.9); hyp.mean = [2;3];
par = {mean,cov,lik,x,y}; mfun = @minimize; % input for GP function

% a) plain marginal likelihood optimisation (maximum likelihood)
im = @infExact;                                  % inference method
hyp_plain = feval(mfun, hyp, @gp, -10, im, par{:});      % optimise

% b) regularised optimisation (maximum a posteriori) with 1d priors
prior.mean = {pg;pc};  % Gaussian prior for first, clamp second par
prior.cov  = {p1;[]}; % box prior for first, nothing for second par
im = {@infPrior,@infExact,prior};                % inference method
hyp_p1 = feval(mfun, hyp, @gp, -10, im, par{:});         % optimise

% c) regularised optimisation (maximum a posteriori) with Nd priors
prior = [];                                   % clear the structure
% multivariate Student's t prior on the first and second mean hyper
prior.multi{1} = {@priorTMulti,[mu;mu],diag([s2,s2]),nu,...
                 struct('mean',[1,2])};          % use hyper struct
% Equivalent shortcut (same mu and s2 for all dimensions)
prior.multi{1} = {@priorTMulti,mu,s2,nu,struct('mean',[1,2])};
% multivariate Gaussian prior jointly on 1st and 3rd hyper
prior.multi{2} = {@priorGaussMulti,[mu;mu],diag([s2,s2]),...
                 [1,3]};               % use unwrapped hyper vector
% Equivalent shortcut (same mu and s2 for all dimensions)
prior.multi{2} = {@priorGaussMulti,mu,s2,[1,3]};
im = {@infPrior,@infExact,prior};                % inference method
hyp_pN = feval(mfun, hyp, @gp, -10, im, par{:});         % optimise

[unwrap2vec(hyp), unwrap2vec(hyp_plain), unwrap2vec(hyp_p1), unwrap2vec(hyp_pN)]

author: acerbilab

repo: bads

save_path: github-repos/MATLAB/acerbilab-bads


text:
classdef StateSpaceModelX < BaseX & dynamicprops
% StateSpaceModelX Abstract class
%
% Summary of StateSpaceModelX:
%  This is the base class for all TrackingX state-space models.
%  State-space model is effectivevely a container that brings together and
%  provides a common interface to all underlying models
%
% StateSpaceModelX Properties:
%   - Transition     Object handle to a TransitionModelX subclass
%   - Measurement    Object handle to an MeasurementModelX subclass
%   - Control        Object handle to a ControlModelX subclass
%   - Clutter        Object handle to a ClutterModelX subclass
%   - Detection      Object handle to an DetectionModelX subclass
%   - Birth          Object handle to a BirthModelX subclass
%
% StateSpaceModelX Methods:
%    StateSpaceModelX    - Constructor method
%
% See also TransitionModelX, MeasurementModelX and ControlModel template classes
%
% January 2018 Lyudmil Vladimirov, University of Liverpool.

    properties
        Transition
        Measurement
        Control
        Clutter
        Birth
        Detection
    end

    methods
        function this = StateSpaceModelX(varargin)
        % STATESPACEMODEL Constructor method
        %
        % Parameters
        % ==========
        % Transition: TransitionModelX subclass object handle
        %   A transition model
        % Measurement: MeasurementModelX subclass object handle
        %   A measurement model
        % Control: ControlModelX subclass object handle, optional
        %
        % Usage
        % =====
        % * StateSpaceModelX(ssm,Transition,Measurement) configures and returns the object
        %   ssm given the preconfigured transition and measurement models.
        % * StateSpaceModelX(config); where "config" is a structure MUST contain
        %   the required fields "config.Transition" and "config.Measurement" as objects
        %   handles, while it can optionally contain "config.Control"
        % * StateSpaceModelX(___,Name,Value) instantiates an object
        %   handle, configured with options specified by one or more Name,Value
        %   pair arguments.
        %
        %  See also TransitionModelX, MeasurementModelX and ControlModel.

            if(nargin==0)
                return;
            end

            % First check to see if a structure was received
            if(nargin==1)
                if(isstruct(varargin{1}))
                    if (isfield(varargin{1},'Transition'))
                        this.Transition = varargin{1}.Transition;
                    end
                    if (isfield(varargin{1},'Measurement'))
                        this.Measurement = varargin{1}.Measurement;
                    end
                    if (isfield(varargin{1},'Control'))
                        this.Control = varargin{1}.Control;
                    end
                    if (isfield(varargin{1},'Clutter'))
                        this.Clutter = varargin{1}.Clutter;
                    end
                    if (isfield(varargin{1},'Birth'))
                        this.Birth = varargin{1}.Birth;
                    end
                    if (isfield(varargin{1},'Detection'))
                        this.Detection = varargin{1}.Detection;
                    end
                end
                return;
            end

            % Otherwise, fall back to input parser
            parser = inputParser;
            parser.KeepUnmatched = true;
            parser.addOptional('Transition',[]);%, @(x) isa(x,'TransitionModelX'));
            parser.addOptional('Measurement',[]); %, @(x) isa(x,'MeasurementModelX'));
            parser.addParameter('Control',[]);
            parser.addParameter('Clutter',[]);
            parser.addParameter('Birth',[]);
            parser.addParameter('Detection',[]);
            parser.parse(varargin{:});

            this.Transition = parser.Results.Transition;
            this.Measurement = parser.Results.Measurement;
            this.Clutter = parser.Results.Clutter;
            this.Birth = parser.Results.Birth;
            this.Detection = parser.Results.Detection;
        end
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function M = callback_blurring( M, dir, options )

% callback_blurring - callback for the deconvolution with wavelets
%
%   M = callback_blurring( M, dir, options );
%
%   Copyright (c) 2007 Gabriel Peyre

if isfield(options, 'eta')
    eta = options.eta;
else
    eta = 5;
end
if isfield(options, 'eta')
    Jmin = options.Jmin;
else
    Jmin = 4;
end

n = sqrt(size(M,1));
M = reshape(M,n,n);

if dir==1
    % wavelet synthesis and then blurring
    M = perform_wavelet_transform(M, Jmin, -1, options );
    M = perform_blurring( M, eta );
else
    % transpose
    M = perform_blurring( M, eta );
    M = perform_wavelet_transform(M, Jmin, 1, options );
end

M = M(:);

author: gpeyre

repo: matlab-toolboxes


text:
function a = schur_block_inverse ( n, x, y )

%*****************************************************************************80
%
%% SCHUR_BLOCK_INVERSE returns the inverse of the SCHUR_BLOCK matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real X( (N+1)/2 ), specifies the diagonal elements
%    of A.
%
%    Input, real Y( N/2 ), specifies the off-diagonal elements
%    of the Schur blocks.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for i = 1 : n
    for j = 1 : n

      k = floor ( ( i + 1 ) / 2 );

      if ( i == j )

        if ( i == n & mod ( n, 2 ) == 1 )
          a(i,j) = 1.0 / x(k);
        else
          a(i,j) = x(k) / ( x(k)^2 + y(k)^2 );
        end

      elseif ( mod ( i, 2 ) == 1 & j == i + 1 )

        a(i,j) = - y(k) / ( x(k)^2 + y(k)^2 );

      elseif ( mod ( i, 2 ) == 0 & j == i - 1 )

        a(i,j) =   y(k) / ( x(k)^2 + y(k)^2 );

      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Signal Processing Q17734
% https://dsp.stackexchange.com/questions/17734
% Estimate the Discrete Fourier Series of a Signal with Missing Samples
% References:
%   1.  aa
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.001     17/01/2021  Royi Avital     RoyiAvital@yahoo.com
%   *   First release.
% - 1.0.000     20/07/2017
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;
run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = OFF;

%% Simulation Parameters

numSamples      = 150;
numGivenSamples = 120; %<! Must be not greater than 'numSamples'
filterNumCoeff  = 15;

%% Generate Data and Samples

vFilterCoeff = rand([filterNumCoeff, 1]);
vFilterCoeff = vFilterCoeff / sum(vFilterCoeff);

% Generating the signal
vX = randn([numSamples, 1]);
vX = conv2(vX, vFilterCoeff, 'same');

% Generating the given data
vGivenIdx       = sort(randperm(numSamples, numGivenSamples));
vXX             = vX(vGivenIdx);

% Generating the DFT Matrix with the respected rows
% Pay attention to scaling in order to make it Unitary
mF = dftmtx(numSamples) / sqrt(numSamples);
mA = mF(:, vGivenIdx);

%% Estimation

% This is equivalent to padding with zeros all the samples which are not
% given.

% vY = (mA * mA') \ (mA * vXX);
vY = pinv(mA') * vXX;

%% Display Results

refDftIdx   = floor((numSamples + 1) / 2);
vDftIdx     = [0:(refDftIdx - 1)];

vXDftEstimated  = vY(1:refDftIdx);
vXDft           = mF * vX;
vXDft           = vXDft(1:refDftIdx);

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
hLineSeries = plot(vDftIdx, abs([vXDft, vXDftEstimated]));
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', ['Estimating the DFT Given Missing Samples'], ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', 'DFT Sample Index [n]', ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', 'Amplitude', ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['The DFT of the Complete Signal'], ['Estimated DFT from Partial Signal']});

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
% Remove bad tracking coordinates (position jumps)
%
% A position "jump" correspond to position samples that imply that the rat is
% moving quicker than physical possible.
%
%
function [x, y] = removePosJumps(x, y, threshold, stdThreshold)
    N = length(x);
    % Indexes to position samples that are to be removed
    remInd = zeros(N,1);
    remCounter = 0;

    diffX = diff(x);
    diffY = diff(y);
    diffR = sqrt(diffX.^2 + diffY.^2);
    ind = find(diffR > threshold);

    if isempty(ind)
        return;
    end

    if ind(end) == length(x)
        offset = 2;
    else
        offset = 1;
    end

    for ii = 1:length(ind)-offset
        if ind(ii+1) == ind(ii)+1
            % A single sample position jump, tracker jumps out one sample and
            % then jumps back to path on the next sample. Remove bad sample.
            remCounter = remCounter + 1;
            remInd(remCounter) = ind(ii)+1;
            ii = ii+1;
            continue
        else
            % Not a single jump. 2 possibilities:
            % 1. Tracker jumps out, and stay out at the same place for several
            % samples and then jumps back.
            % 2. Tracker just has a small jump before path continues as normal,
            % unknown reason for this. In latter case the samples are left
            % untouched.
            idx = find(x(ind(ii)+1:ind(ii+1)+1)==x(ind(ii)+1));
            if length(idx) == length(x(ind(ii)+1:ind(ii+1)+1));
                n = ind(ii+1)+1 - ind(ii);
                remInd(remCounter+1:remCounter+n) = (ind(ii)+1:ind(ii+1)+1)';
                remCounter = remCounter + n;
            end
        end
    end

    remInd = remInd(1:remCounter);

    % Remove the samples
    x(remInd) = NaN;
    y(remInd) = NaN;

    % there could be tracking outliers. They are commonly several values between longer list of
    % NaNs. These values lie far away from 'good' points, so discard everything that is
    % further than 2.5*STD
    med(1) = nanmin(x) + (nanmax(x) - nanmin(x))/2; % roughy middle point of the arena, better
    med(2) = nanmin(y) + (nanmax(y) - nanmin(y))/2; % than median or mean since it's less biased
    std = nanstd([x, y]);
    lowerBound = med - (stdThreshold*std);
    upperBound = med + (stdThreshold*std);

    x(x < lowerBound(1)) = nan;
    x(x > upperBound(1)) = nan;

    y(y < lowerBound(2)) = nan;
    y(y > upperBound(2)) = nan;
end

author: kavli-ntnu


text:
function [ value, ifault ] = xinbta ( p, q, beta, alpha )

%*****************************************************************************80
%
%% XINBTA computes inverse of the incomplete Beta function.
%
%  Discussion:
%
%    The accuracy exponent SAE was loosened from -37 to -30, because
%    the code would not otherwise accept the results of an iteration
%    with p = 0.3, q = 3.0, alpha = 0.2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2014
%
%  Author:
%
%    Original FORTRAN77 version by GW Cran, KJ Martin, GE Thomas.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    GW Cran, KJ Martin, GE Thomas,
%    Remark AS R19 and Algorithm AS 109:
%    A Remark on Algorithms AS 63: The Incomplete Beta Integral
%    and AS 64: Inverse of the Incomplete Beta Integeral,
%    Applied Statistics,
%    Volume 26, Number 1, 1977, pages 111-114.
%
%  Parameters:
%
%    Input, real P, Q, the parameters of the incomplete
%    Beta function.
%
%    Input, real BETA, the logarithm of the value of
%    the complete Beta function.
%
%    Input, real ALPHA, the value of the incomplete Beta
%    function.  0 <= ALPHA <= 1.
%
%    Output, real VALUE, the argument of the incomplete
%    Beta function which produces the value ALPHA.
%
%    Output, integer IFAULT, error flag.
%    0, no error occurred.
%    nonzero, an error occurred.
%
%  Local Parameters:
%
%    Local, real SAE, requests an accuracy of about 10^SAE.
%
  sae = -30.0;

  fpu = 10.0 ^ sae;

  ifault = 0;
  value = alpha;
%
%  Test for admissibility of parameters.
%
  if ( p <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'XINBTA - Fatal error!\n' );
    fprintf ( 1, '  P <= 0.0\n' );
    ifault = 1;
    error ( 'XINBTA - Fatal error!' );
  end

  if ( q <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'XINBTA - Fatal error!\n' );
    fprintf ( 1, '  Q <= 0.0\n' );
    ifault = 1;
    error ( 'XINBTA - Fatal error!' );
  end

  if ( alpha < 0.0 || 1.0 < alpha )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'XINBTA - Fatal error!\n' );
    fprintf ( 1, '  ALPHA not between 0 and 1.\n' );
    ifault = 2;
    error ( 'XINBTA - Fatal error!' );
  end
%
%  If the answer is easy to determine, return immediately.
%
  if ( alpha == 0.0 )
    value = 0.0;
    return
  end

  if ( alpha == 1.0 )
    value = 1.0;
    return
  end
%
%  Change tail if necessary.
%
  if ( 0.5 < alpha )
    a = 1.0 - alpha;
    pp = q;
    qq = p;
    indx = 1;
  else
    a = alpha;
    pp = p;
    qq = q;
    indx = 0;
  end
%
%  Calculate the initial approximation.
%
  r = sqrt ( - log ( a * a ) );

  y = r - ( 2.30753 + 0.27061 * r ) ...
    / ( 1.0 + ( 0.99229 + 0.04481 * r ) * r );

  if ( 1.0 < pp && 1.0 < qq )

    r = ( y * y - 3.0 ) / 6.0;
    s = 1.0 / ( pp + pp - 1.0 );
    t = 1.0 / ( qq + qq - 1.0 );
    h = 2.0 / ( s + t );
    w = y * sqrt ( h + r ) / h - ( t - s ) ...
      * ( r + 5.0 / 6.0 - 2.0 / ( 3.0 * h ) );
    value = pp / ( pp + qq * exp ( w + w ) );

  else

    r = qq + qq;
    t = 1.0 / ( 9.0 * qq );
    t = r * ( 1.0 - t + y * sqrt ( t ) )^3;

    if ( t <= 0.0 )
      value = 1.0 - exp ( ( log ( ( 1.0 - a ) * qq ) + beta ) / qq );
    else

      t = ( 4.0 * pp + r - 2.0 ) / t;

      if ( t <= 1.0 )
        value = exp ( ( log ( a * pp ) + beta ) / pp );
      else
        value = 1.0 - 2.0 / ( t + 1.0 );
      end

    end

  end
%
%  Solve for X by a modified Newton-Raphson method,
%  using the function BETAIN.
%
  r = 1.0 - pp;
  t = 1.0 - qq;
  yprev = 0.0;
  sq = 1.0;
  prev = 1.0;

  if ( value < 0.0001 )
    value = 0.0001;
  end

  if ( 0.9999 < value )
    value = 0.9999;
  end

  iex = max ( - 5.0 / pp / pp - 1.0 / a ^ 0.2 - 13.0, sae );

  acu = 10.0 ^ iex;
%
%  Iteration loop.
%
  while ( 1 )

    [ y, ifault ] = betain ( value, pp, qq, beta );

    if ( ifault ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'XINBTA - Fatal error!\n' );
      fprintf ( 1, '  BETAIN returns IFAULT = %d\n', ifault );
      error ( 'XINBTA - Fatal error!' );
      ifault = 3;
      return
    end

    xin = value;
    y = ( y - a ) * exp ( beta + r * log ( xin ) + t * log ( 1.0 - xin ) );

    if ( y * yprev <= 0.0 )
      prev = max ( sq, fpu );
    end

    g = 1.0;

    while ( 1 )

      while ( 1 )

        adj = g * y;
        sq = adj * adj;

        if ( sq < prev )

          tx = value - adj;

          if ( 0.0 <= tx && tx <= 1.0 )
            break
          end

        end

        g = g / 3.0;

      end
%
%  Check whether current estimate is acceptable.
%  The change "VALUE = TX" was suggested by Ivan Ukhov.
%
      if ( prev <= acu && y * y <= acu )
        value = tx;
        if ( indx )
          value = 1.0 - value;
        end
        return
      end

      if ( tx ~= 0.0 && tx ~= 1.0 )
        break
      end

      g = g / 3.0;

    end

    if ( tx == value )
      break
    end

    value = tx;
    yprev = y;

  end

  if ( indx )
    value = 1.0 - value;
  end

  return
end


text:
function changed = Changed(Problem,Population)
% Detect whether the problem changes

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    RePop1  = Population(randperm(end,ceil(end/10)));
    RePop2  = Problem.Evaluation(RePop1.decs);
    changed = ~isequal(RePop1.objs,RePop2.objs) || ~isequal(RePop1.cons,RePop2.cons);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function determ = skew_circulant_determinant ( n, x )

%*****************************************************************************80
%
%% SKEW_CIRCULANT_DETERMINANT returns the determinant of the SKEW_CIRCULANT matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Input, real X(N), the values in the first row of A.
%
%    Output, real DETERM, the determinant.
%
  determ = 1.0;

  j_hi = floor ( ( n + 1 ) / 2 );

  for j = 1 : j_hi

    lambda = 0.0;

    for k = 1 : n
      angle = ( 2 * j - 1 ) * ( k - 1 ) * pi / n;
      lambda = lambda + x(k) * complex ( cos ( angle ), sin ( angle ) );
    end

    if ( 2 * j <= n )
      determ = determ * ( abs ( lambda ) )^2;
    else
      determ = determ * real ( lambda );
    end

  end

  return
end


text:
function r = sqr(a)
%SQR          Hessian (elementwise) square
%

% written  04/04/04     S.M. Rump
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

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

  K = prod(size(a.x));
  if K==1                   % scalar hessian

    r.x = sqr(a.x);
    r.dx = (2*a.x) * a.dx;
    r.hx = a.hx * (2*a.x) + reshape(a.dx * a.dx.',size(a.hx));

  else                      % matrix hessian

    N = getappdata(0,'INTLAB_HESSIAN_NUMVAR');
    N2 = N^2;

    r.x = sqr(a.x);
    if issparse(a.hx)               % input sparse

      ax = 2*a.x(:);
      sizeax = length(ax);
      [ia,ja,sa] = find(a.dx);
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if isempty(ia)
        r.dx = sparse([],[],[],N,sizeax);
        r.hx = sparse([],[],[],N2,sizeax);
      else
        if isa(a.x,'intval')          % sparse intval
          rdx = times(ax(ja),sa(:),0);
          if rdx.complex
            r.dx = intval( sparse(ia,ja,rdx.mid,N,sizeax) , sparse(ia,ja,rdx.rad,N,sizeax) , 'midrad' );
          else
            r.dx = intval( sparse(ia,ja,rdx.inf,N,sizeax) , sparse(ia,ja,rdx.sup,N,sizeax) , 'infsup' );
          end
        else                          % sparse point
          r.dx = sparse(ia,ja,ax(ja).*sa(:),N,sizeax);
        end
        r.hx = adx2rhx(N,sizeax,a.dx);
      end
      [ia,ja,sa] = find(a.hx);        % sparse point or intval
      % check for emptyness: cures Matlab bug
      % a=sparse([],[],[],2,1), [i,j,s]=find(a), s(i).*s(:)  yields error
      if ~isempty(ia)
        if isa(a.x,'intval')
          rhx = times(ax(ja),sa(:),0);
          if rhx.complex
            r.hx = r.hx + intval( sparse(ia,ja,rhx.mid,N2,sizeax) , sparse(ia,ja,rhx.rad,N2,sizeax) , 'midrad' );
          else
            r.hx = r.hx + intval( sparse(ia,ja,rhx.inf,N2,sizeax) , sparse(ia,ja,rhx.sup,N2,sizeax) , 'infsup' );
          end
        else
          r.hx = r.hx + sparse(ia,ja,ax(ja).*sa(:),N2,sizeax);
        end
      end

    else                            % input full

      ax = 2*a.x(:).';
      ax = ax(ones(N2,1),:);
      r.dx = a.dx .* ax(1:N,:);
      r.hx = a.hx .* ax + a.dx(repmat(1:N,N,1),:) .* a.dx(repmat(1:N,1,N),:);

    end

  end

  r = class(r,'hessian');

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
% SERVER Write a message over the specified port
%
% Usage - server(message, output_port, number_of_retries)
function server(message, output_port, number_of_retries)

    import java.net.ServerSocket
    import java.io.*

    if (nargin < 3)
        number_of_retries = 20; % set to -1 for infinite
    end
    retry             = 0;

    server_socket  = [];
    output_socket  = [];

    while true

        retry = retry + 1;

        try
            if ((number_of_retries > 0) && (retry > number_of_retries))
                fprintf(1, 'Too many retries\n');
                break;
            end

            fprintf(1, ['Try %d waiting for client to connect to this ' ...
                        'host on port : %d\n'], retry, output_port);

            % wait for 1 second for client to connect server socket
            server_socket = ServerSocket(output_port);
            server_socket.setSoTimeout(1000);

            output_socket = server_socket.accept;

            fprintf(1, 'Client connected\n');

            output_stream   = output_socket.getOutputStream;
            d_output_stream = DataOutputStream(output_stream);

            % output the data over the DataOutputStream
            % Convert to stream of bytes
            fprintf(1, 'Writing %d bytes\n', length(message))
            d_output_stream.writeBytes(char(message));
            d_output_stream.flush;

            % clean up
            server_socket.close;
            output_socket.close;
            break;

        catch
            if ~isempty(server_socket)
                server_socket.close
            end

            if ~isempty(output_socket)
                output_socket.close
            end

            % pause before retrying
            pause(1);
        end
    end
end

author: Sable

repo: mcbench-benchmarks


text:
function [H,p,ci,stats] = ttest2_printout(sc1,sc2,varargin)
% :Usage:
% ::
%
%     [H,p,ci,stats] = ttest2_printout(sc1,sc2,[doplot],[covts])
%
% one or two sample t-test printout and plot
% covariates are not done yet!
%
% :Inputs:
%
%   **sc1:**
%        data from first group
%
%   **sc2:**
%        data from second group (if missing or empty, performs one-sample
%        t-test)
%
% ..
%    tor wager, last updated Sept 2007 (cosmetic update)
% ..

if length(varargin) > 2
    % covariates of no interest
    X = varargin{2};
    % NOT DONE YET.
end

if nargin == 1 || isempty(sc2)
    [H,p,ci,stats] = ttest(sc1,0,.05,'both');
    fprintf(1,'u1 = %3.2f, t(%3.1f) = %3.2f, p = %3.4f\n',nanmean(sc1),stats.df,stats.tstat,p);
else
    [H,p,ci,stats] = ttest2(sc1,sc2,.05,'both','unequal');
    fprintf(1,'u1 = %3.2f, u2 = %3.2f, udiff = %3.2f, t(%3.1f) = %3.2f, p = %3.4f\n',nanmean(sc1),nanmean(sc2),nanmean(sc1) - nanmean(sc2),stats.df,stats.tstat, p);
end

means = [nanmean(sc1) nanmean(sc2)];
stats.means = means;

if length(varargin) > 0 && varargin{1}

    hh = bar(means);
    set(hh,'FaceColor',[.7 .7 .7]);

    se = stats.sd ./ sqrt(length(sc1)+length(sc2));

    se = se';
    tor_bar_steplot(means,se,{'k'});

    set(gca,'XTick',[1 2],'XTickLabel',{'High' 'Low'});

end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function p = normal2(mu, Sigma, dom)
%NORMAL2   Bivariate normal distribution.
%   NORMAL2(MU, SIGMA, DOMAIN) returns the joint probability density for two
%   variables whose means are in the vector MU and whose covariance is the 2x2
%   positive definite matrix SIGMA. The result is a chebfun2 defined on the
%   (finite) domain DOMAIN.
%
%   If DOMAIN is not supplied, the domain is taken to be large enough to make
%   the density essentially zero at the boundary.
%
%   Example:
%      p = cheb.normal2([1,0], [1 1;1 5]);
%      surf(p)

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

if ( any(eig(Sigma) < 0) || norm(Sigma - Sigma') > 0 )
    error('CHEB:NORMAL2:covariance:nonSymPosDef', ...
        'Covariance matrix must be symmetric positive definite.')
end

if nargin < 3
    sig = svd(Sigma);
    dom = [ mu(1)+4*[-sig(1) sig(1)], mu(2)+4*[-sig(1) sig(1)] ];
end

x = chebfun2(@(x,y) x,dom);
y = chebfun2(@(x,y) y,dom);
invSig = inv(Sigma);
x = x - mu(1);
y = y - mu(2);
z = invSig(1,1)*x.^2 + (invSig(1,2)+invSig(2,1))*x.*y + invSig(2,2)*y.^2;
const = 2*pi*sqrt(det(Sigma));

p = exp(-z/2)/const;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function outfig=view(Tree,varargin)
%VIEW View classification or regression tree graphically.
%   VIEW(T) displays the decision tree T as computed by the CLASSREGTREE
%   constructor in a figure window.  Each branch in the tree is labeled
%   with its decision rule, and each terminal node is labeled with the
%   predicted value for that node.  Click on any node to get more
%   information about it.  The information displayed is specified by
%   leftmost pop-up menu at the top of the figure.
%
%   VIEW(T,'PARAM1',val1,'PARAM2',val2,...) specifies optional
%   parameter name/value pairs:
%
%      'names'       A cell array of names for the predictor variables,
%                    in the order in which they appear in the X matrix
%                    from which the tree was created (see TREEFIT)
%      'prunelevel'  Initial pruning level to display
%
%   For each branch node, the left child node corresponds to the points
%   that satisfy the condition, and the right child node corresponds to
%   the points that do not satisfy the condition.
%
%   Example:  Create and graph classification tree for Fisher's iris data.
%             The names are abbreviations for the column contents (sepal
%             length, sepal width, petal length, and petal width).
%      load fisheriris;
%      t = classregtree(meas, species);
%      view(t,'names',{'SL' 'SW' 'PL' 'PW'});
%
%   See also CLASSREGTREE, CLASSREGTREE/EVAL, CLASSREGTREE/TEST, CLASSREGTREE/PRUNE.

%   Copyright 2006-2009 The MathWorks, Inc.
%   $Revision: 1.1.6.7 $  $Date: 2009/05/07 18:32:53 $

% Process inputs
nvars = max(abs(Tree.var));
names = Tree.names;
okargs =   {'names' 'prunelevel'};
defaults = {names   0};
[eid,emsg,names,curlevel] = getargs(okargs,defaults,varargin{:});
if ~isempty(eid)
   error(sprintf('stats:treedisp:%s',eid),emsg);
end
if isempty(names)
   names = cell(nvars,1);
   for j=1:nvars
      names{j} = sprintf('x%d',j);
   end
elseif length(names)<nvars
   error('stats:treedisp:InvalidNames',...
         'Incorrect number of variable names.');
end

doclass = isequal(Tree.method,'classification');

% Create empty figure and axes to receive tree display
fig = setupfigure(doclass);
try
   adjustmenu(fig);
catch
end

% Draw tree
[X,Y,Tree] = drawtree(Tree,fig,names,curlevel);

% Save information for call-backs
fulltree = Tree;
set(fig,'ButtonDownFcn',@removelabels, ...
        'UserData',{X Y 0 names doclass fulltree curlevel});

% Update ui elements
updateenable(fig);
updatelevel(fig,curlevel,fulltree);

if nargout>0
    outfig = fig;
end

% ----------------------------------------------
function [X,Y,Tree] = drawtree(Tree,fig,names,curlevel)
%DRAWTREE Draw decision tree into specified figure.

ax = get(fig,'CurrentAxes');
splitvar = Tree.var;
cutoff = Tree.cut;
parentnode = Tree.parent;
nonroot = parentnode~=0;

% Get coordinates of nodes
isleaf = splitvar==0;
[X,Y] = layouttree(Tree,isleaf);

% Which leaves and branches are pruned?
if ~isempty(Tree.prunelist)
   isbranch = (Tree.prunelist>curlevel);
else
   isbranch = ~isleaf;
end
if any(isbranch)
   isleaf = false(size(isbranch));
   c = Tree.children(isbranch,:);
   c = c(c>0);
   isleaf(c) = 1;
   isleaf = isleaf & ~isbranch;
else
   isleaf = ~nonroot;
end
pruned = ~(isleaf | isbranch);

branchnodes = find(isbranch);
leafnodes = find(isleaf);
doclass = isequal(Tree.method,'classification');

% Get coordinates of connecting lines
p = parentnode(nonroot & ~pruned);
x = [X(nonroot & ~pruned)'; X(p)'; NaN+p'];
y = [Y(nonroot & ~pruned)'; Y(p)'; NaN+p'];

% Plot nodes and connections for unpruned nodes, but stop listening to axis
% and remember some things that may get changed during the plotting
axislistener(ax,false);
xlim = get(ax,'XLim');
ylim = get(ax,'YLim');
ud = get(ax,'UserData');
h = plot(X(branchnodes),Y(branchnodes),'b^', ...
         X(leafnodes),Y(leafnodes),'b.', ...
         x(:),y(:),'b-','Parent',ax);
set(ax,'UserData',ud,'Visible','off','XLim',xlim,'YLim',ylim);
axislistener(ax,true);

% Same for pruned nodes
t = nonroot & pruned;
p = parentnode(t);
x = [X(t)'; X(p)'; NaN+p'];
y = [Y(t)'; Y(p)'; NaN+p'];
line('Parent',ax,'XData',X(pruned),'YData',Y(pruned),'Tag','prunednode',...
     'Marker','o','Color',[.2 .2 .2],'Linestyle','none','HitTest','off');
line('Parent',ax,'XData',x(:),'YData',y(:),'Tag','prunedconnection',...
     'Marker','none','LineStyle',':','Color',[.2 .2 .2],'HitTest','off');
if length(h)==3
   set(h(1),'ButtonDownFcn',@labelpoint,'Tag','branch','MarkerSize',10);
   set(h(2),'ButtonDownFcn',@labelpoint,'Tag','leaf','MarkerSize',20);
   set(h(end),'HitTest','off','Tag','connection');
else
   set(h,'ButtonDownFcn',@labelpoint,'Tag','leaf','MarkerSize',20);
end

% Label leaf nodes with class, branch nodes with split rule
if doclass
   cnames = Tree.classname;
   ctext = cnames(Tree.class(leafnodes));
else
   ctext = num2str(Tree.class(leafnodes));
end

h = findobj(fig,'Tag','menuleaf');
vis = get(h,'Checked');
text(X(leafnodes),Y(leafnodes),ctext,'HitTest','off','Parent',ax,...
         'VerticalAlignment','top','HorizontalAlignment','center',...
          'Tag','leaflabel','Clipping','on','Visible',vis,'Interpreter','none');

lefttext = cell(length(branchnodes),1);
righttext = cell(length(branchnodes),1);
for j=1:length(branchnodes)
   k = branchnodes(j);
   cut = cutoff{k};
   if splitvar(k)>0
      varname = names{splitvar(k)};
      lefttext{j} = sprintf('%s < %g   ',varname,cut);
      righttext{j} = sprintf('  %s >= %g',varname,cut);
   else
      varname = names{-splitvar(k)};
      cats = cut{1};
      if length(cats)==1
         lefttext{j} = sprintf('%s = %s   ',varname,num2str(cats,'%g '));
      else
         lefttext{j} = sprintf('%s in (%s)   ',varname,num2str(cats,'%g '));
      end
      cats = cut{2};
      if length(cats)==1
         righttext{j} = sprintf('   %s = %s',varname,num2str(cats,'%g '));
      else
         righttext{j} = sprintf('   %s in (%s)',varname,num2str(cats,'%g '));
      end
   end
end

h = findobj(fig,'Tag','menubr');
vis = get(h,'Checked');
text(X(branchnodes),Y(branchnodes),lefttext,'HitTest','off','Parent',ax,...
     'Tag','branchlabel','Clipping','on','Visible',vis,'Interpreter','none',...
     'HorizontalAlignment','right');
text(X(branchnodes),Y(branchnodes),righttext,'HitTest','off','Parent',ax,...
     'Tag','branchlabel','Clipping','on','Visible',vis,'Interpreter','none',...
     'HorizontalAlignment','left');

% Show pruned nodes or not as desired
doprunegraph(fig);

% Adjust axes contents
dozoom(fig);

% Adjust layout of controls to fit figure
layoutfig(fig);

% ---------------------------------------------
function growprune(varargin)
%GROWPRUNE Expand or contract tree using optimal pruning sequence.

% Fetch stored information
h = gcbo;
fig = gcbf;
ud = get(fig,'UserData');
names    = ud{4};
doclass  = ud{5};
curlevel = ud{7};
fulltree = ud{6};

% Get optimal pruning sequence and current pruning level
prunelist = fulltree.prunelist;
if isequal(get(h,'Tag'),'prune')
   curlevel = min(max(prunelist),curlevel+1);
else
   curlevel = max(0,curlevel-1);
end

% Clear axes, then draw at new pruning level
ax = get(fig,'CurrentAxes');
delete(get(ax,'Children'));
[X,Y] = drawtree(fulltree,fig,names,curlevel);

% Remember everything
set(fig,'ButtonDownFcn',@removelabels, ...
        'UserData',{X Y 0 names doclass fulltree curlevel});

updateenable(fig);
updatelevel(fig,curlevel,fulltree);

% ----------------------------------------------
function updatelevel(fig,curlevel,Tree)
%UPDATELEVEL Update text display of current pruning level.

if ~isempty(Tree.prunelist)
   maxlevel = max(Tree.prunelist);
else
   maxlevel = 0;
end
h = findobj(fig,'Tag','prunelev');
set(h,'String',sprintf('%d of %d',curlevel,maxlevel));
e = get(h,'Extent');
p = get(h,'Position');
p(3) = e(3);
set(h,'Position',p);

% ----------------------------------------------
function updateenable(fig)
%UPDATEENABLE Update enabled/disabled status of buttons in figure

ud = get(fig,'UserData');
curlevel = ud{7};
fulltree = ud{6};
enableg = 'on';     % enable grow button?
enablep = 'on';     % enable prune button?
if isempty(fulltree.prunelist)
   enableg = 'off';
   enablep = 'off';
else
   maxlevel = max(fulltree.prunelist);
   if curlevel >= maxlevel
      enablep = 'off';
   end
   if curlevel <= 0;
      enableg = 'off';
   end
end
set(findobj(fig,'tag','prune'),'Enable',enablep);
set(findobj(fig,'tag','grow'),'Enable',enableg);

% ----------------------------------------------
function labelpoint(varargin)
%LABELPOINT Label point on tree in response to mouse click.

h = gcbo;
f = gcbf;
stype = get(f,'SelectionType');
if ~isequal(stype,'alt') && ~isequal(stype,'extend')
   removelabels;
end
t = get(h,'Tag');
if isequal(t,'branch') || isequal(t,'leaf')
   ud = get(f,'UserData');
   X = ud{1};       % x coordinates
   Y = ud{2};       % y coordinates
   names = ud{4};   % variable names
   doclass = ud{5}; % classification, not regression
   Tree = ud{6};    % complete tree

   if doclass
      cnames = Tree.classname;
   end
   splitvar = Tree.var;
   cutoff = Tree.cut;

   % Find closest point
   ax = get(f,'CurrentAxes');
   cp = get(ax,'CurrentPoint');
   [ignore,node] = min(abs(X-cp(1,1)) + abs(Y-cp(1,2)));

   uih = findobj(f,'Tag','clicklist');
   labeltype = get(uih,'Value');

   if isequal(labeltype,4) && doclass
      % Show fitted class probabilities
      P = Tree.classprob;
      txt = 'Class probabilities:';
      for j=1:size(P,2);
         txt = sprintf('%s\n%s = %.3g',txt,cnames{j},P(node,j));
      end

   elseif isequal(labeltype,3) && ~doclass
      % Show node statistics
      xbar = Tree.class(node);
      Nk = Tree.nodesize(node);
      if Nk > 1
         s = sqrt((Tree.nodeerr(node) * Nk) / (Nk - 1));
         txt = sprintf('N = %d\nMean = %g\nStd. dev. = %g',Nk,xbar,s);
      else
         txt = sprintf('N = %d\nMean = %g',Nk,xbar);
      end

   elseif isequal(labeltype,3) && doclass
      % Show class membership in data
      C = Tree.classcount;
      N = Tree.nodesize(node);
      txt = sprintf('Total data points = %d',N);
      for j=1:size(C,2);
         txt = sprintf('%s\n%d %s',txt,C(node,j),cnames{j});
      end

   elseif isequal(labeltype,1)
      % Get a display of the split rule at branch nodes,
      % or the majority class at leaf nodes
      if ~isequal(t,'branch')
         if doclass
            txt = sprintf('Node %d (leaf)\nClass: %s',node,cnames{Tree.class(node)});
         else
            txt = sprintf('Node %d (leaf)\nPrediction: %g',node,Tree.class(node));
         end
      elseif splitvar(node)>0
         txt = sprintf('Node %d (branch)\nRule:  %s < %g',...
                       node,names{splitvar(node)},cutoff{node});
      else
         cut = cutoff{node};
         cats = cut{1};
         if length(cats)==1
            txt = sprintf('Node %d (branch)\nRule:  %s = %s',node,...
                          names{-splitvar(node)},num2str(cats,'%g '));
         else
            txt = sprintf('Node %d (branch)\nRule:  %s in (%s)',node,...
                          names{-splitvar(node)},num2str(cats,'%g '));
         end
      end
   elseif isequal(labeltype,2)
      % Get the conditions satisfied by points in this node
      if node==1
         txt = 'Root of tree';
      else
         % Find limits for each variable along the path to this node
         nvars = max(abs(splitvar(:)));
         lims = cell(nvars,3);
         lims(:,1) = {-Inf};               % lower limit
         lims(:,2) = {Inf};                % upper limit
         lims(:,3) = num2cell((1:nvars)'); % variable number
         p = Tree.parent(node);
         c = node;
         while(p>0)
            leftright = 1 + (Tree.children(p,2)==c);
            vnum = splitvar(p);
            vcut = cutoff{p};
            if vnum>0
               if isinf(lims{vnum,3-leftright})
                  lims{vnum,3-leftright} = vcut;
               end
            else
               if ~iscell(lims{-vnum,1})
                  lims{-vnum,1} = vcut{leftright};
               end
            end
            c = p;
            p = Tree.parent(p);
         end

         % Create label listing any variable with finite limits
         txt = 'At this node:';
         for j=1:size(lims,1);
            L1 = lims{j,1};
            L2 = lims{j,2};
            if ~iscell(L1) && isinf(L1) && isinf(L2)
               continue
            end
            vnum = lims{j,3};

            if iscell(L1)
               cats = L1{1};
               if length(cats)==1
                  txt = sprintf('%s\n%s = %s',txt,names{vnum},num2str(cats,'%g '));
               else
                  txt = sprintf('%s\n%s \\in (%s)',txt,names{vnum},num2str(cats,'%g '));
               end
            elseif isinf(L1)
               txt = sprintf('%s\n%s < %g',txt,names{vnum},L2);
            elseif isinf(L2)
               txt = sprintf('%s\n%g < %s',txt,L1,names{vnum});
            else
               txt = sprintf('%s\n%g < %s < %g',txt,L1,names{vnum},L2);
            end
         end
      end
   else
      txt = '';
   end

   % Add label
   if ~isempty(txt)
      x = X(node);
      y = Y(node);
      xlim = get(gca,'xlim');
      ylim = get(gca,'ylim');
      if x<mean(xlim)
         halign = 'left';
         dx = 0.02;
      else
         halign = 'right';
         dx = -0.02;
      end
      if y<mean(ylim)
         valign = 'bottom';
         dy = 0.02;
      else
         valign = 'top';
         dy = -0.02;
      end
      h = text(x+dx*diff(xlim),y+dy*diff(ylim),txt,'Interpreter','none');
      yellow = [1 1 .85];
      set(h,'backgroundcolor',yellow,'margin',3,'edgecolor','k',...
            'HorizontalAlignment',halign,'VerticalAlignment',valign,...
            'tag','LinescanText','ButtonDownFcn',@startmovetips);
      line(x,y,'Color',yellow,'Marker','.','MarkerSize',20,...
           'Tag','LinescanMarker');
   end
end

% ----------------------------------------------
function startmovetips(varargin)
%STARTMOVETIPS Start movement of node tips.

f = gcbf;
set(f,'WindowButtonUpFcn',@donemovetips,...
      'WindowButtonMotionFcn',@showmovetips,...
      'Interruptible','off','BusyAction','queue');

o = gcbo;
p1 = get(f,'CurrentPoint');
a = get(f,'CurrentAxes');
ud = get(a,'UserData');
ud(1:2) = {o p1};
set(a,'UserData',ud);

% ----------------------------------------------
function showmovetips(varargin)
%SHOWMOVETIPS Show current movement of node tips.
domovetips(0,varargin{:});

% ----------------------------------------------
function donemovetips(varargin)
%DONEMOVETIPS Finish movement of node tips.
domovetips(1,varargin{:});

% ----------------------------------------------
function domovetips(alldone,varargin)
%DOMOVETIPS Carry out movement of node tips.

f = gcbf;
if alldone
   set(f,'WindowButtonUpFcn','','WindowButtonMotionFcn','');
end
a = get(f,'CurrentAxes');
ud = get(a,'UserData');
o = ud{1};
p1 = ud{2};
p2 = get(f,'CurrentPoint');
p0 = get(a,'Position');
pos = get(o,'Position');
dx = (p2(1)-p1(1)) * diff(get(a,'XLim')) / p0(3);
dy = (p2(2)-p1(2)) * diff(get(a,'YLim')) / p0(4);
pos(1) = pos(1) + dx;
pos(2) = pos(2) + dy;
set(o,'Position',pos);
ud{2} = p2;
set(a,'UserData',ud);

% ----------------------------------------------
function resize(varargin)
%RESIZE Resize figure showing decision tree.
layoutfig(gcbf)

% ----------------------------------------------
function layoutfig(f)
%LAYOUTFIG Layout figure contents

set(f,'Units','points');
fpos = get(f,'Position');

% Resize frame
h = findobj(f,'Tag','frame');
frpos = get(h,'Position');
frpos(2) = fpos(4) - frpos(4);
frpos(3) = fpos(3);
set(h,'Position',frpos);

% Resize controls inside frame
tags = {'clicktext'  'clicklist'  'magtext' 'maglist' ...
        'pruneframe' 'prunelabel' 'prunelev'};
mult = [1.6          1.35         1.6       1.35      ...
        1.7          1.6          1.6];
for j=1:length(tags)
   h = findobj(f,'Tag',tags{j});
   p = get(h,'Position');
   if j==1, theight = p(4); end
   p(2) = fpos(4) - mult(j)*theight;
   set(h,'Position',p);
end

h = findobj(f,'Tag','grow');
p = get(h,'Position');
p(2) = frpos(2)+2;
set(h,'Position',p);
h = findobj(f,'Tag','prune');
p(2) = p(2)+p(4);
set(h,'Position',p);

% Resize sliders
hh = findobj(f,'Tag','hslider');
hv = findobj(f,'Tag','vslider');
p1 = get(hh,'Position');
sw = p1(4);
p1(3) = frpos(3) - sw - 1;
set(hh,'Position',p1);
p2 = get(hv,'Position');
p2(1) = frpos(3) - sw - 1;
p2(4) = frpos(2) - sw - 1;
set(hv,'Position',p2);
if isequal(get(hh,'Visible'),'off')
   sw = 0;
end

% Resize axes
h = get(f,'CurrentAxes');
p = [0, sw, frpos(3)-sw, frpos(2)-sw];
set(h,'Position',p);

% ----------------------------------------------
function removelabels(varargin)
%REMOVELABELS Remove any labels remaining on the graph.

f = gcbf;
delete(findall(f,'Tag','LinescanMarker'));
delete(findall(f,'Tag','LinescanText'));

% ----------------------------------------------
function [X,Y] = layouttree(Tree,isleaf)
%LAYOUTTREE Select x,y coordinates of tree elements.
n = length(Tree.node);
X = zeros(n,1);
Y = X;
layoutstyle = 1;

% Plot top node on one level, its children at next level, etc.
for j=1:n
   p = Tree.parent(j);
   if p>0
      Y(j) = Y(p)+1;
   end
end
if layoutstyle==1
   % Layout style 1
   % Place terminal nodes, then put parents above them

   % First get preliminary placement, used to position leaf nodes
   for j=1:n
      p = Tree.parent(j);
      if p==0
         X(j) = 0.5;
      else
         dx = 2^-(Y(j)+1);
         if j==Tree.children(p,1)
            X(j) = X(p) - dx;
         else
            X(j) = X(p) + dx;
         end
      end
   end

   % Now make leaf nodes equally spaced, preserving their order
   leaves = find(isleaf);
   nleaves = length(leaves);
   [ignore,b] = sort(X(leaves));
   X(leaves(b)) = (1:nleaves) / (nleaves+1);

   % Position parent nodes above and between their children
   for j=max(Y):-1:0
      a = find(~isleaf & Y==j);
      c = Tree.children(a,:);
      X(a) = (X(c(:,1))+X(c(:,2)))/2;
   end
else
   % Layout style 2
   % Spread out the branch nodes, somewhat under their parent nodes
   X(Y==0) = 0.5;
   for j=1:max(Y)
      vis = (Y==j);                  % real nodes at this level
      invis = (Y==(j-1) & isleaf);   % invisible nodes to improve layout
      visrows = find(vis);
      nvis = sum(vis);
      nboth = nvis + sum(invis);
      x = [X(Tree.parent(vis))+1e-10*(1:nvis)'; X(invis)];
      [xx,xidx] = sort(x);
      xx(xidx) = 1:nboth;
      X(visrows) = (xx(1:nvis) / (nboth+1));
   end
end

k = max(Y);
Y = 1 - (Y+0.5)/(k+1);

% ---------------------------------------
function fig = setupfigure(doclass)
%SETUPFIGURE Set up uicontrols on decision tree figure.

fig = figure('IntegerHandle','off', 'NumberTitle','off', ...
             'Units','points','PaperPositionMode','auto',...
             'Tag','tree viewer');
ax = axes('Parent',fig,'UserData',cell(1,4),'XLim',0:1,'YLim',0:1);

% Set default print options
pt = printtemplate;
pt.PrintUI = 0;
set(fig,'PrintTemplate',pt)

if doclass
   figtitle = 'Classification tree viewer';
else
   figtitle = 'Regression tree viewer';
end

% Arrange figure contents
pos = [0 0 1 1];
set(ax,'Visible','off','XLim',0:1,'YLim',0:1,'Position',pos);
set(ax,'Units','points');
apos = get(ax,'Position');
fpos = get(fig,'Position');
hframe = uicontrol(fig,'Units','points','Style','frame',...
                   'Position',[0 0 1 1],'Tag','frame');

% Tip-related items
h=uicontrol(fig,'units','points','Tag','clicktext',...
            'String','Click to display:', 'style','text',...
            'HorizontalAlignment','left','FontWeight','bold');
extent = get(h,'Extent');
theight = extent(4);
aheight = apos(4);
tbottom = aheight - 1.5*theight;
posn = [2, tbottom, 150, theight];
set(h,'Position',posn);
textpos = posn;
e = get(h,'Extent');
if doclass
   choices = 'Identity|Variable ranges|Class membership|Estimated probabilities';
else
   choices = 'Identity|Variable ranges|Node statistics';
end
posn = [e(1)+e(3)+2, aheight-1.25*theight, 100, theight];
uicontrol(fig,'units','points','position',posn,'Tag','clicklist',...
            'String',choices, 'Style','pop','BackgroundColor',ones(1,3),...
            'Callback',@removelabels);
set(ax,'Position',[0 0 apos(3) tbottom]);
set(fig,'Toolbar','figure', 'Name',figtitle,'HandleVisibility','callback');

% Magnification items
textpos(1) = posn(1) + posn(3) + 10;
h=uicontrol(fig,'units','points','Tag','magtext','Position',textpos,...
            'String','Magnification:', 'style','text',...
            'HorizontalAlignment','left','FontWeight','bold');
e = get(h,'Extent');
textpos(3) = e(3);
set(h,'Position',textpos);
posn = [textpos(1)+textpos(3)+2, posn(2), 55, posn(4)];
h=uicontrol(fig,'units','points','position',posn,'Tag','maglist',...
            'String','x', 'Style','pop','BackgroundColor',ones(1,3),...
            'Callback',@domagnif);
adjustcustomzoom(h,false);

% Prune-related items
textpos(1) = posn(1) + posn(3) + 10;
h = uicontrol(fig,'units','points','position',textpos,'Tag','prunelabel',...
              'Style','text','HorizontalAlignment','left',...
              'FontWeight','bold','String','Pruning level:');
e = get(h,'Extent');
textpos(3) = e(3);
set(h,'Position',textpos);

posn(1) = textpos(1) + textpos(3) + 5;
posn(2) = posn(2) - 0.25*e(4);
posn(3) = 60;
posn(4) = 1.5*e(4);
textpos(1) = posn(1) + 3;
textpos(3) = posn(3) - 6;
uicontrol(fig,'Style','frame','Units','points','Position',posn,...
              'Tag','pruneframe');
uicontrol(fig,'units','points','position',textpos,'Tag','prunelev',...
              'Style','text','HorizontalAlignment','left',...
              'FontWeight','bold','String','1234 of 9999');

% Create an arrow for labeling button controls
fcolor = get(fig,'Color');
ar = ...
[1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 0 0 0 0 0 0 0 0 0 0 0 1
 1 1 0 0 0 0 0 0 0 0 0 1 1
 1 1 1 0 0 0 0 0 0 0 1 1 1
 1 1 1 1 0 0 0 0 0 1 1 1 1
 1 1 1 1 1 0 0 0 1 1 1 1 1
 1 1 1 1 1 1 0 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1
 1 1 1 1 1 1 1 1 1 1 1 1 1];
ar = repmat(ar,[1 1 3]);
ar(:,:,1) = min(ar(:,:,1),fcolor(1));
ar(:,:,2) = min(ar(:,:,2),fcolor(2));
ar(:,:,3) = min(ar(:,:,3),fcolor(3));

posn(1) = posn(1) + posn(3);
posn(4) = posn(4)/2;
posn(3) = posn(4);
uicontrol(fig,'units','points','position',posn,'Tag','prune',...
              'CData',ar(end:-1:1,:,:),...
              'Style','pushbutton','Callback',@growprune);
posn(2) = posn(2) + posn(4);
uicontrol(fig,'units','points','position',posn,'Tag','grow',...
              'CData',ar,...
              'Style','pushbutton','Callback',@growprune);

if posn(1)+posn(3) > fpos(3)
   fpos(3) = posn(1)+posn(3);
   set(fig,'Position',fpos);
   apos(3) = fpos(3);
   set(ax,'Position',apos);
end

% Adjust frame position
lowest = min(posn(2),textpos(2))-2;
frpos = apos;
frpos(4) = 1.1*(apos(4)-lowest);
frpos(2) = apos(4) - frpos(4);
set(hframe,'Position',frpos);

% Add scroll bars, start out invisible
h1 = uicontrol(fig,'Style','slider','Tag','hslider','Visible','off',...
               'Units','points','Callback',@dopan);
p1 = get(h1,'Position');
sw = p1(4);               % default slider width
p1(1:2) = 1;
p1(3) = fpos(3)-sw;
set(h1,'Position',p1);
p2 = [fpos(3)-sw, sw, sw, frpos(2)-sw];
uicontrol(fig,'Style','slider','Tag','vslider','Visible','off',...
               'Units','points','Position',p2,'Callback',@dopan);

% Add new menu before the Window menu
hw = findall(fig,'Type','uimenu', 'Label','&Window');
h0 = uimenu(fig,'Label','T&ree','Position',get(hw,'Position'));
uimenu(h0, 'Label','Show &Full Tree', 'Position',1,'Tag','menufull',...
           'Checked','on','Callback',@domenu);
uimenu(h0, 'Label','Show &Unpruned Nodes', 'Position',2,'Tag','menuunpr',...
           'Checked','off','Callback',@domenu);
uimenu(h0, 'Label','Label &Branch Nodes', 'Position',3,'Tag','menubr',...
           'Checked','on','Callback',@domenu,'Separator','on');
uimenu(h0, 'Label','Label &Leaf Nodes', 'Position',4,'Tag','menuleaf',...
           'Checked','on','Callback',@domenu);

set(fig,'ResizeFcn',@resize);

% ------------------------------------------
function domenu(varargin)
%DOMENU Carry out menu actions for tree viewer.

o = gcbo;
f = gcbf;
t = get(o,'Tag');
switch(t)
 % Change display from full tree to unpruned nodes or vice versa
 case {'menufull' 'menuunpr'}
   ischecked = isequal(get(o,'Checked'),'on');
   isfull = isequal(t,'menufull');
   if isfull
      dofull = ~ischecked;
   else
      dofull = ischecked;
   end
   mfull = findobj(f,'Type','uimenu','Tag','menufull');
   munpr = findobj(f,'Type','uimenu','Tag','menuunpr');
   if dofull
      set(mfull,'Checked','on');
      set(munpr,'Checked','off');
   else
      set(mfull,'Checked','off');
      set(munpr,'Checked','on');
   end
   doprunegraph(f,dofull);
   dozoom(f);

 % Turn on/off branch labels
 case 'menubr'
   curval = get(o,'Checked');
   if isequal(curval,'on')
      set(o,'Checked','off');
      h = findobj(f,'Type','text','Tag','branchlabel');
      set(h,'Visible','off');
   else
      set(o,'Checked','on');
      h = findobj(f,'Type','text','Tag','branchlabel');
      set(h,'Visible','on');
   end

 % Turn on/off leaf labels
 case 'menuleaf'
   curval = get(o,'Checked');
   if isequal(curval,'on')
      set(o,'Checked','off');
      h = findobj(f,'Type','text','Tag','leaflabel');
      set(h,'Visible','off');
   else
      set(o,'Checked','on');
      h = findobj(f,'Type','text','Tag','leaflabel');
      set(h,'Visible','on');
   end
end

% ------------------------------------------
function doprunegraph(f,dofull)
%DOPRUNEGRAPH Adjust graph to show full/pruned setting

a = get(f,'CurrentAxes');
h1 = findobj(a,'Type','line','Tag','prunednode');
h2 = findobj(a,'Type','line','Tag','prunedconnection');

% Figure out whether to show full tree
if nargin<2
   o = findobj(f,'Type','uimenu','Tag','menufull');
   dofull = isequal(get(o,'Checked'),'on');
end

% Adjust graph
if dofull
   set(h1,'Visible','on');
   set(h2,'Visible','on');
   xlim = get(a,'XLim');
   ylim = get(a,'YLim');
   bigxlim = 0:1;
   bigylim = 0:1;
else
   set(h1,'Visible','off');
   set(h2,'Visible','off');
   h1 = findobj(f,'Type','line','Tag','leaf');
   h2 = findobj(f,'Type','line','Tag','branch');
   x1 = get(h1,'XData');
   y1 = get(h1,'YData');
   y2 = get(h2,'YData');
   dx = 1 / (1+length(x1));
   ally = sort(unique([y1(:); y2(:)]));
   if length(ally)>1
      dy = 0.5 * (ally(2)-ally(1));
   else
      dy = 1-ally;
   end
   xlim = [min(x1)-dx, max(x1)+dx];
   ylim = [min(ally)-dy, max(ally)+dy];
   bigxlim = 0:1;
   bigylim = [ylim(1) 1];
end
axislistener(a,false);
set(a,'XLim',xlim,'YLim',ylim);
axislistener(a,true);
hh = findobj(f,'Tag','hslider');
set(hh,'UserData',bigxlim);
hv = findobj(f,'Tag','vslider');
set(hv,'UserData',bigylim);

% ------------------------------------------
function domagnif(varargin)
%DOMAGNIF React to magnification level change

f = gcbf;
o = gcbo;

% We need sliders if the magnification level is > 100%
h = [findobj(f,'Tag','hslider'), findobj(f,'Tag','vslider')];
maglevel = get(o,'Value');
if maglevel==1
   set(h,'Visible','off');
else
   set(h,'Visible','on');
end

% Adjust layout if necessary
resize;

% Adjust axes contents
dozoom(f);

% Remove custom zoom amount from list if not in use
if maglevel<=4
   adjustcustomzoom(o,false);
end

% Turn off any manual zooming
zoom(f,'off');

% --------------------------------------------------
function adjustcustomzoom(o,add)
%ADJUSTCUSTOMZOOM Add or remove special custom magnification level
nchoices = size(get(o,'String'),1);
choices = '100%|200%|400%|800%';
if ~add && nchoices~=4
   set(o,'String',choices);
elseif add && nchoices~=5
   choices = [choices '|Custom'];
   set(o,'String',choices);
end

% ------------------------------------------
function dozoom(f)
%DOZOOM Adjust axes contents to match magnification settings

a = get(f,'CurrentAxes');
hh = findobj(f,'Tag','hslider');
hv = findobj(f,'Tag','vslider');
hm = findobj(f,'Tag','maglist');

% Get information about x/y ranges and current midpoint
bigxlim = get(hh,'UserData');
bigylim = get(hv,'UserData');
xlim = get(a,'XLim');
ylim = get(a,'YLim');
currx = (xlim(1)+xlim(2))/2;
curry = (ylim(1)+ylim(2))/2;

% How much are we magnifying each axis?
magfact = [1 2 4 8];
mag = get(hm,'Value');
if mag<=4
   magfact = magfact(mag)*ones(1,2);
else
   magfact = [diff(bigxlim)/diff(xlim), diff(bigylim)/diff(ylim)];
end
magfact = max(magfact,1);

if all(magfact==1)                 % no magnification
   xlim = bigxlim;
   ylim = bigylim;
else                               % magnify by showing a subset of range
   magfact = max(magfact,1.01);
   dx = diff(bigxlim)/magfact(1);
   dy = diff(bigylim)/magfact(2);
   xval = max(bigxlim(1), min(bigxlim(2)-dx, currx-dx/2));
   xlim = xval + [0,dx];
   yval = max(bigylim(1), min(bigylim(2)-dy, curry-dy/2));
   ylim = yval + [0,dy];
   set(hh,'Min',bigxlim(1),'Max',bigxlim(2)-dx,'Value',xval);
   set(hv,'Min',bigylim(1),'Max',bigylim(2)-dy,'Value',yval);
end
axislistener(a,false);
set(a,'XLim',xlim,'YLim',ylim);
axislistener(a,true);

% ------------------------------------------
function dopan(varargin)
%DOPAN Pan around magnified tree display

f = gcbf;
a = get(f,'CurrentAxes');
o = gcbo;
val = get(o,'Value');

axislistener(a,false);
if isequal(get(o,'Tag'),'hslider')
   xlim = get(a,'XLim');
   xlim = xlim + (val-xlim(1));
   set(a,'XLim',xlim);
else
   ylim = get(a,'YLim');
   ylim = ylim + (val-ylim(1));
   set(a,'YLim',ylim);
end
axislistener(a,true);

% -----------------------------------------
function axislistener(a,enable)
%AXISLISTENER Enable or disable listening to axis limit changes

f = get(a,'Parent');
ud = get(a,'UserData');
if enable
   % Start listening to axis limit changes
   list1 = addlistener(a, 'XLim', 'PostSet', @(src,evt) customzoom(f));
   list2 = addlistener(a, 'YLim', 'PostSet', @(src,evt) customzoom(f));
   ud(3:4) = {list1 list2};
else
   % Delete listeners
   for j=3:4
      lstnr = ud{j};
      if ~isempty(lstnr) && ishandle(lstnr), delete(lstnr); end
   end
   ud(3:4) = {[]};
end
set(a,'UserData',ud);

% -----------------------------------------
function customzoom(f)
%CUSTOMPAN Deal with panning of a zoomed tree view

a = get(f,'CurrentAxes');
xlim = get(a,'XLim');
ylim = get(a,'YLim');

hh = findobj(f,'Tag','hslider');
hv = findobj(f,'Tag','vslider');
hm = findobj(f,'Tag','maglist');

bigxlim = get(hh,'UserData');
bigylim = get(hv,'UserData');
magfact = [1 2 4 8];

% Figure out if we have a standard zoom amount (100%, 200%, etc.) or
% a custom zoom amount
xratio = diff(bigxlim) / diff(xlim);
yratio = diff(bigylim) / diff(ylim);
standard = abs(xratio-yratio)<=0.02 & abs(xratio-round(xratio))<=0.02 ...
                                    & abs(yratio-round(yratio))<=0.02;
if standard
   xratio = round(xratio);
   standard = ismember(xratio,magfact);
end

% Update the magnification setting
if standard
   set(hm,'Value',find(magfact==xratio));
   adjustcustomzoom(hm,false);
   if xratio==1
      h = [findobj(f,'Tag','hslider'), findobj(f,'Tag','vslider')];
      set(h,'Visible','off');
   end
else
   adjustcustomzoom(hm,true);
   set(hm,'Value',5);
   h = [findobj(f,'Tag','hslider'), findobj(f,'Tag','vslider')];
   set(h,'Visible','on');
end

dozoom(f);

% ---------------------- helper to fix menu contents
function adjustmenu(fig)
%ADJUSTMENU Adjust contents of curve fit plot menus and toolbar

% Remove some menus entirely]
badmenus = {'&Edit' '&View' '&Insert'};
h = findall(fig, 'Type','uimenu', 'Parent',fig);
for j=1:length(badmenus)
   h0 = findall(h,'flat', 'Label',badmenus{j});
   if (~isempty(h0))
      j = find(h==h0);
      delete(h0);
      h(j) = [];
   end
end

% Add or remove some items from other menus
% Fix FILE menu
h0 = findall(h,'flat', 'Label','&File');
h1 = findall(h0, 'Type','uimenu', 'Parent',h0);
for j=length(h1):-1:1
   mlabel = get(h1(j),'Label');
   if isempty(findstr(mlabel,'Print...'))
      delete(h1(j));
      h1(j) = [];
   end
end

% Fix TOOLS menu
h0 = findall(h,'flat', 'Label','&Tools');
h1 = findall(h0, 'Type','uimenu', 'Parent',h0);
for j=length(h1):-1:1
   mlabel = get(h1(j),'Label');
   if (isempty(findstr(mlabel,'Zoom'))&& isempty(findstr(mlabel,'Pan')))
     delete(h1(j));
     h1(j) = [];
   end
end

% Fix HELP menu
h0 = findall(h,'flat', 'Label','&Help');
h1 = findall(h0, 'Type','uimenu', 'Parent',h0);
delete(h1);
uimenu(h0, 'Label', '&Help Tree Viewer', 'Position',1,'Callback',@helpviewer);

% Remove icons that don't apply here.  Keep zoom, PAN and print only.
h0 = findall(fig,'Type','uitoolbar');
h1 = findall(h0,'Parent',h0);
for j=length(h1):-1:1
   mlabel = get(h1(j),xlate('TooltipString'));
   if isempty(findstr(mlabel,'Zoom')) && isempty(findstr(mlabel,'Print')) && isempty(findstr(mlabel,'Pan'))
      delete(h1(j));
   end
end

% ---------------------- callback to display help
function helpviewer(varargin)
mapfilename = [docroot '/toolbox/stats/stats.map'];
helpview(mapfilename, 'view_classregtree_ref');

author: Cloud-CV

repo: object-proposals


text:
[afilename, apathname] = uigetfile('*.*', 'Input for Beads');
filename=fullfile(apathname,afilename);
%filename='E:\2009-05-12_GlassesBeads\tif\tryout\163m_glass40.tif'

s=40;
%a=ZEN_LSMLoad(filename);
a=readtimeseries(filename);

a=a{1};
%%Original version 20/05/2009
% ma=max(a{1},[],3);
% ExtractResults(ma);

%%Marie, 21/05/2009
% ma=max(a{1},[],3);
[pxy,bestslice]=ExtractCoordinates(a,s);

number=size(bestslice,2);
for b=1:number
    slice=squeeze(extract(a,[s s 1],[pxy(b,:) bestslice(b)]));

    FWHM(b)=FitPSF(slice);
end
sumFWHM=0;
for b=1:number
   sumFWHM=sumFWHM+FWHM(b);
end
fprintf('The average FWHM in pxl of the selected beads is\n');
averageFWHM=sumFWHM/number
fprintf('The average FWHM in nm of the selected beads is\n');
ConvertZENpxl_nm(averageFWHM)

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function res = imTangentCrop(img, pos, boxSize)
% Crop an image around a point based on local orientation.
%
%   RES = imTangentCrop(IMG, BOXCENTER, BOXSIZE)
%   Computes and orientated crop of the input image IMG, by considering all
%   pixels within a oriented box with centered given by BOXCENTER, size
%   given by BOXSIZE, and orientation evaluated from local gradient of the
%   image at the point POS.
%
%
%   Example
%
%   See also
%     imCropOrientedBox, imCropBox, imTangentCrop3d
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
% Created: 2022-06-01,    using Matlab 9.9.0.1570001 (R2020b) Update 4
% Copyright 2022 INRAE.

%% Compute transform

% evaluate gradient
grad = imLocalGradient(img, pos, 3);

% convert gradient to rotation angle
angle = atan2(grad(2), grad(1)) - pi/2;

% create the transform matrix that maps from box coords to global coords
transfo = createTranslation(pos) * createRotation(angle);

%% Sample points within box

% generate point coords along each box axis:
% * number of values equals to round(boxSize)
% * use single pixel spacing
radius = round(boxSize) / 2;
lx = -radius(1)+0.5:radius(1)-0.5;
ly = -radius(2)+0.5:radius(2)-0.5;

% map into global coordinate space
[x, y] = meshgrid(lx, ly);
[x, y] = transformPoint(x, y, transfo);

% evaluate within image
res = imEvaluate(img, x, y);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function GdDTPA
% GdDTPA reaction advection steady state
%    using MATLAB fzero for solution of implicit equation
%
%   $Ekkehard Holzbecher  $Date: 2006/05/03 $
%--------------------------------------------------------------------------
xdata = [0.22 0.42 0.84 1.66 3.33 5 10 15 20 25 30];
ydata = [54.7 54.3 54.2 54.0 52.8 52.4 51.0 49.9 48.7 47.6 46.2];
cin = 55;
aeta = 4304;
beta = 5937;
for i=1:size(xdata,2)
    yfunc(i) = fzero (@GD,cin,odeset,cin,beta,aeta,xdata(i));
end
plot (xdata,ydata,'o',xdata,yfunc);
xlabel ('distance [m]'); ylabel ('Gd concentration [ng/ml]');

function y = GD(c,cin,beta,aeta,x);
y = x+(c+2*beta*log(c)-beta*beta/c-cin-2*beta*log(cin)+beta*beta/cin)/aeta;


text:
function a = cv(algo,hyper)

%=====================================================================================
%  Cross validation object
%=====================================================================================
%  a=cv(algo,hyperParam) Returns a cv object on algorithm algo using with
%                        given Hyperparameters.
%
%  Possible hyperparameters (with defaults):
%   folds=5               - number of folds
%   repeats=1             - can do n*CV for reduced variance
%   balanced=0            - determines if cv shall be balanced  (same number of
%                            positives in each fold)
%   store_all=1           - determines if models trained in all folds shall be stored
%   output_train_error=0  - determines wheter to output training error on each fold
%                            (else cv error, i.e  test left out fold as default)
%   train_on_fold=0       - determines to  test on left out fold and train on the
%                            rest (set to true for the opposite).
%   store_trialbytrial=0  - whether to store the field a.trialbytrial, the rows of
%                            which are [i, foldnumber, y_i, f(x_i)] for each index
%                            i of a data point on which the algorithm was tested
%
%  Model:
%   child                 - stored in child algorithm
%
%  Methods:
%   train                 - selfexplanatory
%   test                  - selfexplanatory
%=====================================================================================
% Reference :
% Author    :
% Link      :
%=====================================================================================

  % <<---- Initialisation of Hyperparams ---->>
  a.folds=5;
  a.repeats=1;
  a.balanced=0;
  a.store_all=1;
  a.output_train_error=0;
  a.train_on_fold=0;
  a.store_trialbytrial=0;
  a.trialbytrial=[];

  % <<---- Initialisation of model ---->>
  if nargin>0
   a.child{1}=algo;
  else
   a.child{1}=svm;
  end

  % <<---- convert {} to group({}) ---->>
  if isa(a.child{1},'cell')
      a.child{1}=group(a.child{1});
  end;

  p=algorithm('cv');
  a= class(a,'cv',p);
  if nargin==2,
    eval_hyper;
  end;

  if a.store_all
    for i=2:a.folds*a.repeats
     a.child{i}=a.child{1};
    end
  end


text:
function mColormap = gui_Colormap_Rastafari(nSize)

% If size is not specified, set it to 256
if nargin < 1
  nSize = 256;
end

mColormap = [];

nSteps = floor(nSize / 2);
mColormap = [mColormap; gui_Interpolate([0.024169, 0.287879, 0.012572], [0.206233, 1.000000, 0.148478], nSteps)];

mColormap = [mColormap; [0.75 0.75 0.75]];

nSteps = floor(nSize / 2);
mColormap = [mColormap; gui_Interpolate([1.000000, 0.235294, 0.235294], [0.621212, 0.017977, 0.035956], nSteps)];

author: CelsoReyes


text:
function blas1_z_test ( )

%*****************************************************************************80
%
%% BLAS1_Z_TEST tests the BLAS1_Z library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 May 2006
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../blas0' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BLAS1_Z_PRB:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the BLAS1_Z library.\n' );

  blas1_z_test01 ( );
  blas1_z_test02 ( );
  blas1_z_test03 ( );
  blas1_z_test04 ( );
  blas1_z_test05 ( );
  blas1_z_test06 ( );
  blas1_z_test07 ( );
  blas1_z_test08 ( );
  blas1_z_test09 ( );

  blas1_z_test10 ( );
  blas1_z_test11 ( );
  blas1_z_test12 ( );
  blas1_z_test13 ( );
  blas1_z_test14 ( );
  blas1_z_test15 ( );
  blas1_z_test16 ( );
  blas1_z_test17 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BLAS1_Z_PRB:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../blas0' );

  return
end


text:
function dtscale = dtscale2D;

% function dtscale = dtscale2D;
% Purpose : Compute inscribed circle diameter as characteristic
%           for grid to choose timestep

Globals2D;

% Find vertex nodes
vmask1   = find( abs(s+r+2) < NODETOL)';
vmask2   = find( abs(r-1)   < NODETOL)';
vmask3   = find( abs(s-1)   < NODETOL)';
vmask  = [vmask1;vmask2;vmask3]';
vx = x(vmask(:), :); vy = y(vmask(:), :);

% Compute semi-perimeter and area
len1 = sqrt((vx(1,:)-vx(2,:)).^2+(vy(1,:)-vy(2,:)).^2);
len2 = sqrt((vx(2,:)-vx(3,:)).^2+(vy(2,:)-vy(3,:)).^2);
len3 = sqrt((vx(3,:)-vx(1,:)).^2+(vy(3,:)-vy(1,:)).^2);
sper = (len1 + len2 + len3)/2.0;
Area = sqrt(sper.*(sper-len1).*(sper-len2).*(sper-len3));

% Compute scale using radius of inscribed circle
dtscale = Area./sper;
return;

repo: Aero-matlab


text:
function fid = openLogFile(filename,mode)

% OPENLOGFILE  Open log file.
%   FID = OPENLOGFILE(FILENAME,MODE) opens the file FILENAME in MODE mode.
%   The mode MODE is a string in {'r','read','w','write'}. FID is the file
%   identifer used in other IO functions.
%
%   See also WRITELOGHEADER, WRITELOGDATA, CLOSELOGFILE, FOPEN.

switch mode
    case {'r','read'}
        m = 'r';
    case {'w','write'}
        m = 'w';
end

fid =fopen(filename,m);

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function outargs = hlp_resolve(x,default,context)
% Look up a symbol from the current hlp_scope scope.
% Data = hlp_resolve(Symbol,Default,Context)
%
% More precisely, perform a lookup for a symbol in the current (dynamic) scope formed by (possibly
% nested) hlp_scope(s). Returns a cell array of all the values associated with the symbol. Used to
% implement dynamic scoping across functions, used by parts of the expression & online system.
%
% In:
%   Symbol : a symbol (string or function handle) to look up; if omitted, the entire symbol frame
%            will be returned as a struct
%
%   Default : the default value, if the symbol does not exist (note: you may also omit the default
%             if you prefer to get an error in case the symbol does not exist)
%
%   Context : Optionally the execution context (struct with field .stack) at this point, if known.
%             Can be obtained via dbstack or try/catch.
%
% Out:
%   Data : the data associated with the symbol, or the symbol itself, if the lookup failed
%
% Examples:
%   % resolve the value of the symbol 'test' against the current dynamic scope
%   value = hlp_resolve('test')
%
% See also:
%   hlp_scope
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-05-03

% Copyright (C) Christian Kothe, SCCN, 2010, christian@sccn.ucsd.edu
%
% This program is free software; you can redistribute it and/or modify it under the terms of the GNU
% General Public License as published by the Free Software Foundation; either version 2 of the
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without
% even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with this program; if not,
% write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
% USA

global tracking;

% get stack entries
if nargin < 3
    stack = dbstack;
    entries = {stack.name};
else
    entries = {context.stack.name};
end

% check for all entries whether they refer to a stack frame...
if hlp_matlab_version < 806
    marker = 'make_func/@(f,a,frame__';
else
    marker = '@(f,a,frame__';
end
frames = find(strncmp(entries,marker,length(marker)));

% go through all stack frames...
symbol = char(x);
for k=frames
    frameid = entries{k}(length(marker)+1:end-14);
    % and check if the symbol is present there, then look up
    if isfield(tracking.stack.frames.(frameid),symbol)
        outargs = tracking.stack.frames.(frameid).(symbol);
        return;
    end
end

try
    if isfield(tracking.stack.base,symbol)
        % symbol is fetched from the base workspace
        outargs = tracking.stack.base.(symbol);
    else
        % symbol is not in the base workspace: remains unsubstituted
        outargs = default;
    end
catch
    % base workspace doesn't exist yet: create it...
    tracking.stack.base = struct();
    outargs = default;
end

repo: REST


% ==============================================================================
% 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: shows how to visulaize 2D data with FAIR.
%
% - load data             (here 'US.jpg')
% - setup viewer          (viewImage2D)
% - view  data
%==============================================================================

clear, close all, help(mfilename);
echo on

% load data, get size and specify physical domain
dataT = double(imread('US.jpg'));
m     = size(dataT);
omega = [0,m(1),0,m(2)];

% setup image viewer
viewImage('reset','viewImage','viewImage2D','colormap','gray(256)','axis','off');

% view data
viewImage(dataT,omega,m,'title','my first image');

echo off
%==============================================================================


text:
%-----------------------------------------------------------------------
% Job saved on 18-Jul-2014 13:13:47 by cfg_util (rev $Rev$)
% spm SPM - SPM12b (5672)
% cfg_basicio BasicIO - Unknown
%-----------------------------------------------------------------------
matlabbatch{1}.spm.spatial.coreg.write.ref = '<UNDEFINED>';
matlabbatch{1}.spm.spatial.coreg.write.source = '<UNDEFINED>';
matlabbatch{1}.spm.spatial.coreg.write.roptions.interp = 7;
matlabbatch{1}.spm.spatial.coreg.write.roptions.wrap = [0 0 0];
matlabbatch{1}.spm.spatial.coreg.write.roptions.mask = 0;
matlabbatch{1}.spm.spatial.coreg.write.roptions.prefix = 'r';

author: translationalneuromodeling

repo: tapas


text:
function [w b err] = myFitPlane(pc, reg, typ, inlierThresh)
% function [w b err] = myFitPlane(pc, reg, typ, inlierThresh)
% Output: Returns w'x + b = 0;
% Input: pc is Nx3 matrix.

	if(~exist('inlierThresh','var'))
		inlierThresh = 3;
	end
	X = pc(:,:,1);
	Y = pc(:,:,2);
	Z = pc(:,:,3);

	A(:,1) = X(reg);
	A(:,2) = Y(reg);
	A(:,3) = Z(reg);
	ind = ~isnan(A(:,1));
	A = A(ind,:);

	if(size(A,1) < 10),
		w = NaN(3,1); b = NaN; err = NaN(3,1);
		return;
	end

	switch typ,
		case 'hack',
				B = A;
				b = 100*ones(size(B,1),1);
				nTmp = linsolve(B,b);
				w = nTmp./norm(nTmp);
				b = -b(1)./norm(nTmp);

		case 'disparity',
				% Following XYZ paper which wants to treat disparity as the random variable!
				B = A;
				b = 100*ones(size(B,1),1)./B(:,3);
				B = bsxfun(@times,B,1./B(:,3));
				nTmp = linsolve(B,b);
				w = nTmp./norm(nTmp);
				b = -100./norm(nTmp);

		case 'ransachack',
				B = A(randperm(size(A,1)),:);
				B = A(1:min(500,end),:);
				[wb, P, inliers] = ransacfitplane(B', inlierThresh, false);
				w = wb(1:3)./norm(wb(1:3));
				b = wb(4)./norm(wb(1:3));
	end

	%Reorienting the normals of the floor.
	sg = sign(sum(sign(A*w)));
	if(isnan(sg) || sg == 0)
		sg = 1;
	end
	w = w.*sg;
	b = b*sg;

	%Find the errors - fraction of inliers, mean distance from all points, mean distance from inliers.
	dist = A*w+b;
	inliers = find(abs(dist) < inlierThresh);
	err(1) = length(inliers)/size(A,1);
	err(2) = mean(abs(dist));
	err(3) = mean(abs(dist(inliers)));
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function [ lambda, v, it_num ] = power_method2 ( n, a, x_init, it_max, tol )

%*****************************************************************************80
%
%% POWER_METHOD2 applies the power method for possibly complex eigenvalues.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 July 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Eric VanDeVelde,
%    Concurrent Scientific Programming,
%    Springer, 1994,
%    ISBN: 0-387-94195-9,
%    LC: QA76.58.V35.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(N,N), the matrix.
%
%    Input, real X(N), the initial estimate for the eigenvector.
%
%    Input, integer IT_MAX, the maximum number of iterations to take.
%    1 <= IT_MAX.
%
%    Input, real TOL, an error tolerance.
%
%    Output, complex LAMBDA, the estimate for the eigenvalue.
%
%    Output, complex V(N), the estimate for the eigenvector.
%
%    Output, integer IT_NUM, the number of iterations taken.
%
  it_num = 0;
%
%  Compute data necessary to start the iteration.
%
  x(1:n,1) = x_init(1:n);
  pi_xx = x(1:n)' * x(1:n);
  x(1:n) = x(1:n) / pi_xx;
  y(1:n,1) = a(1:n,1:n) * x(1:n);
  pi_xy = x(1:n)' * y(1:n);
  pi_yy = y(1:n)' * y(1:n);

  for it = 1 : it_max

    if ( pi_yy - pi_xy * pi_xy < tol * tol * pi_yy )
      lambda = pi_xy;
      v(1:n) = y(1:n) / sqrt ( pi_yy );
      return
    end

    z(1:n,1) = a(1:n,1:n) * y(1:n);

    pi_xz = x(1:n)' * z(1:n);
    pi_yz = y(1:n)' * z(1:n);
    pi_zz = z(1:n)' * z(1:n);

    alpha = - ( pi_yz - pi_xy * pi_xz ) / ( pi_yy - pi_xy * pi_xy );
    beta = ( pi_xy * pi_yz - pi_yy * pi_xz ) / ( pi_yy - pi_xy * pi_xy );
    gamma = pi_zz + alpha * alpha * pi_yy + beta * beta ...
      + 2.0 * ( alpha * pi_yz + beta * pi_xz + alpha * beta * pi_xy );

    if ( gamma < tol * tol * pi_zz & alpha * alpha < 4.0 * beta )

      lambda_real = - alpha / 2.0;
      lambda_imag = sqrt ( 4.0D+00 * beta - alpha * alpha ) / 2.0;
      lambda = lambda_real + lambda_imag * sqrt ( - 1.0 );

      v(1:n) = ( lambda * y(1:n) - z(1:n) ) ...
       / sqrt ( beta * pi_yy + alpha * pi_yz + pi_zz );

      return
    end

    x(1:n) = y(1:n) / sqrt ( pi_yy );
    y(1:n) = z(1:n) / sqrt ( pi_yy );

    pi_xy = pi_yz / pi_yy;
    pi_yy = pi_zz / pi_yy;

    it_num = it;

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POWER_METHOD2 - Fatal error!\n' );
  fprintf ( 1, '  Convergence was not reached.\n' );

  error ( 'POWER_METHOD2 - Fatal error!' );
end


text:
function distVal=GaussianSimilarity(mu1,Sigma1,mu2,Sigma2)
%%GAUSSIANSIMILARITY This is the normalized similarity measure used between
%               two Gaussian distributions in the Appendix of [1]. It
%               ranges from 0 (completely different) to 1 (same
%               distribution).
%
%INPUTS: mu1 The xDimX1 mean of the first Gaussian distribution.
%     Sigma1 The xDimXxDim covariance matrix of the first Gaussian
%            distribution.
%        mu2 The xDimX1 mean of the second Gaussian distribution.
%     Sigma2 The xDimXxDim covariance matrix of the second Gaussian
%            distribution.
%
%OUTPUTS: distVal The scalar distance between the two distributions.
%
%EXAMPLE:
%We show the similairty of identical distirbutions, extremely different
%distributions, and distributions that just differ a little bit.
% mu1=[12;3];
% mu2=[2000;8000];
% Sigma1=[10, 6;
%          6, 9];
% Sigma2=[8, -3;
%        -3, 14];
% %Identical Distributions (Similarity is 1).
% GaussianSimilarity(mu1,Sigma1,mu1,Sigma1)
% %Distributions that are very different (Similarity is 0).
% GaussianSimilarity(mu1,Sigma1,mu2,Sigma2)
% %Distributions that only differ a little bit (Similarity is about
% %0.2163).
% GaussianSimilarity(mu1,Sigma1,mu1*0.5,Sigma1*0.5)
%
%REFERENCES:
%[1] F. Faubel, J. McDonough, and D. Klakow, "The split and merge unscented
%    Gaussian mixture filter," IEEE Signal Processing Letters, vol. 16, no.
%    9, pp. 786-789, Sep. 2009.
%
%January 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Equation 16.
Sigma12Inv=inv(Sigma1)+inv(Sigma2);
%Equation 15.
mu12=Sigma12Inv\(Sigma1\mu1+Sigma2\mu2);

numDim=size(mu1,1);
z=zeros(numDim,1);

p1=GaussianD.PDF(z,mu1,Sigma1/2);
p2=GaussianD.PDF(z,mu2,Sigma2/2);
p12=GaussianD.PDFI(z,mu12,Sigma12Inv);

if(p1==0||p2==0||p12==0||p1*p2==0)
    distVal=0;
else
    distVal=sqrt(p1*p2)/p12;
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 rr = PRC(data, up)
%PRC fuses RR estimates using the pole ranking criterion.

%% Find times at which all three RRs were estimated
mutual_times = intersect(intersect(data.bw.t, data.fm.t), data.am.t);
quit_log = 0;
for mod = {'bw', 'am', 'fm'}
    eval(['temp_data = data.' mod{1,1} ';']);
    [rel_times,rel_els,~] = intersect(temp_data.t,mutual_times);
    if ~sum(strcmp(fieldnames(temp_data), 'pole_mag'))
        quit_log = 1;
    else
        eval([mod{1,1} '.est = data.' mod{1,1} '.v(rel_els);']);
        eval([mod{1,1} '.mag = data.' mod{1,1} '.pole_mag(rel_els);']);
        eval([mod{1,1} '.ang = data.' mod{1,1} '.pole_angle(rel_els);']);
    end
end

rr.v = nan(length(rel_times),1);
rr.t = rel_times;

% quit if one of the modulations wasn't present
if quit_log
    return
end

%% Find PRC
for win_no = 1 : length(mutual_times)
    [rel_ests, rel_mags, rel_angs] = deal([]);
    rel_mods = cell(0);
    for mod = {'bw', 'am', 'fm'}
        eval(['rel_ests = [rel_ests, ' mod{1,1} '.est(win_no)];']);
        eval(['rel_mags = [rel_mags, ' mod{1,1} '.mag(win_no)];']);
        eval(['rel_angs = [rel_angs, ' mod{1,1} '.ang(win_no)];']);
        rel_mods{length(rel_mods)+1,1} = mod{1,1};
    end
    %% Find RR
    comps = {'bwam', 'bwfm', 'amfm'};
    [prc.v, prc.rr] = deal(nan(length(comps),1));
    for comp_no = 1:length(comps)
        comp = comps(comp_no);
        % identify corresponding els
        temp1 = strfind(rel_mods, comp{1,1}(1:2));
        temp2 = strfind(rel_mods, comp{1,1}(3:4));
        rel_els = sort([find(~cellfun(@isempty,temp1)), find(~cellfun(@isempty,temp2))]);
        % find PRC for this pair
        prc.v(comp_no) = mean(rel_mags(rel_els))/abs(diff(rel_angs(rel_els)));
        prc.rr(comp_no) = mean(rel_ests(rel_els));
    end

    % identify relevant RR est
    [~,rel_comp_no] = max(prc.v);
    rr.v(win_no) = prc.rr(rel_comp_no);
end

end

author: peterhcharlton

repo: RRest


text:
function [spectrum] = specest_nanfft(dat, time, varargin)

% SPECEST_NANFFT computes a fast Fourier transform in the presence of NaNs
% in the data
%
% Use as
%   [spectrum] = specest_nanfft(dat, ...)
% where
%   dat      = matrix of chan*sample
%   time     = vector, containing time in seconds for each sample
%   spectrum = matrix of taper*chan*foi*toi of fourier coefficients
%
% Optional arguments should be specified in key-value pairs and can include:
%   basis      = precomputes set of basis functions (sines/cosines)
%   datataype  = 0, 1, 2
%
% FIXME: FFT speed not yet optimized, e.g. MATLAB version, transpose or not, ...
% FIXME: function is recursive, should be avoided in favor of transparancy
%
% See also SPECEST_MTMFFT, SPECEST_CONVOL, SPECEST_HILBERT, SPECEST_MTMCONVOL, SPECEST_MVAR, SPECEST_WAVELET

% Copyright (C) 2008, 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$

% get the optional arguments
basis     = ft_getopt(varargin, 'basis');
datatype  = ft_getopt(varargin, 'datatype');

% datatype = 0, no missing data
% datatype = 1, the missing data is at the same location for all channels
% datatype = 2, the missing data is at different timepoints for different channels

% determine the data characteristics
[nchan, nsample] = size(dat);
fsample = 1./mean(diff(time));

if mod(nsample,2)==0
  % the number of samples is even
  k = nsample/2+1;
else
  % the number of samples is odd
  k = floor(nsample/2+1);
end

% determine the type of data and thereby the most suitable algorithm to use
nancount = sum(isnan(dat), 1);
if isempty(datatype)
  if all(nancount==0)
    % there is no missing data
    datatype = 0;
  elseif all(nancount==0 | nancount==nchan)
    % the missing data is at the same location for all channels
    datatype = 1;
  else
    % the missing data is at different timepoints for different channels
    datatype = 2;
  end
end

if datatype==0
  % no basis functions are needed, because the standard FFT routine will be used

elseif datatype~=0 && isempty(basis)
  % create a separate set of basis functions for the cosine and sine
  basis_c = zeros(k, nsample);
  basis_s = zeros(k, nsample);

  % create the time axis
  t = linspace(0, 2*pi, nsample+1);
  t = t(1:end-1);

  for w=1:k
    c = cos((w-1)*t);
    s = sin((w-1)*t);
    if w==1 || (w==(k) && mod(nsample,2)==0)
      % the normalization for the lowest (DC) and the highest frequency component is different
      s = s/(nsample);
      c = c/(nsample);
    else
      s = s/(nsample/2);
      c = c/(nsample/2);
    end
    basis_c(w,:) = c;
    basis_s(w,:) = s;
  end
  % concatenate the sine and cosine basis functions
  % leaving the first and last sine functions out, since those are all zero
  if mod(nsample,2)==0
    % the number of samples is even -> the last sine wave basis function is zero
    basis = cat(1, basis_c, basis_s(2:end-1, :));
  else
    % the number of samples is odd -> also include the last sine wave basis function
    basis = cat(1, basis_c, basis_s(2:end, :));
  end
end

switch datatype
  case 0
    % there is no missing data
    % use the standard FFT implementation
    y = fft(dat, [], 2);

  case 1
    % the missing data is at the same location for all channels
    % remove that piece from the data and from the basis functions and use linear estimation

    keep = ~isnan(dat(1,:));

    if all(~keep)
      % the data is all NaN, no reason to try to estimate the basis
      % functions
      y = nan(size(dat));

    else

      basis = basis(:,keep);
      dat = dat(:,keep);

      % do the linear estimation based on dat=y*basis
      % y = dat / basis;
      y = dat * pinv(basis);

      % disentagle the estimated components

      if mod(nsample,2)==0
        % the number of samples is even -> the last sine wave basis function is zero
        sel1 = 1;       % lowest cosine, i.e. DC
        sel2 = 2:(k-1); % all cosines in between
        sel3 = k;       % highest cosine
        sel4 = (k+1):nsample; % all sines

        est1 = y(:,sel1);
        est2 = y(:,sel2);
        est3 = y(:,sel3);
        est4 = y(:,sel4);

        % combine the various estimates into a complex representation compatible with standard FFT
        y_real = cat(2, est1,        est2, est3,          fliplr(est2));
        y_imag = cat(2, zeros(nchan,1), -est4, zeros(nchan,1),    fliplr(est4));
        y = y_real + i*y_imag;

      else
        % the number of samples is odd -> also include the last sine wave basis function
        sel1 = 1;       % lowest cosine, i.e. DC
        sel2 = 2:k;     % all other cosines
        sel3 = (k+1):nsample; % all sines

        est1 = y(:,sel1);
        est2 = y(:,sel2);
        est3 = y(:,sel3);

        % combine the various estimates into a complex representation compatible with standard FFT
        y_real = cat(2, est1,        est2, fliplr(est2));
        y_imag = cat(2, zeros(nchan,1), -est3, fliplr(est3));
        y = y_real + i*y_imag;
      end

    end % if all(~keep)

  case 2
    % the missing data is at different timepoints for different channels
    % use recursion to compute the nanfft for each channel
    y = zeros(size(dat));
    for k=1:nchan
      y(k,:) = specest_nanfft(dat(k,:), time, 'basis', basis);
    end

  otherwise
    ft_error('unsupported configuration of NaNs in the data');
end

% set output
spectrum = y;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% =========================================================================
% NCSR for image denoising, Version 1.0
% Copyright(c) 2013 Weisheng Dong, Lei Zhang, Guangming Shi, and Xin Li
% All Rights Reserved.
%
% ----------------------------------------------------------------------
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is here
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
%----------------------------------------------------------------------
%
% This is an implementation of the algorithm for image interpolation
%
% Please cite the following paper if you use this code:
%
% Weisheng Dong, Lei Zhang, Guangming Shi, and Xin Li.,"Nonlocally
% centralized sparse representation for image restoration", IEEE Trans. on
% Image Processing, vol. 22, no. 4, pp. 1620-1630, Apr. 2013.
%
%--------------------------------------------------------------------------
function [pos_arr, wei_arr]  =  Block_matching(im, par)
S         =  30;
f         =  par.win;
f2        =  f^2;
nv        =  par.nblk;
s         =  par.step;
hp        =  max(12*par.nSig, par.hp);

N         =  size(im,1)-f+1;
M         =  size(im,2)-f+1;
r         =  [1:s:N];
r         =  [r r(end)+1:N];
c         =  [1:s:M];
c         =  [c c(end)+1:M];
L         =  N*M;
X         =  zeros(f*f, L, 'single');

k    =  0;
for i  = 1:f
    for j  = 1:f
        k    =  k+1;
        blk  =  im(i:end-f+i,j:end-f+j);
        X(k,:) =  blk(:)';
    end
end

I     =   (1:L);
I     =   reshape(I, N, M);
N1    =   length(r);
M1    =   length(c);
pos_arr   =  zeros(nv, N1*M1 );
wei_arr   =  zeros(nv, N1*M1 );
X         =  X';

for  i  =  1 : N1
    for  j  =  1 : M1

        row     =   r(i);
        col     =   c(j);
        off     =  (col-1)*N + row;
        off1    =  (j-1)*N1 + i;

        rmin    =   max( row-S, 1 );
        rmax    =   min( row+S, N );
        cmin    =   max( col-S, 1 );
        cmax    =   min( col+S, M );

        idx     =   I(rmin:rmax, cmin:cmax);
        idx     =   idx(:);
        B       =   X(idx, :);
        v       =   X(off, :);

        dis     =   (B(:,1) - v(1)).^2;
        for k = 2:f2
            dis   =  dis + (B(:,k) - v(k)).^2;
        end
        dis   =  dis./f2;
        [val,ind]   =  sort(dis);
        dis(ind(1))  =  dis(ind(2));

        wei         =  exp( -dis(ind(1:nv))./hp );
        wei         =  wei./(sum(wei)+eps);
        indc        =  idx( ind(1:nv) );
        pos_arr(:,off1)  =  indc;
        wei_arr(:,off1)  =  wei;
    end
end
pos_arr  = pos_arr';
wei_arr  = wei_arr';

author: lbasek

repo: image-denoising-benchmark


text:
function [x,gpd,nu]=gdpower(N,k,c)
%GDPOWER Signal with power-law group delay.
%	[X,GPD,F]=GDPOWER(N,K,C) generates a signal with a
%	power-law group delay of the form tx(f) = T0 + C*f^(K-1).
%	The output signal is of unit energy.
%
%	N   : number of points in time		(must be even)
%	K   : degree of the power-law		(default : 0)
%	C   : rate-coefficient of the power-law group delay.
%	      C must be non-zero. 		(default : 1)
%	X   : time row vector containing the modulated signal samples
%	GPD : group delay	      (length : round(N/2)-1)
%	F   : frequency bins
%
%	Examples :
%	 sig=gdpower(64,0); tfrbert(sig,1:64,0.01,0.32,256,1);
%	 [sig,gpd,f]=gdpower(256,1/2); plot(gpd,f);
%
%	See also FMPOWER.

%	O. Lemoine - April, July 1996
%	Copyright (c) 1996 by CNRS (France).
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin < 1),
 error ( 'At least one parameter.' );
elseif nargin==1,
 k=0; c=1;
elseif nargin==2,
 c=1;
end

if (N <= 0),
 error ('The signal length N must be strictly positive' );
end

t0=0;

Lnu=round(N/2);
nu=linspace(0,0.5,Lnu+1); nu=nu(2:Lnu+1);
AM=nu.^((k-2)/6);

if c==0,
 error('C must be non-zero');
end

t=1:N;
TFx=zeros(1,N);

if (k<1 & k~=0),
 d=N^k*c; t0=N/10;
 TFx(1:Lnu)=exp(-j*2*pi*(t0*(nu)+d*(nu).^k/k)).*AM;
 x=ifft(TFx).';
elseif k==0,
 d=c; t0=N/10;
 TFx(1:Lnu)=exp(-j*2*pi*(t0*(nu)+d*log(nu))).*AM;
 x=ifft(TFx).';
elseif k==1,
 t0=N;
 x=anapulse(N,t0);
elseif (k>1),
 d=N*2^(k-1)*c;
 TFx(1:Lnu)=exp(-j*2*pi*(t0*(nu)+d*(nu).^k/k)).*AM;
 x=ifft(TFx).';
end

if k~=1,
 gpd=t0+abs(sign(c)-1)/2*(N+1)+d*nu.^(k-1);
else
 gpd=t0*ones(1,N/2);
end

x=x-mean(x);
x=x/norm(x);

author: HeLiangHIT


text:
% GAINSANDREFERENCES initializes the controller parameters: gains,
%                    regularization parameters and references.
%

%% --- Initialization ---

% References for the demo with joints movements
if Config.MOVE_JOINTS

    % Postural task gains
    KP = 20*diag(ones(1,ROBOT_DOF));
    KD = 2*sqrt(KP)*0;

    Ref.Amplitude = 7.5*ones(1,ROBOT_DOF);
    Ref.Frequency = 0.5*ones(1,ROBOT_DOF);
else
    % Postural task gains
    KP = 0*diag(ones(1,ROBOT_DOF));
    KD = 0*2*sqrt(KP);

    Ref.Amplitude = 0*ones(1,ROBOT_DOF);
    Ref.Frequency = 0*ones(1,ROBOT_DOF);
end

if size(KP,1) ~= ROBOT_DOF

    error('Dimension of KP different from ROBOT_DOF')
end

author: robotology

repo: whole-body-controllers