text:
function [ value, ifault ] = ppchi2 ( p, v, g )

%*****************************************************************************80
%
%% PPCHI2 evaluates the percentage points of the Chi-squared PDF.
%
%  Discussion
%
%    Incorporates the suggested changes in AS R85 (vol.40(1),
%    pages 233-5, 1991) which should eliminate the need for the limited
%    range for P, though these limits have not been removed
%    from the routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2008
%
%  Author:
%
%    Original FORTRAN77 version by Donald Best, DE Roberts.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Donald Best, DE Roberts,
%    Algorithm AS 91:
%    The Percentage Points of the Chi-Squared Distribution,
%    Applied Statistics,
%    Volume 24, Number 3, 1975, pages 385-390.
%
%  Parameters:
%
%    Input, real P,  value of the chi-squared cumulative
%    probability density function.
%    0.000002 <= P <= 0.999998.
%
%    Input, real V, the parameter of the chi-squared probability
%    density function.
%    0 < V.
%
%    Input, real G, the value of log ( Gamma ( V / 2 ) ).
%
%    Output, real VALUE, the value of the chi-squared random
%    deviate with the property that the probability that a chi-squared random
%    deviate with parameter V is less than or equal to PPCHI2 is P.
%
%    Output, integer IFAULT, is nonzero if an error occurred.
%    0, no error.
%    1, P is outside the legal range.
%    2, V is not positive.
%    3, an error occurred in GAMMAD.
%    4, the result is probably as accurate as the machine will allow.
%
  aa = 0.6931471806;
  c1 = 0.01;
  c2 = 0.222222;
  c3 = 0.32;
  c4 = 0.4;
  c5 = 1.24;
  c6 = 2.2;
  c7 = 4.67;
  c8 = 6.66;
  c9 = 6.73;
  c10 = 13.32;
  c11 = 60.0;
  c12 = 70.0;
  c13 = 84.0;
  c14 = 105.0;
  c15 = 120.0;
  c16 = 127.0;
  c17 = 140.0;
  c18 = 175.0;
  c19 = 210.0;
  c20 = 252.0;
  c21 = 264.0;
  c22 = 294.0;
  c23 = 346.0;
  c24 = 420.0;
  c25 = 462.0;
  c26 = 606.0;
  c27 = 672.0;
  c28 = 707.0;
  c29 = 735.0;
  c30 = 889.0;
  c31 = 932.0;
  c32 = 966.0;
  c33 = 1141.0;
  c34 = 1182.0;
  c35 = 1278.0;
  c36 = 1740.0;
  c37 = 2520.0;
  c38 = 5040.0;
  e = 0.5E-06;
  maxit = 20;
  pmax = 0.999998;
  pmin = 0.000002;
%
%  Test arguments and initialize.
%
  value = - 1.0;

  if ( p < pmin || pmax < p )
    ifault = 1;
    return
  end

  if ( v <= 0.0 )
    ifault = 2;
    return
  end

  ifault = 0;
  xx = 0.5 * v;
  c = xx - 1.0;
%
%  Starting approximation for small chi-squared
%
  if ( v < - c5 * log ( p ) )

    ch = ( p * xx * exp ( g + xx * aa ) )^( 1.0 / xx );

    if ( ch < e )
      value = ch;
      return
    end
%
%  Starting approximation for V less than or equal to 0.32
%
  elseif ( v <= c3 )

    ch = c4;
    a = log ( 1.0 - p );

    while ( 1 )

      q = ch;
      p1 = 1.0 + ch * ( c7 + ch );
      p2 = ch * ( c9 + ch * ( c8 + ch ) );

      t = - 0.5 + (c7 + 2.0 * ch ) / p1 ...
        - ( c9 + ch * ( c10 + 3.0 * ch ) ) / p2

      ch = ch - ( 1.0 - exp ( a + g + 0.5 * ch + c * aa ) * p2 / p1) / t;

      if ( abs ( q / ch - 1.0 ) <= c1 )
        break
      end

    end

  else
%
%  Call to algorithm AS 111 - note that P has been tested above.
%  AS 241 could be used as an alternative.
%
   [ x, ifault ] = ppnd ( p );
%
%  Starting approximation using Wilson and Hilferty estimate
%
    p1 = c2 / v;
    ch = v * ( x * sqrt ( p1 ) + 1.0 - p1)^3;
%
%  Starting approximation for P tending to 1.
%
    if ( c6 * v + 6.0 < ch )
       ch = - 2.0 * ( log ( 1.0 - p ) - c * log ( 0.5 * ch ) + g );
    end

  end
%
%  Call to algorithm AS 239 and calculation of seven term
%  Taylor series
%
  for i = 1 : maxit

    q = ch;
    p1 = 0.5 * ch;
    [ temp, if1 ] = gammad ( p1, xx );
    p2 = p - temp;

    if ( if1 ~= 0 )
      ifault = 3;
      return
    end

    t = p2 * exp ( xx * aa + g + p1 - c * log ( ch ) );
    b = t / ch;
    a = 0.5 * t - b * c;
    s1 = ( c19 + a * ( c17 + a * ( c14 + a * ( c13 + a * ( c12 + ...
      c11 * a ))))) / c24;
    s2 = ( c24 + a * ( c29 + a * ( c32 + a * ( c33 + c35 * a )))) / c37;
    s3 = ( c19 + a * ( c25 + a * ( c28 + c31 * a ))) / c37;
    s4 = ( c20 + a * ( c27 + c34 * a) + c * ( c22 + a * ( c30 + c36 * a ))) ...
      / c38;
    s5 = ( c13 + c21 * a + c * ( c18 + c26 * a )) / c37;
    s6 = ( c15 + c * ( c23 + c16 * c )) / c38;
    ch = ch + t * ( 1.0 + 0.5 * t * s1 - b * c * ( s1 - b * ...
      ( s2 - b * ( s3 - b * ( s4 - b * ( s5 - b * s6 ))))));

    if ( e < abs ( q / ch - 1.0 ) )
       value = ch;
       return
    end

  end

  ifault = 4;
  value = ch;

  return
end


text:
function OUT = TabPrint(DATA,hlabel,vlabel,approx)
% =======================================================================
% Prints a numerical table with labels, with specified numbers of decimal
% digits
% =======================================================================
% OUT = TabPrint(DATA,hlabel,vlabel,approx)
% -----------------------------------------------------------------------
% INPUT
%   - DATA = a (TxN) matrix of numbers
% -----------------------------------------------------------------------
% OPTIONAL INPUT
%	- hlabel = a (Tx1) vector of horizontal labels
%	- vlabel = a (1xT) vector of vertical labels
%   - approx = number of decimal digits. Default = 2
% -----------------------------------------------------------------------
% OUTPUT
%   - OUT = a cell array with the formatted table
% -----------------------------------------------------------------------
% EXAMPLE
%   x = [1 2; 3 4; 5 6; 7 8; 9 10];
%   hlab = {'a';'b';'c';'d';'e';}
%   vlab = {'A','B'}
%   OUT = TabPrint(x,hlab,vlab)
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

[n, m] = size(DATA);

if ~exist('hlabel','var') || isempty(hlabel)
    hlabel = cell(1,m);
end

if ~exist('vlabel','var') || isempty(vlabel)
    vlabel = cell(n,1);
end

if ~exist('approx','var')
    approx = 2;
end

if length(hlabel)~=m
    error('ERROR: horizontal label has wrong dimension')
end

if length(vlabel)~=n
    error('ERROR: vertical label has wrong dimension')
end

OUT = roundnum2cell(DATA,approx); % Trasform matrix in cell wih approx decimal digits

OUT = [hlabel ; OUT];
aux = cell(rows(hlabel),cols(vlabel));
vlabel = [aux; vlabel];
OUT = [vlabel, OUT];

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function [xformVAnatToAcpc] = dtiXformVanatCompute(dtiT1, dtiT1ToAcpcXform, vAnatomy, vAnatMm, vAnatTal, scaleFlag)
%
% [xformVAnatToAcpc] = dtiXformVanatCompute(dtiT1, dtiT1ToAcpcXform, vAnatomy, vAnatMm, [vAnatTal], [scaleFlag=false])
%
% Computes the affine transform that will convert coordinates from the DTI
% ac-pc space to the vAnatomy space. Ie.
%   xform * [vAnatCoords 1]' = dtiImgCoords
%   inv(xform) * [dtiImgCoords 1]' = vAnatCoords
%
%
% HISTORY:
% 2004.04.30 RFD (bob@white.stanford.edu) wrote it.
% 2006.08.11 RFD: cleaned up a bit and changed calling
% convention. vAnatTal is now optional. A reasonable guess is used
% if it isn't passed. One more thing- the output xform is
% completely different! It now does something much more logical- it
% converts vAnatomy coords to acpc coords.

if(~exist('vAnatTal','var')), vAnatTal = []; end
if(~exist('scaleFlag','var')||isempty(scaleFlag)), scaleFlag = false; end

disp('Coregistering (using SPM tools)...');
%figure(999); hist(t1.cannonical_img(:),255); [clip,y] = ginput(2);
dtiT1 = double(dtiT1);
clip = [0,max(dtiT1(:))];
VG.uint8 = uint8(round((dtiT1-min(clip))*(255/(max(clip)-min(clip)))));
VG.mat = dtiT1ToAcpcXform;
%[imgSlice,x,y,z] = dtiGetSlice(VG.mat, double(VG.uint8), 3, 0);
%figure; imagesc(imgSlice); axis image xy; colormap gray;
%VG.mat = diag([sqrt(sum(dtiAcpcXform(1:3,1:3).^2)) 1]);
%figure(999); hist(b0.cannonical_img(:),255); [x,y] = ginput(2);

VF.uint8 = uint8(round(vAnatomy));
if(~isempty(vAnatTal))
  % GET XFORM FROM TAL FILE
  % We want a transform that will go into ac-pc space (ie.
  % Talairach, but without any scaling). We can compute that by
  % adjusting the scale factors in the talairach transform.
  swapXY = [0 1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1];
  vAnatAcpcXform = vAnatTal.vol2Tal.transRot'*swapXY;
  [trans,rot,scale,skew] = affineDecompose(vAnatAcpcXform);
  % We also have to rescale the translations.
  scaleDiff = vAnatMm./scale;
  trans = trans.*abs(scaleDiff);
  scale = vAnatMm.*sign(scale);
  vAnatAcpcXform = affineBuild(trans,rot,scale,skew);
  VF.mat = vAnatAcpcXform;
  %VF.mat = diag([vAnatMm 1]);
else
  % Guess a reasonable vAnatomy xform
  origin = size(VF.uint8)./2;
  % Gleaned from trial and error:
  VF.mat = [fliplr(diag(vAnatMm([3 2 1]).*[1 -1 -1])), -origin([3 2 1])'.*[1 -1 -1]' ; [0 0 0 1]];
  %VF.mat = affineBuild([0 0 0],[0 0 pi],[1 1 1])*VF.mat;
end

f.cost_fun = 'nmi';
f.sep = [4 2];
f.fwhm = [7 7];
f.params = [0 0 0 0 0 0];
if(scaleFlag) f.params(7:9) = [1 1 1]; end
rotTrans = spm_coreg(VG,VF,f);
xformVAnatToAcpc = spm_matrix(rotTrans(end,:));
xformVAnatToAcpc = xformVAnatToAcpc\VF.mat;
%xformVAnatToAcpc = xform*VG.mat\VF.mat;
%xform = inv(VF.mat\spm_matrix(rotTrans(:)'));

return;

[f,p] = uigetfile('*.dat','Select vAnatomy file...');
[vAnatomy,vAnatMm] = readVolAnat(fullfile(p,f));
h = guidata(gcf);
dtiT1 = h.bg(4).img;
dtiT1ToAcpcXform = dtiGet(h,'t1toacpcxform');
xformVAnatToAcpc = dtiXformVanatCompute(dtiT1, dtiT1ToAcpcXform, vAnatomy, vAnatMm);
h.xformVAnatToAcpc = xformVAnatToAcpc;
guidata(gcf, h);

figure; imagesc(dtiGetSlice(dtiT1ToAcpcXform,dtiT1,3,0)); axis equal tight
figure; imagesc(dtiGetSlice(xformVAnatToAcpc,vAnatomy,3,0)); axis equal tight

figure; imagesc(dtiGetSlice(VF.mat,vAnatomy,1,0)); axis equal tight

[dti,x,y,z] = dtiGetSlice(diag([sqrt(sum(dtiAcpcXform(1:3,1:3).^2)) 1]), dtiT1, 3, 20, 0);
figure; imagesc(dti); axis equal tight
[vanat,x2,y2,z2] = dtiGetSlice(inv(xform), vAnatomy, 3, 20, 0);
figure; imagesc(vanat); axis equal tight

figure; imagesc(vAnatomy(:,:,60)); axis equal tight
[dti,x,y,z] = dtiGetSlice(xform, dtiT1, 3, 60, 0);
figure; imagesc(dti); axis equal tight

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function data_num = p01_data_num ( )

%*****************************************************************************80
%
%% P01_DATA_NUM returns the number of data points for problem p01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 August 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer DATA_NUM, the number of data points.
%
  data_num = 18;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function g = orderedNoiseGradientParam(noise, mu, varsigma, y)

% ORDEREDNOISEGRADIENTPARAM Gradient of ORDERED noise's parameters.
% FORMAT
% DESC computes the gradient of the log Z of the ordered categorical noise model with respect to the of functions with respect to the
% ordered categorical
% noise's parameters.
% ARG noise : the noise structure for which the gradients are being
% computed.
% ARG mu : the input means for which the gradients are being computed.
% ARG varSigma : the input variances for which the gradients are being computed.
% ARG y : the target values for the noise model.
% RETURN g : gradients of the log Z with respect to
% the noise parameters. The ordering of the vector should match
% that provided by the function noiseExtractParam.
%
%
% SEEALSO orderedNoiseParamInit, orderednoiseGradVals, noiseGradientParam
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

D = size(y, 2);
c = 1./sqrt(noise.variance + varsigma);
gnoise.bias = zeros(1, D);
gnoise.widths = zeros(noise.C-2, 1);
for j = 1:D
  % Do lowest category first
  index = find(y(:, j)==0);
  if ~isempty(index)
    mu(index, j) = mu(index, j) + noise.bias(j) ;
    mu(index, j) = mu(index, j).*c(index, j);
    gnoise.bias(j) = gnoise.bias(j) - sum(c(index, j).*gradLogCumGaussian(-mu(index, j)));
  end

  % Intermediate categories
  index = find(y(:, j)>0 & y(:, j) <noise.C-1);
  if ~isempty(index)
    for i = index'
      mu(i, j) = mu(i, j) + noise.bias(j) - sum(noise.widths(1:y(i, j)-1));
    end
    u = mu(index, j).*c(index, j);
    uprime = (mu(index, j) - noise.widths(y(index, j))).*c(index, j);
    B1 = gaussOverDiffCumGaussian(u, uprime, 1);
    B2 = gaussOverDiffCumGaussian(u, uprime, 2);
    gnoise.bias(j) = gnoise.bias(j) + sum(c(index, j).*(B1 - B2));
    for cat = 1:noise.C-2

      subIndex = find(y(index, j) == cat);
      if ~isempty(subIndex)
        addpart = sum(c(index(subIndex), j)...
                      .*B2(subIndex));
        gnoise.widths(1:cat) = gnoise.widths(1:cat) ...
            + repmat(addpart, cat, 1);
        if(cat > 1)
          addpart = sum(c(index(subIndex), j)...
                        .*B1(subIndex));
          gnoise.widths(1:cat-1) = gnoise.widths(1:cat-1) ...
              - repmat(addpart, cat-1, 1);
        end
      end
    end
  end

  % Highest category
  index = find(y(:, j) == noise.C-1);
  if ~isempty(index)
    for i = index'
      mu(i, j) = mu(i, j) + noise.bias(j) - sum(noise.widths(1:y(i, j)-1));
    end
    mu(index, j) = mu(index, j).*c(index, j);
    addpart = sum(c(index, j).*gradLogCumGaussian(mu(index, j)));
    gnoise.bias(j) = gnoise.bias(j) + addpart;
    if length(noise.widths > 0)
      gnoise.widths = gnoise.widths ...
          - repmat(addpart, noise.C-2, 1);
    end
  end
end
if length(noise.widths>0)
  g = [gnoise.bias gnoise.widths(:)'];
else
  g = gnoise.bias;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function  [Z] =  MCWNNM_ADMM2( Y, NSig, Par )
% This routine solves the following weighted nuclear norm optimization problem with column weights,
%
% min_{X, Z} ||W(Y-X)||_F^2 + ||Z||_w,*  s.t.  X = Z
%
% Inputs:
%        Y      -- 3p^2 x M dimensional noisy matrix, D is the data dimension, and N is the number of image patches.
%        NSig -- 3p^2 x 1 dimensional vector of weights
%        Par   -- structure of parameters
% Output:
%        Z      -- 3p^2 x M dimensional denoised matrix

% tol = 1e-8;
if ~isfield(Par, 'maxIter')
    Par.maxIter = 10;
end
if ~isfield(Par, 'rho')
    Par.rho = 1;
end
if ~isfield(Par, 'mu')
    Par.mu = 1;
end
if ~isfield(Par, 'display')
    Par.display = true;
end
% Initializing optimization variables
% Intialize the weight matrix W
W = 1 ./ (NSig+eps);
% Initializing optimization variables
X = zeros(size(Y));
Z = zeros(size(Y));
A = zeros(size(Y));
%% Start main loop
iter = 0;
PatNum       = size(Y,2);
TempC  = Par.Constant * sqrt(PatNum);
while iter < Par.maxIter
    iter = iter + 1;

    % update X, fix Z and A
    % min_{X} ||W * Y - W * X||_F^2 + 0.5 * rho * ||X - Z + 1/rho * A||_F^2
    X = diag(1 ./ (W.^2 + 0.5 * Par.rho)) * (diag(W.^2) * Y + 0.5 * Par.rho * Z - 0.5 * A);

    % update Z, fix X and A
    % min_{Z} ||Z||_*,w + 0.5 * rho * ||Z - (X + 1/rho * A)||_F^2
    Temp = X + A/Par.rho;
    [U, SigmaTemp, V] =   svd(full(Temp), 'econ');
    [SigmaZ, svp] = ClosedWNNM(diag(SigmaTemp), 2/Par.rho*TempC, eps);
    Z =  U(:, 1:svp) * diag(SigmaZ) * V(:, 1:svp)';
    %     % check the convergence conditions
    %     stopC = max(max(abs(X - Z)));
    %     if Par.display && (iter==1 || mod(iter,10)==0 || stopC<tol)
    %         disp(['iter ' num2str(iter) ',mu=' num2str(Par.mu,'%2.1e') ...
    %             ',rank=' num2str(rank(Z,1e-4*norm(Z,2))) ',stopALM=' num2str(stopC,'%2.3e')]);
    %     end
    %     if stopC < tol
    %         break;
    %     else
    % update the multiplier A, fix Z and X
    A = A + Par.rho * (X - Z);
    Par.rho = min(1e4, Par.mu * Par.rho);
    %     end
end
return;

author: csjunxu


text:
% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
xg = [5 5]; % goal position
clf
axis([0 10 0 10]);
hold on
xyzlabel
grid on
xc = 5; yc = 5;
N = 8;
radius = 3;

for i=1:N
    th = (i-1)*2*pi/N;
    x0 = [xc+radius*cos(th) yc+radius*sin(th) th+pi/2];

    plot_vehicle(x0, 'r');
    r = sim('sl_drivepoint');
    y = r.find('yout');
    plot(y(:,1), y(:,2));
end
plot(xg(1), xg(2), '*')

author: Allopart

repo: rbpf-gmapping


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function [S, u, d, df, p] = createbintree(S0, T, n, r, sigma)
% creates a simple binomial tree using repmat
    dt = T / n;                  % length of one period
    u = exp(sigma * sqrt(dt));   % up move
    d = 1 / u;                   % down move
    df = exp(-r * dt);           % discount
    p = (1/df - d) / (u - d);    % probability up

    S = zeros(2^n, n+1);
    S(:,1) = S0*ones(2^n,1);

    for i=1:1:n
        a = [u * ones(2^(n-i),1); d * ones(2^(n-i),1)];
        S(:,i+1) = S(:,i).*repmat(a, 2^(i-1),1);
    end

end

author: Sable

repo: mcbench-benchmarks


% gps.m
%   Compute the output of gps sensor
%
%  Revised:
%   3/5/2010 - RB
%   5/14/2010 - RB

function y = gps(uu, P)

    % relabel the inputs
    Va      = uu(1);
%    alpha   = uu(2);
%    beta    = uu(3);
    wn      = uu(4);
    we      = uu(5);
%    wd      = uu(6);
    pn      = uu(7);
    pe      = uu(8);
    pd      = uu(9);
%    u       = uu(10);
%    v       = uu(11);
%    w       = uu(12);
%    phi     = uu(13);
%    theta   = uu(14);
    psi     = uu(15);
%    p       = uu(16);
%    q       = uu(17);
%    r       = uu(18);
    t       = uu(19);

    persistent v_n;
    persistent v_e;
    persistent v_h;

    if t==0
        v_n = 0;
        v_e = 0;
        v_h = 0;
    else
        v_n = exp(-P.k_gps * P.Ts_gps) * v_n + P.sigma_gps_n * randn;
        v_e = exp(-P.k_gps * P.Ts_gps) * v_e + P.sigma_gps_e * randn;
        v_h = exp(-P.k_gps * P.Ts_gps) * v_h + P.sigma_gps_altitude * randn;
    end

    % construct North, East, and altitude GPS measurements
    y_gps_n = pn + v_n;
    y_gps_e = pe + v_e;
    y_gps_h = -pd + v_h;

    % construct groundspeed and course measurements
    V_n = Va * cos(psi) + wn;
    V_e = Va * sin(psi) + we;
    V_g = sqrt(V_n^2 + V_e^2);

    sigma_course = P.sigma_gps_V_g / V_g;

    y_gps_Vg     = sqrt(V_n^2 + V_e^2) + P.sigma_gps_V_g * randn;
    y_gps_course = atan2(V_e, V_n) + sigma_course * randn;

    % construct total output
    y = [...
        y_gps_n;...
        y_gps_e;...
        y_gps_h;...
        y_gps_Vg;...
        y_gps_course;...
        ];

end


text:
% [K, dK_logtheta, dK_x2] = gpcov(x1, x2, logtheta)
%
function [K, dK_logtheta, dK_x2] = gpcovDot(logtheta, x1, x2)

warning('Under development, not ready')

% Covariance matrix
K = x1' * x2;

% Gradient for hyperparameters
if nargout >= 2
  dK_logtheta = nan;
end

% Gradients for inputs x2
if nargout >= 3
  if isempty(x2)
    error('Can''t calculate gradient: x2 not given');
  end
  d = rows(x2); % dimensionality of inputs
  m = cols(x1); % number of other inputs
  n = cols(x2); % number of inputs
  dK_x2 = zeros([d,m,n]);
  for j=1:n
    dK_x2(:,:,j) = x1;
  end
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function s = gen(mode, len, frequency, samplerate, amplitude, startphase, yunit)

% signal/gen
%
% generate various kind of signals with time axis
%
% si = gen(mode, len, frequency, samplerate, amplitude, startphase, yunit)
%
% len is measured in seconds
% frequency, samplerate in Hertz
% amplitude is measured in the current yunit (default is dimensionless)
% startphase in radians
% yunit is an optional argument
%
% mode may be one of the following :
% sine              - sinusodial signal between [-amplitude +amplitude]
% unoise            - uniform distrubuted random numbers in interval [-amplitude +amplitude]
% gnoise            - gaussian random numbers with zero mean
% fnoise            - spectral flat noise with zero mean (DC value = 0)
% rect				- rectangular signal with values discret values 0 and amplitude
% symmrect          - rectangular signal with values discret values -amplitude and amplitude
% comb              - train of (one sample wide) spikes of height amplitude, otherwise signal is zero
%
%
% cmerk Feb. 1998

if nargin < 1
	mode = 'sine';		% sine, uniform noise, gaussian noise, rect, puls ...
end
if nargin < 2
	len = 1;		% one second
end
if nargin < 3
	frequency = 1000;	% 1000 Hz
end
if nargin < 4
	samplerate = 8000;	% Hertz
end
if nargin < 5
	amplitude = 1.0;
end
if nargin < 6
	startphase = 0; 	% radian, 0..2pi
end
if nargin < 7
	yunit = unit;		% no dimension
end

len = ceil(len*samplerate);
xfirst = 0;
optinaltext = '';
switch mode
case	'sine'
	t = (2*pi*frequency/samplerate)*(0:len-1)' + startphase;
	tmp = amplitude * sin(t);
	optinaltext = [' with frequency : ' num2str(frequency) 'Hz'];
case	'unoise'
	tmp = (2* amplitude) * (rand(len,1)-0.5);
case	'gnoise'
	tmp = amplitude * gauss(len);
case	'fnoise'
	if mod(len, 2)
		N = (len-1)/2;
		r = exp(2*pi*i*rand(N,1));
		tmp = [0; r; conj(r(end:-1:1))];
	else							% even length
		N = (len/2)-1;
		r = exp(2*pi*i*rand(N,1));
		tmp = [0; r; 1; conj(r(end:-1:1))];
	end
	tmp = (amplitude/2) * real(ifft(tmp));
case	'rect'
	tolerance = 1e-8;
	t = (2*pi*frequency/samplerate)*(0:len-1)' + startphase;
	tmp = sin(t);
	ind = find(tmp>tolerance);
	tmp(ind) = amplitude;
	ind = find(abs(tmp)<=tolerance);
	tmp(ind(1:2:end)) = amplitude;	% zeros values are set to amplitude zero amplitude zero ...
	tmp(ind(2:2:end)) = 0;
	ind = find(tmp<tolerance);
	tmp(ind) = 0;
	optinaltext = [' with frequency : ' num2str(frequency) 'Hz'];
case	'symmrect'
	tolerance = 1e-8;
	t = (2*pi*frequency/samplerate)*(0:len-1)' + startphase;
	tmp = sin(t);
	ind = find(tmp>tolerance);
	tmp(ind) = amplitude;
	ind = find(abs(tmp)<=tolerance);
	tmp(ind(1:2:end)) = amplitude;	% zeros values are set to amplitude zero amplitude zero ...
	tmp(ind(2:2:end)) = 0;
	ind = find(tmp<tolerance);
	tmp(ind) = -amplitude;
	optinaltext = [' with frequency : ' num2str(frequency) 'Hz'];
case	'comb'
	t = (pi*frequency/samplerate)*(0:len-1)' + startphase;
	tmp = sin(t);
	zerocrossing = tmp(1:end-1) .* tmp(2:end);
	ind = find(zerocrossing<0);
	tmp = zeros(len,1);
	tmp(ind) = amplitude;
	if mod(startphase,2*pi) == 0
		tmp(1) = amplitude;
	end
	optinaltext = [' with frequency : ' num2str(frequency) 'Hz'];
case	'sinc'
	tolerance = 1e-10;
	t = (2* pi*frequency/samplerate)*(0:len-1)' + startphase;
	ind = find(abs(t)<tolerance);
	t(ind) = 1;
	tmp = amplitude * (sin(t)./t);
	tmp(ind) = amplitude;
	optinaltext = [' with frequency : ' num2str(frequency) 'Hz'];
otherwise
	error(['mode ' mode ' not supported']);
end

s = signal(tmp, ['Generated ' mode ' signal' optinaltext]);
s = setaxis(s, 1,achse(unit('s'), xfirst, 1/samplerate));
s = setyunit(s, yunit);
s = settype(s, 'scalar timeseries');

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function EIN = PQIntNoise (f)
% Generate the internal noise energy vector

% P. Kabal $Revision: 1.1 $  $Date: 2003/12/07 13:34:10 $

N = length (f);
for (m = 0:N-1)
    INdB = 1.456 * (f(m+1) / 1000)^(-0.8);
    EIN(m+1) = 10^(INdB / 10);
end

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
function [resample] = resampledesign(cfg, design)

% RESAMPLEDESIGN returns a resampling matrix, in which each row can be
% used to resample either the original design matrix or the original data.
% The random resampling is done given user-specified constraints on the
% experimental design, e.g. to swap within paired observations but not
% between pairs.
%
% Use as
%   [resample] = randomizedesign(cfg, design)
% where the configuration can contain
%   cfg.resampling       = 'permutation' or 'bootstrap'
%   cfg.numrandomization = number (e.g. 300), can be 'all' in case of two conditions
%   cfg.ivar             = number or list with indices, independent variable(s)
%   cfg.uvar             = number or list with indices, unit variable(s)
%   cfg.wvar             = number or list with indices, within-cell variable(s)
%   cfg.cvar             = number or list with indices, control variable(s)
%
% The "Independent variable" codes the condition number. Since the data is
% assumed to be independent from the condition number any reshuffeling of
% the condition number is allowed and ivar does NOT affect the resampling
% outcome.
%
% The "Unit of observation variable" corresponds to the subject number (in a
% within-subject manipulation) or the trial number (in a within-trial
% manipulation). It is best understood by considering that it corresponds
% to the "pairing" of the data in a paired T-test or repeared measures
% ANOVA. The uvar affects the resampling outcome in the way that only
% resamplings within one unit of observation are returned (e.g. swap
% conditions within a subject, not over subjects).
%
% The "Within-cell variable" corresponds to the grouping of the data in
% cells, where the multiple observations in a groups should not be broken
% apart. This for example applies to multiple tapers in a spectral estimate
% of a single trial of data (the "rpttap" dimension), where different
% tapers should not be shuffled separately. Another example is a blocked
% fMRI design, with a different condition in each block and multiple
% repetitions of the same condition within a block. Assuming that there is
% a slow HRF that convolutes the trials within a block, you can shuffle the
% blocks but not the individual trials in a block.
%
% The "Control variable" can be seen as the opposite from the within-cell
% variable: it allows you to specify blocks in which the resampling should
% be done, at the same time controlling that repetitions are not shuffled
% between different control blocks.
%
% See also FT_STATISTICS_MONTECARLO

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

ft_checkopt(cfg, 'ivar', {'numericscalar', 'numericvector', 'empty'});
ft_checkopt(cfg, 'uvar', {'numericscalar', 'numericvector', 'empty'});
ft_checkopt(cfg, 'wvar', {'numericscalar', 'numericvector', 'empty'});
ft_checkopt(cfg, 'cvar', {'numericscalar', 'numericvector', 'empty'});
ft_checkopt(cfg, 'numrandomization', {'numericscalar', 'char'});
ft_checkopt(cfg, 'resampling', {'char'});

cfg.ivar = ft_getopt(cfg, 'ivar'); % the default is empty
cfg.uvar = ft_getopt(cfg, 'uvar'); % the default is empty
cfg.wvar = ft_getopt(cfg, 'wvar'); % the default is empty
cfg.cvar = ft_getopt(cfg, 'cvar'); % the default is empty

% the default is 'no', consistent with the situation prior to 14 Sept 2011
% FIXME this needs finalization
efficient = ft_getopt(cfg, 'efficient', 'no');

% if size(design,1)>size(design,2)
%   % this function wants the replications in the column direction
%   % the matrix seems to be transposed
%   design = transpose(design);
% end

Nvar  = size(design,1);   % number of factors or regressors
Nrepl = size(design,2);   % number of replications

if ~isempty(intersect(cfg.ivar, cfg.uvar)), ft_warning('there is an intersection between cfg.ivar and cfg.uvar'); end
if ~isempty(intersect(cfg.ivar, cfg.wvar)), ft_warning('there is an intersection between cfg.ivar and cfg.wvar'); end
if ~isempty(intersect(cfg.ivar, cfg.cvar)), ft_warning('there is an intersection between cfg.ivar and cfg.cvar'); end
if ~isempty(intersect(cfg.uvar, cfg.wvar)), ft_warning('there is an intersection between cfg.uvar and cfg.wvar'); end
if ~isempty(intersect(cfg.uvar, cfg.cvar)), ft_warning('there is an intersection between cfg.uvar and cfg.cvar'); end
if ~isempty(intersect(cfg.wvar, cfg.cvar)), ft_warning('there is an intersection between cfg.wvar and cfg.cvar'); end

fprintf('total number of measurements     = %d\n', Nrepl);
fprintf('total number of variables        = %d\n', Nvar);
fprintf('number of independent variables  = %d\n', length(cfg.ivar));
fprintf('number of unit variables         = %d\n', length(cfg.uvar));
fprintf('number of within-cell variables  = %d\n', length(cfg.wvar));
fprintf('number of control variables      = %d\n', length(cfg.cvar));
fprintf('using a %s resampling approach\n', cfg.resampling);

if ~isempty(cfg.cvar)
  % the different levels of the control variable indicate the blocks in which the resampling can be done
  % the replications should not be resampled over the blocks
  blocklevel = unique(design(cfg.cvar,:)', 'rows')';
  for i=1:size(blocklevel,2)
    blocksel{i} = find(all(design(cfg.cvar,:)==repmat(blocklevel(:,i), 1, Nrepl), 1));
    blocklen(i) = length(blocksel{i});
  end
  for i=1:size(blocklevel,2)
    fprintf('------------------------------------------------------------\n');
    if length(cfg.cvar)>1
      fprintf('resampling the subset where the control variable level is [%s]\n', num2str(blocklevel(:,i)));
    else
      fprintf('resampling the subset where the control variable level is %s\n', num2str(blocklevel(:,i)));
    end
    % use recursion to resample the replications within each block
    tmpcfg = cfg;
    tmpcfg.cvar = [];
    blockres{i} = blocksel{i}(resampledesign(tmpcfg, design(:, blocksel{i})));
  end
  % concatenate the blocks and return the result
  resample = cat(2, blockres{:});
  return
end

if isnumeric(cfg.numrandomization) && cfg.numrandomization==0
  % no randomizations are requested, return an empty shuffling matrix
  resample = zeros(0,Nrepl);
  return;
end

if ~isempty(cfg.wvar)
  % keep the elements within a cell together, e.g. multiple tapers in a trial
  % this is realized by replacing the design matrix temporarily with a smaller version
  blkmeas = unique(design(cfg.wvar,:)', 'rows')';
  for i=1:size(blkmeas,2)
    blksel{i} = find(all(design(cfg.wvar,:)==repmat(blkmeas(:,i), 1, Nrepl), 1));
    blklen(i) = length(blksel{i});
  end
  if any(blklen~=blklen(1))
    ft_error('the number of repetitions per block should be constant');
  end
  for i=1:size(blkmeas,2)
    if any(diff(design(:, blksel{i}), 1, 2)~=0)
      ft_error('the design matrix variables should be constant within a block');
    end
  end
  orig_design = design;
  orig_Nrepl  = Nrepl;
  % replace the design matrix by a blocked version
  design = zeros(size(design,1), size(blkmeas,2));
  Nrepl  = size(blkmeas,2);
  for i=1:size(blkmeas,2)
    design(:,i) = orig_design(:, blksel{i}(1));
  end
end

% do some validity checks
if Nvar==1 && ~isempty(cfg.uvar)
  ft_error('A within-units shuffling requires a at least one unit variable and at least one independent variable');
end

if isempty(cfg.uvar) && strcmp(cfg.resampling, 'permutation')
  % reshuffle the colums of the design matrix
  if ischar(cfg.numrandomization) && strcmp(cfg.numrandomization, 'all')
    % systematically shuffle the columns in the design matrix
    Nperm = prod(1:Nrepl);
    fprintf('creating all possible permutations (%d)\n', Nperm);
    resample = perms(1:Nrepl);
  elseif ~ischar(cfg.numrandomization)
    % randomly shuffle the columns in the design matrix the specified number of times
    resample = zeros(cfg.numrandomization, Nrepl);
    for i=1:cfg.numrandomization
      resample(i,:) = randperm(Nrepl);
    end
  end

elseif isempty(cfg.uvar) && strcmp(cfg.resampling, 'bootstrap')
  % randomly draw with replacement, keeping the number of elements the same in each class
  % only the test under the null-hypothesis (h0) is explicitely implemented here
  % but the h1 test can be achieved using a control variable
  resample = zeros(cfg.numrandomization, Nrepl);
  for i=1:cfg.numrandomization
    resample(i,:) = randsample(1:Nrepl, Nrepl, true);
  end

elseif ~isempty(cfg.uvar) && strcmp(cfg.resampling, 'permutation')
  % reshuffle the colums of the design matrix, keep the rows of the design matrix with the unit variable intact
  unitlevel = unique(design(cfg.uvar,:)', 'rows')';
  for i=1:size(unitlevel,2)
    unitsel{i} = find(all(design(cfg.uvar,:)==repmat(unitlevel(:,i), 1, Nrepl), 1));
    unitlen(i) = length(unitsel{i});
  end
  if length(cfg.uvar)==1
    fprintf('repeated measurement in variable %d over %d levels\n', cfg.uvar, length(unitlevel));
    fprintf('number of repeated measurements in each level is '); fprintf('%d ', unitlen); fprintf('\n');
  else
    fprintf('repeated measurement in mutiple variables over %d levels\n', length(unitlevel));
    fprintf('number of repeated measurements in each level is '); fprintf('%d ', unitlen); fprintf('\n');
  end
  fprintf('the maximum number of unique permutations is %d\n', prod(unitlen));

  if ischar(cfg.numrandomization) && strcmp(cfg.numrandomization, 'all')
    % create all possible permutations by systematic assignment
    if any(unitlen~=2)
      % it would be possible to also implement it for other cases
      % but so far ther has not been a concrete need for it
      ft_error('cfg.numrandomization=''all'' is only supported for two repeated measurements');
    end
    Nperm = 2^(length(unitlevel));
    resample = zeros(Nperm, Nrepl);
    for i=1:Nperm
      flip  = dec2bin( i-1, length(unitlevel));
      for j=1:length(unitlevel)
        if     strcmp('0', flip(j))
          resample(i, unitsel{j}(1)) = unitsel{j}(1);
          resample(i, unitsel{j}(2)) = unitsel{j}(2);
        elseif strcmp('1', flip(j))
          resample(i, unitsel{j}(1)) = unitsel{j}(2);
          resample(i, unitsel{j}(2)) = unitsel{j}(1);
        end
      end
    end
    fprintf('generated all %d possible permutations\n', 2^(length(unitlevel)));

  elseif ~ischar(cfg.numrandomization)
    % see https://github.com/fieldtrip/fieldtrip/issues/1313
    if cfg.numrandomization > prod(unitlen)
      ft_warning('the number of randomizations (%d) is larger than the maximum number of unique permutations (%d), better use cfg.numrandomization=''all''', cfg.numrandomization, prod(unitlen))
    elseif cfg.numrandomization/prod(unitlen) > 0.5
      ft_warning('the number of randomizations (%d) is close to the maximum number of unique permutations (%d), better use cfg.numrandomization=''all''', cfg.numrandomization, prod(unitlen))
    end

    % create the desired number of permutations by random shuffling
    resample = zeros(cfg.numrandomization, Nrepl);
    for i=1:cfg.numrandomization
      for j=1:length(unitlevel)
        resample(i, unitsel{j}) = unitsel{j}(randperm(length(unitsel{j})));
      end
    end
    fprintf('generated %d random permutations\n', cfg.numrandomization);
  end

elseif length(cfg.uvar)==1 && strcmp(cfg.resampling, 'bootstrap') && isempty(cfg.cvar)
  % randomly draw with replacement, keeping the number of elements the same in each class
  % only the test under the null-hypothesis (h0) is explicitely implemented here
  % but the h1 test can be achieved using a control variable

  % FIXME allow for length(cfg.uvar)>1, does it make sense in the first place
  % bootstrap the units of observation
  units = design(cfg.uvar,:);
  for k = 1:length(unique(units))
    sel = find(units==k);
    indx(:,k) = sel;
    Nrep(k)   = length(sel);
  end
  resample = zeros(cfg.numrandomization, Nrepl);

  %sanity check on number of repetitions
  if any(Nrep~=Nrep(1)), ft_error('all units of observation should have an equal number of repetitions'); end

  if max(units(:))<20
    ft_warning('fewer than 20 units warrants explicit checking of double occurrences of ''bootstraps''');
    checkunique = 1;
  else
    checkunique = 0;
  end

  if ~checkunique
    for i=1:cfg.numrandomization
      tmp           = randsample(1:Nrepl/Nrep(1), Nrepl/Nrep(1), true);
      for k=1:size(indx,1)
        resample(i,indx(k,:)) = indx(k,tmp);
      end
    end
  else
    tmp = zeros(cfg.numrandomization*10, Nrepl/Nrep(1));
    for i=1:cfg.numrandomization*10
      tmp(i,:) = sort(randsample(1:Nrepl/Nrep(1), Nrepl/Nrep(1), true));
    end

    tmp = unique(tmp, 'rows');
    fprintf('found %d unique rows in bootstrap matrix of %d bootstraps', size(tmp,1), cfg.numrandomization*10);

    if size(tmp,1)<cfg.numrandomization
      fprintf('using only %d unique bootstraps\n', size(tmp,1));
      cfg.numrandomization = size(tmp,1);
      index = 1:size(tmp,1);
    else
      index = randperm(size(tmp,1));
      index = index(1:cfg.numrandomization);
    end

    tmp = tmp(index,:);
    for i=1:cfg.numrandomization
      for k=1:size(indx,1)
        resample(i,indx(k,:)) = indx(k,tmp(i,:));
      end
    end

  end

else
  ft_error('Unsupported configuration for resampling.');
end

if ~isempty(cfg.wvar)
  % switch back to the original design matrix and expand the resample ordering
  % matrix so that it reorders all elements within a cell together
  Nrepl  = orig_Nrepl;
  expand = zeros(size(resample,1), Nrepl);
  for i=1:size(resample,1)
    expand(i,:) = cat(2, blksel{resample(i,:)});
  end
  resample = expand;
  % replace the blocked version of the design by the original one
  design = orig_design;
  Nrepl  = orig_Nrepl;
end

% in some cases it is possible to reduce the number of permutations by taking only the ones that result in unique condition sequences
% but important is that the relative requencies of the condition sequences remains the same
if strcmp(efficient, 'yes')
  if numel(cfg.ivar)<1
    ft_error('this reqiures at least one independent variable to be specified (ivar)');
  end
  if numel(cfg.uvar)>0
    ft_error('this is not yet supported in combination with a unit of observation (uvar)');
  end

  original = zeros(size(resample,1), numel(cfg.ivar)*size(resample,2));
  for i=1:size(resample,1)
    for j=1:numel(cfg.ivar)
      sel = (1:Nrepl) + (j-1)*Nrepl;
      original(i,sel) = design(cfg.ivar(j), resample(i,:));
    end
  end

  [reduced, indx1, indx2] = unique(original, 'rows');
  % this returns the reduced condition sequences, where
  % indx1 contains for the reduced row number the original (last) row number
  % indx2 contains for the original condition sequences the reduced row number
  frequency = zeros(size(reduced,1),1);
  for i=1:size(reduced,1)
    frequency(i) = sum(indx2==i);
  end

  if all(frequency==frequency(1))
    % the relative frequency of the reduced condition sequences is the same
    % as the relative frequency of the original condition sequences, which
    % means that the reduced set of permutations is appropriate
    fprintf('using the reduced set of permutations (%d) instead of the original permutations (%d)\n', size(reduced,1), size(resample,1));
    resample = resample(indx1,:);
  else
    % don't use the reduced set of permutations
    fprintf('the reduced set has different relative frequencies of the conditions, retaining the original permutations\n')
  end
end % if efficient

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [delete]=PlotInitialStateMatrix3D(ToPlotOrNot)

if ToPlotOrNot==1
    for k=1:size(state,3)
        for j=1:size(state,2)
            for i=1:size(state,1)
                if state(i,j,k)==1
                    plot3(x(i,j,k),y(i,j,k),z(i,j,k),'k.');hold on;
                end
            end
        end
    end
end
axis square;box on

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/34985-monte-carlo-simulation-of-three-dimensional-grain-growth-code-version-no-1-basic/Monte Carlo Simulation Q-state Potts model 3D square-lattice - microstructure/PlotInitialStateMatrix_3D_QPOTTS.m


text:
function vectors2d
%VECTORS2D Description of functions operating on plane vectors.
%
%   A vector is defined by its two cartesian coordinates, put into a row
%   vector of 2 elements:
%   V = [vx vy];
%
%   Several vectors are stored in a matrix with two columns, one for the
%   x-coordinate, one for the y-coordinate.
%   VS = [vx1 vy1 ; vx2 vy2 ; vx3 vy3];
%
%   See also
%   vectorNorm, vectorAngle, isPerpendicular, isParallel
%   normalizeVector, transformVector, rotateVector
%

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

help('vectors2d');

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function linplus_test534 ( )

%*****************************************************************************80
%
%% TEST534 tests R8S3_WRITE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 100;
  nz_num = 3 * n - 2;
  output_file = 'r8s3_matrix.txt';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST534\n' );
  fprintf ( 1, '  For a R8S3 system,\n' );
  fprintf ( 1, '  R8S3_WRITE writes the matrix to a file.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N =         %d\n', n );
  fprintf ( 1, '  Matrix nonzeros NZ_NUM = %d\n', nz_num );

  isym = 0;
%
%  Set the matrix values.
%
  k = 0;
  for i = 1 : n
    k = k + 1;
    j = i;
    row(k) = i;
    col(k) = j;
    a(k) = 100 * i + j;
  end

  for i = 2 : n
    j = i - 1;
    k = k + 1;
    row(k) = i;
    col(k) = j;
    a(k) = 100 * i + j;
  end

  for i = 1 : n-1
    j = i + 1;
    k = k + 1;
    row(k) = i;
    col(k) = j;
    a(k) = 100 * i + j;
  end

  r8s3_print_some ( n, n, nz_num, isym, row, col, a, 1, 1, ...
    10, 10, '  Initial 10x10 block of R8S3 matrix:' );

  r8s3_write ( n, nz_num, isym, row, col, a, output_file );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  R8S3_WRITE wrote the matrix data to "%s"\n.', ...
    output_file );

  return
end


text:
function hash = hlp_cryptohash(data,fromfile)
% Compute an MD5 hash of a file, string or generic data structure.
% Hash = hlp_cryptohash(Data,FromFile)
%
% In:
%   Data : data to be hashed; can be a filename, a string, or any other MATLAB data structure.
%
%   FromFile : if true, data is interpreted as a file name (default: false)
%
% Out:
%   Hash : MD5 hash (decimal, lowercase) of the Data
%
% Examples:
%   % calculate an md5 hash of a file
%   hlp_cryptohash('myscript.m',true);
%
%   % calculate an md5 hash of a data structure
%   hlp_cryptohash(lastmodel);
%
% See also:
%   hlp_fingerprint, hlp_serialize
%
%					      	     Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-10-10

% 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

% maximum amount of memory used inside the Java VM
max_java_memory = 2^26; % 64 MB

if nargin < 2
    fromfile = false; end
if ~(isequal(fromfile,true) || isequal(fromfile,false))
    error('The given FromFile argument must be true or false.'); end

if fromfile
    % take data as a file name
    if ~ischar(data)
        error('To represent a file name, Data should be a string.'); end
    if ~exist(data,'file')
        error('The file %s does not exist.',data); end
    f = fopen(data,'r');
    try
        data = fread(f,Inf);
        fclose(f);
    catch e
        try
            fclose(f);
        catch %#ok<CTCH>
        end
        rethrow(e);
    end
else
    % take data literally
    if ~ischar(data)
        data = hlp_serialize(data); end
end

% use Java to hash the data (idea from Michael Kleder)
hasher = java.security.MessageDigest.getInstance('MD5');
data = uint8(data);
if length(data) <= max_java_memory
    hasher.update(data);
else
    numsplits = ceil(length(data)/max_java_memory);
    for i=0:numsplits-1
        range = 1+floor(i*length(data)/numsplits) : min(length(data),floor((i+1)*length(data)/numsplits));
        hasher.update(data(range));
    end
end

hash = dec2hex(typecast(hasher.digest,'uint8'),2)';
hash = lower(hash(:)');

repo: REST


text:
function [b a]=get_high_shelving_filter(g,Q,f,Fs)

A=10^(g/40);
w=2*pi*f/Fs;
sn=sin(w);
cs=cos(w);
al=sn/(2*Q);
bt=sqrt(A)/Q;

b=A*[(A+1)+(A-1)*cs+bt*sn,...
   -2*((A-1)+(A+1)*cs),...
   (A+1)+(A-1)*cs-bt*sn];

a=[(A+1)-(A-1)*cs+bt*sn,...
   2*((A-1)-(A+1)*cs),...
   (A+1)-(A-1)*cs-bt*sn];

b=b/a(1);
a=a/a(1);


text:
 close all;clear all;clc
tic;
%% inputs

options=odeset('RelTol',1e-9,'AbsTol',1e-9);

solinit = bvpinit(linspace(0,6,3),5*ones(3,1));

[T,Y]=bvpMSM('bvpme',solinit,@bcsme,@dbcsme,options);
figure(2);plot(T,Y);axis tight
toc


text:
function varargout = gallery2(name)
%CHEB.GALLERY2   Chebfun2 example functions.
%   F = CHEB.GALLERY2(NAME) returns a chebfun2 corresponding to NAME. See the
%   listing below for available names.
%
%   For example,  plot(cheb.gallery2('peaks'))  plots the classic MATLAB peaks
%   function represented as a chebfun2. For details of how each function is
%   constructed, try type +cheb/gallery2 or edit cheb.gallery2.
%
%   [F,FA] = CHEB.GALLERY2(NAME) also returns the anonymous function FA used to
%   define the function.
%
%   CHEB.GALLERY2 with no input argument returns a function selected at
%   random from the gallery.
%
%   CHEB.GALLERY2 with no output argument creates a plot of the selected
%   function.
%
%   airyreal     Real part of Airy Ai function
%   airycomplex  Imaginary part of Airy Ai function
%   bump         2D C-infinity function with compact support
%   challenge    Function from SIAM 100-digit challenge
%   peaks        Classic MATLAB peaks function
%   rosenbrock   Challenging test function in optimization
%   roundpeg     Approx characteristic function of a disk, rank 45
%   smokering    A halo, hoop, or hole
%   squarepeg    Approx characteristic function of a square, rank 1
%   tiltedpeg    A tilted variant of squarepeg, rank 100
%   waffle       Function with horizontal and vertical ridges
%
%   Gallery functions are subject to change in future releases of Chebfun.
%
% See also CHEB.GALLERY, CHEB.GALLERYTRIG, CHEB.GALLERY3, CHEB.GALLERYDISK, CHEB.GALLERYSPHERE.

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

% If the user did not supply an input, return a function selected at
% random from the gallery.
if ( nargin == 0 )
    names = {'airyreal', 'airycomplex', 'challenge', 'bump', 'peaks', ...
        'rosenbrock', 'roundpeg', 'smokering', 'squarepeg', 'tiltedpeg', ...
        'waffle'};
    name = names{randi(length(names))};
end

% The main switch statement.
switch lower(name)

    % Real part of Airy Ai function:
    case 'airyreal'
        fa = @(z) real(airy(z));
        f = chebfun2(fa, [-10 10 0 1]);

    % Imaginary part of Airy Ai function:
    case 'airycomplex'
        fa = @(z) airy(z);
        f = chebfun2(fa, [-10 10 -5 5]);

    % Two-dimensional bump function:
    case 'bump'
        fa = @(x,y) (x.^2 + y.^2 < 1).*exp(-(x.^2 + y.^2 < 1)./ ...
            (1 - x.^2 - y.^2));
        f = chebfun2(fa, [-1 1 -1 1]*2);

    % Function from SIAM 100-digit challenge:
    case 'challenge'
        fa = @(x,y) exp(sin(50*x)) + sin(60*exp(y)) + sin(70*sin(x)) + ...
            sin(sin(80*y)) - sin(10*(x+y)) + (x.^2+y.^2)./4;
        f = chebfun2(fa);

    % The classic MATLAB peaks function:
    case 'peaks'
        fa = @(x,y) 3*(1-x).^2.*exp(-x.^2 - (y+1).^2) ...
            - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2 - y.^2) ...
            - 1/3*exp(-(x+1).^2 - y.^2);
        f = chebfun2(fa, [-3 3 -3 3]);

    % A challenging test function in optimization:
    case 'rosenbrock'
        fa = @(x,y) (1-x).^2 + 100*(y-x.^2).^2;
        f = chebfun2(fa, [-2 2 -1 3]);

    % Approximately the characteristic function of a disk:
    case 'roundpeg'
        fa = @(x,y) 1./(1+((2*x).^2+(2*y).^2).^10);
        f = chebfun2(fa);

    % A halo, hoop, or hole:
    case 'smokering'
        fa = @(x,y) exp(-100*(x.^2 - x.*y + 2*y.^2 - 1/2).^2);
        f = chebfun2(fa);

    % Approximately the characteristic function of a square:
    case 'squarepeg'
        fa = @(x,y) 1./((1+(2*x).^20).*(1+(2*y).^20));
        f = chebfun2(fa);

    % A tilted version of squarepeg:
    case 'tiltedpeg'
        fa = @(x,y) 1./((1+(2*x+.4*y).^20).*(1+(2*y-.4*x).^20));
        f = chebfun2(fa);

    % A function with horizontal and vertical ridges:
    case 'waffle'
        fa = @(x,y) 1./(1 + 1e3*((x.^2-.25).^2.*(y.^2-.25).^2));
        f = chebfun2(fa);

    % Raise an error if the input is unknown.
    otherwise
        error('CHEB:GALLERY2:unknown:unknownFunction', ...
            'Unknown function.')
end

% Only return something if there is an output argument.
if ( nargout > 0 )
    varargout = {f, fa};
else
    % Otherwise, plot the function.
    plot(f)
    title([name ', rank = ' num2str(length(f))])
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [M] = spm_DEM_M_set(M)
% sets indices and performs checks on hierarchical models
% FORMAT [M] = spm_DEM_M_set(M)
%
% for each level (i); required fields
%
%   M(i).g  = y(t)  = g(x,v,P)    {inline function, string or m-file}
%   M(i).f  = dx/dt = f(x,v,P)    {inline function, string or m-file}
%
% and
%
%   M(i).m  = number of inputs v(i + 1);
%   M(i).n  = number of states x(i);
%   M(i).l  = number of output v(i);
%
% or
%
%   M(i).x  = hidden states;
%   M(i).v  = causal states;
%
% for each level (i); optional fields
%
%   M(i).pE = prior expectation of p model-parameters
%   M(i).pC = prior covariances of p model-parameters
%   M(i).hE = prior expectation of h log-precision (cause noise)
%   M(i).hC = prior covariances of h log-precision (cause noise)
%   M(i).gE = prior expectation of g log-precision (state noise)
%   M(i).gC = prior covariances of g log-precision (state noise)
%   M(i).xC = prior covariances of states
%   M(i).Q  = precision components (input noise)
%   M(i).R  = precision components (state noise)
%   M(i).V  = fixed precision (input noise)
%   M(i).W  = fixed precision (state noise)
%
%
% sets fields, checks internal consistency of model specification and sets
% estimation parameters.  If a single hyperparameter is supplied i.i.d
% components are assumed (i.e., Q = I, R = I)
%--------------------------------------------------------------------------
%
%   M(1).E.s;     = smoothness (s.d. in time bins)
%   M(1).E.d;     = embedding order q(v)  (i.e., number of derivatives)
%   M(1).E.n;     = embedding order q(x)
%
% If the highest level involves any dynamic or static transformation
% of its inputs a further level is added with flat priors
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_DEM_M_set.m 7322 2018-05-31 09:47:15Z karl $

% order
%--------------------------------------------------------------------------
g      = length(M);

% set missing fields
%==========================================================================

% check for specification of hidden states
%--------------------------------------------------------------------------
if isfield(M,'f') && ~isfield(M,'n') && ~isfield(M,'x')
    msgbox('please specify hidden states or their number')
end

% check supra-ordinate level and add one (with flat priors) if necessary
%--------------------------------------------------------------------------
try
    fcnchk(M(g).g);
    g      = g + 1;
    M(g).l = M(g - 1).m;
end
M(g).m = 0;
M(g).n = 0;

% default fields for static models (hidden states)
%--------------------------------------------------------------------------
if ~isfield(M,'f')
    [M.f] = deal(inline('sparse(0,1)','x','v','P'));
    [M.x] = deal(sparse(0,1));
    [M.n] = deal(0);
end
for i  = 1:g
    try
        fcnchk(M(i).f);
    catch
        M(i).f = inline('sparse(0,1)','x','v','P');
        M(i).x = sparse(0,1);
        M(i).n = 0;
    end
end

% consistency and format check on states, parameters and functions
%==========================================================================

% prior expectation of parameters M.pE
%--------------------------------------------------------------------------
try
    M.pE;
catch
    % Assume fixed parameters
    %----------------------------------------------------------------------
    for i = 1:g
        M(i).pE = sparse(0,0);
    end
end

% and priors covariances - p
%--------------------------------------------------------------------------
try
    M.pC;
catch

    % Assume fixed (zero variance) parameters
    %----------------------------------------------------------------------
    for i = 1:g
        p       = length(spm_vec(M(i).pE));
        M(i).pC = sparse(p,p);
    end
end

% check pC, if user specified
%--------------------------------------------------------------------------
for i = 1:g

    % number of parameters
    %----------------------------------------------------------------------
    np  = length(spm_vec(M(i).pE));

    % Assume fixed parameters if not specified
    %----------------------------------------------------------------------
    if isempty(M(i).pC)
        M(i).pC = sparse(np,np);
    end

    % convert variances to covariances if necessary
    %----------------------------------------------------------------------
    if isvector(M(i).pC) || isstruct(M(i).pC)
        M(i).pC = sparse(diag(spm_vec(M(i).pC)));
    end

    % convert variance to covariances if necessary
    %----------------------------------------------------------------------
    if isscalar(M(i).pC)
        M(i).pC = speye(np,np)*M(i).pC;
    end

    % check size
    %----------------------------------------------------------------------
    if length(M(i).pC) ~= np
        error('please check: M(%i).pC',i)
    end

end

% get inputs
%--------------------------------------------------------------------------
try
    v  = M(g).v;
catch
    v  = sparse(0,0);
end
if isempty(v)
    try
        v = sparse(M(g - 1).m,1);
    end
end
if isempty(v)
    try
        v = sparse(M(g).l,1);
    end
end
M(g).l    = length(spm_vec(v));
M(g).v    = v;

% check functions
%--------------------------------------------------------------------------
for i = (g - 1):-1:1
    try
        x = M(i).x;
    catch
        x = sparse(M(i).n,1);
    end
    if isempty(x) && M(i).n
        x = sparse(M(i).n,1);
    end

    % check f(x,v,P)
    %----------------------------------------------------------------------
    try
        M(i).f  = fcnchk(M(i).f,'x','v','P');
    end
    try
        f       = feval(M(i).f,x,v,M(i).pE);
        if length(spm_vec(x)) ~= length(spm_vec(f))
            error('please check: M(%i).f(x,v,P)',i)
        end

    catch
        sprintf('??? evaluation failure: M(%i).f(x,v,P)',i)
        error(lasterror)
    end
    try M(i).fx = fcnchk(M(i).fx,'x','v','P'); end
    try M(i).fv = fcnchk(M(i).fv,'x','v','P'); end
    try M(i).fp = fcnchk(M(i).fp,'x','v','P'); end

    % check g(x,v,P)
    %----------------------------------------------------------------------
    try
        M(i).g = fcnchk(M(i).g,'x','v','P');
    end
    try
        M(i).m = length(spm_vec(v));
        v      = feval(M(i).g,x,v,M(i).pE);
        M(i).l = length(spm_vec(v));
        M(i).n = length(spm_vec(x));

        M(i).v = v;
        M(i).x = x;

    catch
        sprintf('??? evaluation failure: M(%i).g(x,v,P)',i)
        error(lasterror)
    end
    try M(i).gx = fcnchk(M(i).gx,'x','v','P'); end
    try M(i).gv = fcnchk(M(i).gv,'x','v','P'); end
    try M(i).gp = fcnchk(M(i).gp,'x','v','P'); end

end

% full priors on states
%--------------------------------------------------------------------------
try, M.xP; catch, M(1).xP = []; end
try, M.vP; catch, M(1).vP = []; end
for i = 1:g

    % hidden states
    %----------------------------------------------------------------------
    if isvector(M(i).xP), M(i).xP = diag(M(i).xP); end
    if length(M(i).xP) ~= M(i).n
        try
            M(i).xP = speye(M(i).n,M(i).n)*M(i).xP(1);
        catch
            M(i).xP = sparse(M(i).n,M(i).n);
        end
    end

    % hidden states
    %----------------------------------------------------------------------
    if isvector(M(i).vP), M(i).vP = diag(M(i).vP); end
    if length(M(i).vP) ~= M(i).l
        try
            M(i).vP = speye(M(i).l,M(i).l)*M(i).vP(1);
        catch
            M(i).vP = sparse(M(i).l,M(i).l);
        end
    end
end

% number of x (hidden states)
%--------------------------------------------------------------------------
nx     = sum(spm_vec(M.n));

% Hyperparameters and components (causes: Q V and hidden states R, W)
%==========================================================================
try, M.Q;  catch, M(1).Q  = []; end
try, M.R;  catch, M(1).R  = []; end
try, M.V;  catch, M(1).V  = []; end
try, M.W;  catch, M(1).W  = []; end
try, M.hE; catch, M(1).hE = []; end
try, M.gE; catch, M(1).gE = []; end
try, M.ph; catch, M(1).ph = []; end
try, M.pg; catch, M(1).pg = []; end

% check hyperpriors hE - [log]hyper-parameters and components
%--------------------------------------------------------------------------
pP    = 1;               % prior precision on log-precisions
for i = 1:g

    % make sure components are cell arrays
    %----------------------------------------------------------------------
    if ~isempty(M(i).Q) && ~iscell(M(i).Q), M(i).Q = {M(i).Q}; end
    if ~isempty(M(i).R) && ~iscell(M(i).R), M(i).R = {M(i).R}; end

    % check hyperpriors
    %======================================================================

    % vectorise
    %----------------------------------------------------------------------
    M(i).hE = spm_vec(M(i).hE);
    M(i).gE = spm_vec(M(i).gE);

    % check hyperpriors (expectations)
    %----------------------------------------------------------------------
    if isempty(M(i).hE), M(i).hE = sparse(length(M(i).Q),1); end
    if isempty(M(i).gE), M(i).gE = sparse(length(M(i).R),1); end

    % check hyperpriors (covariances)
    %----------------------------------------------------------------------
    try, M(i).hC*M(i).hE; catch, M(i).hC = speye(length(M(i).hE))/pP; end
    try, M(i).gC*M(i).gE; catch, M(i).gC = speye(length(M(i).gE))/pP; end

    if isempty(M(i).hC), M(i).hC = speye(length(M(i).hE))/pP; end
    if isempty(M(i).gC), M(i).gC = speye(length(M(i).gE))/pP; end

    % check Q and R (precision components)
    %======================================================================

    % check components and assume i.i.d if not specified
    %----------------------------------------------------------------------
    if length(M(i).Q) > length(M(i).hE)
        M(i).hE = sparse(length(M(i).Q),1) + M(i).hE(1);
    end
    if length(M(i).Q) < length(M(i).hE)
        M(i).Q  = {speye(M(i).l,M(i).l)};
        M(i).hE = M(i).hE(1);
    end
    if length(M(i).hE) > length(M(i).hC)
        M(i).hC = speye(length(M(i).Q))*M(i).hC(1);
    end
    if length(M(i).R) > length(M(i).gE)
        M(i).gE = sparse(length(M(i).R),1) + M(i).gE(1);
    end
    if length(M(i).R) < length(M(i).gE)
        M(i).R  = {speye(M(i).n,M(i).n)};
        M(i).gE = M(i).gE(1);
    end
    if length(M(i).gE) > length(M(i).gC)
        M(i).gC = speye(length(M(i).R))*M(i).gC(1);
    end

    % check consistency and sizes (Q)
    %----------------------------------------------------------------------
    for j = 1:length(M(i).Q)
        if length(M(i).Q{j}) ~= M(i).l
            error('wrong size; M(%d).Q{%d}',i,j)
        end
    end

    % check consistency and sizes (R)
    %----------------------------------------------------------------------
    for j = 1:length(M(i).R)
        if length(M(i).R{j}) ~= M(i).n
            error('wrong size; M(%d).R{%d}',i,j)
        end
    end

    % check V and W (lower bound on precisions)
    %======================================================================

    % check V and assume unit precision if improperly specified
    %----------------------------------------------------------------------
    if isvector(M(i).V), M(i).V = diag(M(i).V); end
    if length(M(i).V) ~= M(i).l
        try
            M(i).V = speye(M(i).l,M(i).l)*M(i).V(1);
        catch
            if isempty(M(i).hE) && isempty(M(i).ph)
                M(i).V = speye(M(i).l,M(i).l);
            else
                M(i).V = sparse(M(i).l,M(i).l);
            end
        end
    end

    % check W and assume unit precision if improperly specified
    %----------------------------------------------------------------------
    if isvector(M(i).W), M(i).W = diag(M(i).W); end
    if length(M(i).W) ~= M(i).n
        try
            M(i).W = speye(M(i).n,M(i).n)*M(i).W(1);
        catch
            if isempty(M(i).gE) && isempty(M(i).pg)
                M(i).W = speye(M(i).n,M(i).n);
            else
                M(i).W = sparse(M(i).n,M(i).n);
            end
        end
    end
end

% estimation parameters M(1).E.s, n,...
%==========================================================================
% E.s;                               % smoothness (seconds)
% E.dt;                              % time step
% E.d;                               % approximation order of q(x,v)
% E.n;                               % order of embedding (n >= d)

% temporal smoothness - s.d. of kernel
%--------------------------------------------------------------------------
try M(1).E.s;  catch, if nx, M(1).E.s = 1/2; else M(1).E.s = 0; end, end

% time step
%--------------------------------------------------------------------------
try M(1).E.dt; catch, M(1).E.dt = 1; end

% embedding orders
%--------------------------------------------------------------------------
try M(1).E.d;  catch, if nx, M(1).E.d = 2; else M(1).E.d = 0;  end, end
try M(1).E.n;  catch, if nx, M(1).E.n = 6; else M(1).E.n = 0;  end, end

M(1).E.d = min(M(1).E.d,M(1).E.n);

% number of iterations
%--------------------------------------------------------------------------
try M(1).E.nD; catch, if nx, M(1).E.nD = 1; else M(1).E.nD = 8; end, end
try M(1).E.nE; catch,        M(1).E.nE = 8; end
try M(1).E.nM; catch,        M(1).E.nM = 8; end
try M(1).E.nN; catch,        M(1).E.nN = 8; end

% checks on smoothness hyperparameter
%==========================================================================
try, M = rmfield(M,'sv'); end
try, M = rmfield(M,'sw'); end
for i = 1:g

    try, M(i).sv; catch, M(i).sv = M(1).E.s;   end
    try, M(i).sw; catch, M(i).sw = M(1).E.s;   end

    if ~isscalar(M(i).sv), M(i).sv = M(1).E.s; end
    if ~isscalar(M(i).sw), M(i).sw = M(1).E.s; end
end

% check on linear approximation scheme
%==========================================================================
try
    M(1).E.linear;
catch
    M(1).E.linear = 0;
end

% checks on estimability
%==========================================================================

% check that there are informative priors on the states or the causes
%--------------------------------------------------------------------------
Q     = ~norm(M(end).V,1);
for i = 1:(g - 1)
    P = norm(M(i).pC,1) > exp(8);
    if P && Q
        warndlg('please use informative priors on causes or parameters')
    end
end

author: spm


text:
function [F,C] = bfs_orient(F);
  % BFS_ORIENT Consistently orient faces in orientable patches using BFS
  %
  % F = bfs_orient(F);
  %
  % Inputs:
  %  F  #F by 3 list of faces
  % Outputs:
  %  F  #F by 3 list of faces
  %  C  #F list of component ids
  %
  % Example:
  %  RF = F;
  %  I = rand(size(F,1),1)<0.5;
  %  RF(I,:) = fliplr(RF(I,:));
  %  [FF,C] = bfs_orient(RF);
  %  FF = orient_outward(V,F,C);
  %  I = randperm(max(C))';
  %  meshplot(V,RF,'ScalarFieldF',I(C));
  %  meshplot(V,FF,'ScalarFieldF',I(C));
  %
  % See also: orient_outward, manifold_patches
  %

  [C,A] = manifold_patches(F);
  % No self matches
  A = A-diag(diag(A));
  % loop over components
  for c = 1:max(C)
    F(C==c,:) = bfs_orient_patch(F(C==c,:),A(C==c,C==c));
  end

  function FF = bfs_orient_patch(FF,AA)
    % short circuit if already oriented
    E = [FF(:,[2 3]); FF(:,[3 1]); FF(:,[1 2])];
    % Direct all edges so sortE(:,1) < sortE(:,2)
    sortE = sort(E,2);
    OA = sparse(sortE(:,1),sortE(:,2),1-2*(E(:,1)<E(:,2)));
    DA = sparse(sortE(:,1),sortE(:,2),1);
    % Don't count boundary edges (where "redirected" edge only occured once).
    OA(DA==1) = 0;
    if nnz(OA) == 0
      return;
    end

    m = size(FF,1);
    seen = false(m,1);

    es = [2 3;3 1;1 2];
    Q = 1;
    while ~isempty(Q)
      % pop
      f = Q(1);
      Q = Q(2:end);
      if seen(f)
        continue;
      end
      seen(f) = true;
      % only consider manifold-edge neighbors
      neighbors = find(AA(:,f));
      % number of neighbors
      nn = numel(neighbors);
      % loop over edges
      for ei = 1:3
        % this directed edge
        e = FF(f,es(ei,:));
        % All directed edges of neighbors
        E = [FF(neighbors,[2 3]);FF(neighbors,[3 1]);FF(neighbors,[1 2])];
        % Neighbors with matching directed edge ...
        match = neighbors(mod(find(all(bsxfun(@eq,e,E),2))-1,nn)+1);
        % ... must be flipped
        FF(match,:) = fliplr(FF(match,:));
        % Append neighbors to queue
        Q = [Q;neighbors];
      end

    end
    assert(all(seen));
  end

end

author: alecjacobson

repo: gptoolbox


text:
function payload = aip_detector_concatenate(plA, plB)

    if( ~isempty(plA) )
        if isfield(plA,'result')
            plA = plA.result;
        end
    end

    if( ~isempty(plB) )
        if isfield(plB,'result')
            plB = plB.result;
        end
    end

    payload = ConcatenateQRSdetectionPayloads([], plA, plB);

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [ temp ] = gbSegment( img, sigma, k, minSz )
%GBSEGMENT Summary of this function goes here
%   Detailed explanation goes here
[height, width, ~] = size(img);
img_smooth = smooth(img, sigma);
edges = zeros(width*height*4,3);
[gridX, gridY] = meshgrid(1:width,1:height);
num = 0;

x = gridX(:);
y = gridY(:);
vector = [1 0; 0 1; 1 1; 1 -1];
inda = sub2ind([height width], y, x);

for vid = 1:size(vector,1)
    xv = x+vector(vid,1);
    yv = y+vector(vid,2);
    valid = xv>=1 & xv<=width & yv>=1 & yv<=height;
    indav = inda(valid);
    indbv = sub2ind([height width], yv(valid), xv(valid));
    diff = (img_smooth(indav)-img_smooth(indbv)).^2 ...
         + (img_smooth(indav+width*height)-img_smooth(indbv+width*height)).^2 ...
         + (img_smooth(indav+2*width*height)-img_smooth(indbv+2*width*height)).^2;
    edges(num+1:num+sum(valid),:) = [x(valid)-1 + (y(valid)-1)*width ...
                                     xv(valid)-1 + (yv(valid)-1)*width ...
                                     sqrt(diff)];
    num = num + sum(valid);
end

edges = edges';
segment = segmentGraphMex(width, height, num, edges, k, minSz);
L = unique(segment);

temp = zeros(size(segment));
for i = 1:length(L)
    temp(segment==L(i)) = i;
end

end

author: zouchuhang

repo: LayoutNet


text:
function [lineStyle] = getLineStyleStrFromText(lineStyleText)
%getLineStyleFromString Summary of this function goes here
%   Detailed explanation goes here

    switch lineStyleText
        case 'Solid Line'
            lineStyle = '-';
        case 'Dashed Line'
            lineStyle = '--';
        case 'Dotted Line'
            lineStyle = ':';
        case 'Dashed-dot Line'
            lineStyle = '-.';
        otherwise
            lineStyle = '-';
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function DEM_FEP_Least_Action
%--------------------------------------------------------------------------
% This routine uses a Lorenz system to show that the most likely autonomous
% path (or deterministic path) is uniquely identified by its initial and
% final states. In other words, if we knew the end state of an autonomous
% trajectory, then we would implicitly know the path taken from the initial
% particular state, even if we did not know the external states. This point
% is demonstrated using numerical analyses of the Lorenz system; treating
% the first state and an active state and the third as an external state.
% In this example, 1024 solutions are obtained from the same initial
% particular (i.e., sensory and active) states but sampling from a Gaussian
% distribution over external states. The ensuing trajectories over 128 time
% bins of 1/128 seconds are shown in the left panels. The sample
% distribution over active states is shown as a (scaled) histogram along
% the x-axis. Paths that end within 1/8 of an arbitrary active state (here,
% ?? = -4) are shown in red. The corresponding autonomous (i.e., active)
% paths are shown as a function of time in the right panels. one can repeat
% this analysis for different levels of random fluctuations; e.g.,log
% precisions of 2 and 16. The key thing to observe is that as the amplitude
% of random fluctuations decreases (i.e., its precision increases) the
% paths that begin and end in the same place collapse to a single
% trajectory of least action. This is the most likely or deterministic
% path. Clearly, this behaviour rests upon a diffeomorphic mapping between
% the initial and final states: for example, a final active state of -8 has
% the least two paths of least action (xT in the code below).
%__________________________________________________________________________
% Copyright (C) 2005 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEM_FEP_Least_Action.m 7512 2019-01-05 21:15:16Z karl $

% generative model
%==========================================================================                       % switch for demo
spm_figure('GetWin','DEM'); clf

% flow
%--------------------------------------------------------------------------
dt    = 1/128;
G.f   = @(x,v,P,G) v(:) + [-P(1) P(1) 0; P(3) -1 -x(1); x(2) 0 P(2)]*x*dt;

% set up
%--------------------------------------------------------------------------
T     = 128;                        % length of trajectory
N     = 1024;                       % number of paths
g     = exp(-8);                    % amplitude of random fluctuations
s     = exp(2);                     % amplitude deviations
P     = [10; -8/3; 28];             % Rayleigh parameter
x0    = [1; 1; 25];
for k = 1:N

    % random fluctuations (and action)
    %----------------------------------------------------------------------
    U.u      = randn(T,3)*g;
    A(k)     = (1/2)*sum(U.u(:).^2)/(g^2)/T;

    % random deviations in external states (and probability)
    %----------------------------------------------------------------------
    dx       = randn*s;
    A(k)     = A(k) + (1/2)*sum(dx.^2)/(s^2);

    % integrate timeseries with random initial hidden state
    %----------------------------------------------------------------------
    G.x      = x0;
    G.x(3)   = x0(3) + dx;
    x(:,:,k) = spm_int_L(P,G,U);

end

% plot ensemble densities
%--------------------------------------------------------------------------
subplot(2,2,1),cla, hold on
td    = fix(linspace(1,N,32));
for t = td
    plot(x(:,1,t),x(:,3,t),':','MarkerSize',1)
    plot(x(T,1,t),x(T,3,t),'r.','MarkerSize',8)
end
axis square, axis([-20 20 0 60])
title('Paths from initial state','FontSize',16)
xlabel('active state'),ylabel('external state'),drawnow

% find trajectories that start and end at the same place
%==========================================================================

% get surprisal of final (active) state using sample density
%--------------------------------------------------------------------------
nb    = 32;
xN    = squeeze(x(end,1,:));
[n,a] = hist(xN,nb);
n     = 2*nb*n/sum(n);
bar(a,n,1)

% identify and plot trajectories with the same endpoints
%--------------------------------------------------------------------------
xT    = -4;
k     = find(abs(xN - xT) < 1/8);
for t = k(:)'
    plot(x(:,1,t),x(:,3,t),'r')
end
plot([x0(1),x0(1)],[0 48],'r-.')
plot([xT   ,xT   ],[0 48],'r-.')

% paths as a function of time
%--------------------------------------------------------------------------
subplot(2,2,2), hold on
pst   = (1:T)*dt;
for t = k(:)'
    plot(pst,x(:,1,t),'r')
end
title('Autonomous paths','FontSize',16)
xlabel('time (seconds)'),ylabel('active state'),axis square

return

% NB: numerical analysis of the action and surprisal
%--------------------------------------------------------------------------
for i = 1:N
    [d,j] = min(abs(xn - xN(i)));
    L(i)  = -log(n(j));
end
L      = log(spm_softmax(L(:)));
A      = log(spm_softmax(A(:)));

subplot(2,2,2)
plot(L,A,'.',L,L,':')
title('Action and surprisal','FontSize',16)
xlabel('Surprisal of final state'),ylabel('Action of path'),axis square


text:
function dres = build_graph(dres)
ov_thresh = 0.3;  %original = 0.5

for fr = 2:max(dres.fr)
  f1 = find(dres.fr == fr);     %% indices for detections on this frame
  f2 = find(dres.fr == fr-1);   %% indices for detections on the previous frame
  if(isfield(dres, 'ft'))
      feat1 = dres.ft(f1,:);
      feat2 = dres.ft(f2,:);
  end
  for i = 1:length(f1)
    ovs1  = calc_overlap(dres, f1(i), dres, f2);
    if(isfield(dres, 'ft'))
        dist = pdist2(feat1(i,:),feat2);
        inds1 = find(ovs1 > ov_thresh & dist < 5);                       %% find overlapping bounding boxes.
    else
        inds1 = find(ovs1 > ov_thresh);
    end

    ratio1 = dres.h(f1(i))./dres.h(f2(inds1));
    inds2  = (min(ratio1, 1./ratio1) > 0.8);          %% we ignore transitions with large change in the size of bounding boxes.

    dres.nei(f1(i),1).inds  = f2(inds1(inds2))';      %% each detction window will have a list of indices pointing to its neighbors in the previous frame.
  end
end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function [Gx,Gy] = gradient2( I )
% Compute numerical gradients along x and y directions.
%
% For 2D arrays identical to Matlab's gradient() with a spacing value of
% h=1 but ~10-20x faster (due to mexed implementation). Like gradient(),
% computes centered derivatives in interior of image and uncentered
% derivatives along boundaries. For 3D arrays computes x and y gradient
% separately for each channel and concatenates the results.
%
% This code requires SSE2 to compile and run (most modern Intel and AMD
% processors support SSE2). Please see: http://en.wikipedia.org/wiki/SSE2.
%
% USAGE
%  [Gx,Gy] = gradient2( I )
%
% INPUTS
%  I      - [hxwxk] input k channel single image
%
% OUTPUTS
%  Gx     - [hxwxk] x-gradient (horizontal)
%  Gy     - [hxwxk] y-gradient (vertical)
%
% EXAMPLE
%  I=single(imread('peppers.png'))/255;
%  tic, [Gx1,Gy1]=gradient(I,1); toc
%  tic, [Gx2,Gy2]=gradient2(I); toc
%  isequal(Gx1,Gx2), isequal(Gy1,Gy2)
%
% See also gradient, gradientMag
%
% Piotr's Image&Video Toolbox      Version 3.00
% Copyright 2012 Piotr Dollar & Ron Appel.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

[Gx,Gy]=gradientMex('gradient2',I);


text:
classdef SumOverDim < dagnn.ElementWise
  % @author: Hakan Bilen
  % SumOverDim is the sum of the elements of inputs{1} over dimension dim
  properties
    dim = 3;
  end

  methods
    function outputs = forward(obj, inputs, params)
      outputs{1} = sum(inputs{1},obj.dim) ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)

      ndims = ones(1,numel(size(inputs{1})));
      ndims(obj.dim) = size(inputs{1},obj.dim);
      derInputs{1} = repmat(derOutputs{1},ndims);

      derParams = {} ;
    end

    function outputSizes = getOutputSizes(obj, inputSizes)
      outputSizes{1} = inputSizes{1} ;
      outputSizes{1}(obj.dim) = 1;
    end

    function obj = SumOverDim(varargin)
      obj.load(varargin) ;
      obj.dim = obj.dim;
    end
  end
end

author: hbilen

repo: WSDDN

save_path: github-repos/MATLAB/hbilen-WSDDN


text:
function [m,v]=v_chimv(n,l,s)
%V_CHIMV approximate mean and variance of non-central chi distribution [m,v]=(n,l,s)
%
%  Inputs:	n = degrees of freedom
%           l = non-centrality parameter = sqrt(sum(mean^2)) [default 0]
%               (can be a vector or matrix to calculate many different values at once)
%           s = standard deviation of Gaussian [default 1]
%
% Outputs:  m = mean of chi distribution
%           v = variance of chi distribution
%
% If x=c+randn(n,1) is a column vector of Gaussian random numbers with mean vector c, then
% z=sqrt(x'*x) has a chi distributon with n degrees of freedom and non-centrality parameter
% l=sqrt(c'*c). The mean and variance of a chi distribution are given precisely by
%
%     m = sqrt(2)*exp(gammaln(0.5*n+0.5)-gammaln(0.5*n))*hypergeom(-0.5,0.5*n,-0.5*l^2)
%       = sqrt(pi/2) L(0.5,0.5*n-1,-0.5*l^2)
%     v = n+l^2-m^2
%
% where L(n,a,x) is the generalized Laguerre polynomial L_n^{(a)}(x) but this is very slow
% to calculate so this routine approximates these expressions.
%
% For n=1, the accuracy is high; for n>1, accuracy improves with increasing n.
% Accuracy is worst when the non-centrality parameter, l, is close to s*sqrt(n).
% Worst case errors as a function of n are:
%                       n:    1       2      3       5      10
%   worst case error in m:  1e-15   0.007  0.004  0.0015  0.0005

%      Copyright (C) Mike Brookes 2014
%      Version: $Id: v_chimv.m 4969 2014-08-05 18:24:30Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent ab pp qq nab
if isempty(ab)
    nab=200; % cache a few low values of n
    pp=[ 0.595336298258636  -1.213013700592756  -0.018016200037799   1.999986150447582 0];
    qq=[ -0.161514114798972   0.368983655790737  -0.136992134476950  -0.499681107630725 2];
    ni=1./(1:nab);
    ab=[polyval(qq,ni);polyval(pp,ni)];
end
if nargin<3
    s=1;
    if nargin<2
        l=0;
    end
end
ls=l/s;
l2=(ls).^2;
s2=s^2;
if n<=nab
    if n==1
        m=l.*(1-2*normcdf(-ls))+2*s*normpdf(-ls);
    else
        m=sqrt(l2+n-1+(ab(1,n)+ab(2,n)*l2).^(-1))*s;
    end
else
    m=sqrt(l2+n-1+(polyval(qq,1/n)+polyval(pp,1/n)*l2).^(-1))*s;
end
v=(n+l2)*s2-m.^2;

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD