text:
function [dat] = ft_preproc_derivative(dat, order)

% FT_PREPROC_DERIVATIVE computes the temporal Nth order derivative of the
% data
%
% Use as
%   [dat] = ft_preproc_derivative(dat, order)
% where
%   dat        data matrix (Nchans X Ntime)
%   order      number representing the Nth derivative (default = 1)
%
% See also PREPROC

% 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$

% set the defaults if options are not specified
if nargin<2 || isempty(order)
  order = 1;
end

% preprocessing fails on channels that contain NaN
if any(isnan(dat(:)))
  ft_warning('FieldTrip:dataContainsNaN', 'data contains NaN values');
end

% compute the derivative
for i=1:order
  dat = gradient(dat);
end


text:
function out = ftimes(a,b)
    out = bsxfun(@times,a,b);
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function a = r8mat_geinverse ( a, n, pivot )

%*****************************************************************************80
%
%% R8MAT_GEINVERSE computes the inverse of a matrix factored by R8MAT_GEFA.
%
%  Discussion:
%
%    R8MAT_GEINVERSE is a modified version of the LINPACK routine DGEDI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input/output, real A(N,N).
%    On input, the factor information computed by R8MAT_GEFA.
%    On output, the inverse matrix.
%
%    Input, integer N, the order of the matrix A.
%
%    Input, integer PIVOT(N), the pivot vector from R8MAT_GEFA.
%

%
%  Compute Inverse(U).
%
  for k = 1 : n

    a(k,k) = 1.0 / a(k,k);

    a(1:k-1,k) = - a(1:k-1,k) * a(k,k);

    for j = k + 1 : n

      temp = a(k,j);
      a(k,j) = 0.0;

      a(1:k,j) = a(1:k,j) + temp * a(1:k,k);

    end

  end
%
%  Form Inverse(U) * Inverse(L).
%
  for k = n - 1 : -1 : 1

    work(k+1:n) = a(k+1:n,k);
    a(k+1:n,k) = 0.0;

    for j = k + 1 : n
      a(1:n,k) = a(1:n,k) + work(j) * a(1:n,j);
    end

    if ( pivot(k) ~= k )

      for i = 1 : n
        temp = a(i,k);
        a(i,k) = a(i,pivot(k));
        a(i,pivot(k)) = temp;
      end

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Copyright (C) Daphne Koller, Stanford Univerity, 2012
%
% Author: Huayan Wang

function VisualizeModels(P, G)
K = length(P.c);

f = figure;
while(1)
    for k=1:K
        subplot(1,K,k);
        if size(G,3) == 1  % same graph structure for all classes

            pose = SamplePose(P,G,k);

        else  % different graph structure for each class

            pose = SamplePose(P,G(:,:,k),k);

        end

        img = ShowPose(pose);
        imshow(img);
        pause(0.3)
        if (~ishandle(f)) return; end; % quit loop when user closes the figure
    end
end

author: anhncs

repo: Probabilistic-Graphical-Models


text:
function flag=deletemeshfile(fname)
%
% flag=deletemeshfile(fname)
%
% delete a given work mesh file under the working directory
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%     fname: specified file name (without path)
%
% output:
%     flag: not used
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

try
    if(~isempty(dir(fname)))
        delete(fname);
    end
catch
    error(['You do not have permission to delete temporary files. If you are working in a multi-user ',...
         'environment, such as Unix/Linux and there are other users using iso2mesh, ',...
         'you may need to define ISO2MESH_SESSION=''yourstring'' to make your output ',...
         'files different from others; if you do not have permission to ',mwpath(''),...
         ' as the temporary directory, you have to define ISO2MESH_TEMP=''/path/you/have/write/permission'' ',...
         'in matlab/octave base workspace.']);
end

author: fangq


text:
function [os,os1] = getosmatrix_bb(boxes,gts)
%% Given two sets of bounding boxes, N1 in the first one, and N2 in
%% the second one, compute a N1xN2 overlap score matrix where the
%% overlap score is the ratio of the intersection to union
%% Tomasz Malisiewicz (tomasz@cmu.edu)

if ~exist('gts','var')
  gts = boxes;
end

if numel(boxes) == 0 || numel(gts) == 0
  os = zeros(size(boxes,1),size(gts,1));

  return;
end

x1 = boxes(:,1);
y1 = boxes(:,2);
x2 = boxes(:,3);
y2 = boxes(:,4);

area = (x2-x1+1) .* (y2-y1+1);

xa1 = gts(:,1);
ya1 = gts(:,2);
xa2 = gts(:,3);
ya2 = gts(:,4);

area2 = (xa2-xa1+1) .* (ya2-ya1+1);

os = zeros(size(boxes,1),size(gts,1));
if size(boxes,1)*size(gts,1)==0
  return;
end

if nargout == 2
  os1 = os;
end

for i = 1:size(boxes,1)

  xx1 = max(boxes(i,1),gts(:,1));
  yy1 = max(boxes(i,2),gts(:,2));
  xx2 = min(boxes(i,3),gts(:,3));
  yy2 = min(boxes(i,4),gts(:,4));

  w = xx2-xx1+1;
  h = yy2-yy1+1;

  o = w.*h;
  o( (w<0) | (h<0) ) = 0;

  os(i,:) = o ./ (eps + (area(i) + area2 - o));

  if nargout == 2
    os1(i,:) = (o ./ (eps + area2));
    %os1(i,:) = o ./ area(i);
  end
end

author: quantombone

repo: exemplarsvm


text:
function [varargout] = cubic_cubic_integrated_distance( ...
    x, ...
    y, ...
    g_C, ...
    h_C, ...
    P, ...
    g_D, ...
    h_D, ...
    Q)
  %
  % [E] = cubic_cubic_integrated_distance( x, y, g_C, h_C, P, g_D, h_D, Q)
  % [H,F,c] = cubic_cubic_integrated_distance( x, y, g_C, h_C, P, g_D, h_D)
  %
  % E = ½ ∫ₓʸ ‖ C( g_C u + h_C ) - D( g_D u + h_D ) ‖² du
  %
  % where C's control points are in rows of P and D's control points are in rows
  % of Q.
  %
  % E = 0.5*trace(Q.'*H*Q) + trace(Q.'*F) + c;
  %
  % That is,
  %   Q★ = argmin_Q E(Q)
  %   Q★ = H⁻¹ F
  %
  % Where
  %
  % C(T) = (1-T)³ Pᵢ⁰ + 3(1-T)²T Pᵢ¹ + 3(1-T)T² Pᵢ² + T³ Pᵢ³
  % D(T) = (1-T)³ Qᵢ⁰ + 3(1-T)²T Qᵢ¹ + 3(1-T)T² Qᵢ² + T³ Qᵢ³
  %
  % Inputs:
  %   x  starting value of integral 0≤x<1
  %   y  ending value of integral x<y≤1
  %   g_C  scalar multiplicative factor in parametric values of fixed curve
  %   h_C  scalar additive factor in parametric values of fixed curve
  %   P  4 by dim list of fixed curve control points
  %   g_D  scalar multiplicative factor in parametric values of fixed curve
  %   h_D  scalar additive factor in parametric values of fixed curve
  %   Q  4 by dim list of unknown curve control points (optional)
  % Outputs
  %   E  integrated energy
  %   H  4 by 4 quadratic form of energy as function of Q
  %   F  4 by dim linear term of energy as function of Q
  %   c  scalar constant of energy as function of Q
  %
  if isfloat(x) && isfloat(y)
    assert(x>=0);
    assert(x<1);
    assert(y>=x);
    assert(y<=1);
  end
  dim = size(P,2);

  % https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature
  % n is number of points
  % 2n-1 order polynomial is integrated exactly
  % ↓
  % n = 4

  % https://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules/quadrature_rules.html
  % n=4 Gauss-Legendre quadrature on [-1,1]
  w = [
    0.347854845137453857373063949222
    0.652145154862546142626936050778
    0.652145154862546142626936050778
    0.347854845137453857373063949222
    ];
  a = [
   -0.861136311594052575223946488893
   -0.339981043584856264802665759103
    0.339981043584856264802665759103
    0.861136311594052575223946488893
  ];
  % Move to [x,y] range.
  u = (0.5*a + 0.5)*(y-x)+x;
  % Adjust weights.
  w = ((y-x)/(1- -1))*w;

  % Evaluate C at quadrature points
  C = cubic_eval( P, g_C*u + h_C);

  % Map quadrature points to D's direct paramter space
  T = g_D * u + h_D;

  if nargout<=1
    % Compute energy directly. This is faster. Better be the same as below.
    %
    % It's ever so slightly different (between 1e-16 and 1e-20)
    E = 0.5*sum(w.*(C - cubic_eval( Q, T)).^2,'all');
    varargout{1} = E;
    return;
  end

  % Stack bezier basis matrices for each evaluation point
  M = [(1-T).^3 3*T.*(1-T).^2 3*T.^2.*(1-T) T.^3];

  W = diag(w);

  H = M.'*W*M;
  assert(all(size(H) == [4 4]));
  F = -M.'*W*C;
  assert(size(F,1) == 4);
  assert(size(F,2) == dim);

  c = 0.5*trace(C.'*W*C);

  if isempty(Q)
    warning("Cant compute energy on empty Q");
    E = [];
  else
    % Actually compute energy
    E = 0.5*trace(Q.'*H*Q) + trace(Q.'*F) + c;
  end

  varargout{1} = H;
  varargout{2} = F;
  varargout{3} = c;
  varargout{4} = E;
end

author: alecjacobson

repo: gptoolbox


text:
function cleanDataType(dataTypeName, query)
%
%  cleanDataType(dataTypeName, [query])
%
% Deletes:
%   Inplane/dataType/*.mat
%   Inplane/dataType/TSeries/Scan*/*.mat
%   Likewise for Volume, Gray, and Flat*
%
% If you change this function make parallel changes in:
%     cleanGray, cleanFlat
%
% djh, 2/2001
% We could remove the directory for that data type, too.  Maybe we should.
% ras, 03/07: agreed. now nukes everything in the deleted data type dir.
global HOMEDIR
if ~exist('query','var')
    query = 0;
end
if query
	q = sprintf('Completely delete data type %s?', dataTypeName);
    confirm = questdlg(q, mfilename);
else
    confirm = 'Yes';
end
if strcmp(confirm,'Yes')
    disp(['Deleting dataType: ',dataTypeName]);
    [nDirs,dirList] = countDirs(fullfile(HOMEDIR,'Flat*'));
    dirList{nDirs+1} = 'Inplane';
    dirList{nDirs+2} = 'Gray';
    dirList{nDirs+3} = 'Volume';
    for d = 1:nDirs+3
        datadir = fullfile(HOMEDIR,dirList{d},dataTypeName);
        if exist(datadir,'dir')
            if isunix
                unix(sprintf('rm -r %s', datadir))
			else
				try
					rmdir(datadir, 's');
				catch
					disp('Couldn''t remove the directory.')
				end
            end
        end
    end
end
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function y = meanfilter(x,r,dir)

if ~exist('dir','var'),
  if size(x,1) == 1 && size(x,2) > size(x,1),
    dir = 2;
  else
    dir = 1;
  end
end
sz = size(x);
notdir = true(1,length(sz)); notdir(dir) = false;
notdiridx = find(notdir);
x = permute(x,[dir,notdiridx]);
m = prod(sz(notdir));
x = reshape(x,[sz(dir),m]);

off1 = floor(r/2);
off2 = r - off1 - 1;
y = imfilter(x,ones(r,1),'same',0);
y(1:off1,:) = y(1:off1)./repmat(2*(1:off1),[m,1]);
y(end-off2+1:end,:) = y(end-off2+1:end,:) ./ repmat(2*(off2:-1:1),[m,1]);
y(off1+1:end-off2,:) = y(off1+1:end-off2,:) / r;

y = reshape(y,[sz(dir),sz(notdir)]);
y = permute(y,[2:dir,1,dir+1:length(sz)]);

author: kristinbranson

repo: JAABA


text:
% sine wave fitting from noisy sine signal
% phase fitting has to be considered
% with a fixed omega!
% Written by Dr Chen YangQuan <yqchen@ieee.org>
% Last modified in 05-11-2000
% DESCRIPTIONS:
% s0: sampled series. 1xNp (note here)
% omega: known freq. (2*pi*f) (rad/sec.)
% Ts: sampling period (in Sec.)
% t0: initial time (in Sec.)
% Ahat: estimated amplitude
% Theta: fitted theta_0 (in rad.)
% RMS: root mean squares.
%
% See also "jomega"
function [Ahat,Theta,RMS]=sinefit2(s0,omega,t0,Ts)
Np=size(s0);
t=t0+[0:Np(2)-1]*Ts;
A11= (sin(omega*t)*(sin(omega*t))');
A12= (sin(omega*t)*(cos(omega*t))');
A22= (cos(omega*t)*(cos(omega*t))');
b1=s0*(sin(omega*t))';
b2=s0*(cos(omega*t))';
A=[A11,A12;A12,A22];
Alpha=inv(A)*[b1,b2]';
% be careful here...
Asintheta=Alpha(2);Acostheta=Alpha(1);
Ahat=sqrt(Asintheta*Asintheta+Acostheta*Acostheta);
Theta=atan2(Asintheta,Acostheta);
RMS=sqrt((s0-Ahat*sin(omega*t+Theta))*...
(s0-Ahat*sin(omega*t+Theta))'/(Np(2)-1.));
if (0)
figure;
plot(t,s0,'k:',t, Ahat*sin(omega*t+Theta),'-r');
legend(['measured (RMS=',num2str(RMS),', f=',...
num2str(omega/2/pi),')'],['fitted: ',...
' A.hat=',num2str(Ahat),' | theta.hat=',num2str(Theta*180/pi)])
end
return;

author: Sable

repo: mcbench-benchmarks


text:
function neg_tmpl =  sample_neg(frame, p, out_size, margin)
    %% Sample 16 negative image samples arround the estimated bounding box with 2 different scales each
    % the parameter margin specifices the distance scale between the original bounding box and the sampled negative boxes
    count = 1;
    p_ori = p;
    len = (p(3) + p(4)) / 2;
    for context_scale = [0.5, 1]
        for i = -1 : 1
            for j = -1 : 1
                if i == 0 && j == 0
                    continue;
                end
                tx = i * (p_ori(3) / 2 + len * (context_scale + min(margin / len, 0.1)));
                ty = j * (p_ori(4) / 2 + len * (context_scale + min(margin / len, 0.1)));
                p = p_ori;
                p(1) = p(1) + tx;
                p(2) = p(2) + ty;
                p(3) = 2 * len * context_scale;
                p(4) = 2 * len * context_scale;
                crop_img = im_crop(frame, round(p), out_size);
                neg_tmpl(:, count) = crop_img(:);
                count = count + 1;

                if i == 0
                    p(1) = p(1) - len * context_scale;
                    crop_img = im_crop(frame, round(p), out_size);
                    neg_tmpl(:, count) = crop_img(:);
                    count = count + 1;

                    p(1) = p(1) + 2 * len * context_scale;
                    crop_img = im_crop(frame, round(p), out_size);
                    neg_tmpl(:, count) = crop_img(:);
                    count = count + 1;
                end

                if j == 0
                    p(2) = p(2) - len * context_scale;
                    crop_img = im_crop(frame, round(p), out_size);
                    neg_tmpl(:, count) = crop_img(:);
                    count = count + 1;

                    p(2) = p(2) + 2 * len * context_scale;
                    crop_img = im_crop(frame, round(p), out_size);
                    neg_tmpl(:, count) = crop_img(:);
                    count = count + 1;
                end
            end
        end
    end

end

author: flyers

repo: drone-tracking


text:
function [S,f,R,varS,zerosp,C,Serr]=mtspectrumsegpt(data,win,params,segave,fscorr)
% Multi-taper segmented spectrum for a univariate binned point process
%
% Usage:
%
% [S,f,R,varS,zerosp,C,Serr]=mtspectrumsegpt(data,win,params,segave,fscorr)
% Input:
% Note units have to be consistent. See chronux.m for more information.
%       data (structure array of one channel of spike times;
%            also accepts 1d vector of spike times) -- required
%       win  (duration of the segments) - required.
%       params: structure with fields tapers, pad, Fs, fpass, err
%       - optional
%           tapers : precalculated tapers from dpss or in the one of the following
%                    forms:
%                   (1) A numeric vector [TW K] where TW is the
%                       time-bandwidth product and K is the number of
%                       tapers to be used (less than or equal to
%                       2TW-1).
%                   (2) A numeric vector [W T p] where W is the
%                       bandwidth, T is the duration of the data and p
%                       is an integer such that 2TW-p tapers are used. In
%                       this form there is no default i.e. to specify
%                       the bandwidth, you have to specify T and p as
%                       well. Note that the units of W and T have to be
%                       consistent: if W is in Hz, T must be in seconds
%                       and vice versa. Note that these units must also
%                       be consistent with the units of params.Fs: W can
%                       be in Hz if and only if params.Fs is in Hz.
%                       The default is to use form 1 with TW=3 and K=5
%
%	        pad		    (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%                    -1 corresponds to no padding, 0 corresponds to padding
%                    to the next highest power of 2 etc.
%			      	 e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			      	 to 512 points, if pad=1, we pad to 1024 points etc.
%			      	 Defaults to 0.
%           Fs   (sampling frequency) - optional. Default 1.
%           fpass    (frequency band to be used in the calculation in the form
%                                   [fmin fmax])- optional.
%                                   Default all frequencies between 0 and Fs/2
%           err  (error calculation [1 p] - Theoretical error bars; [2 p] - Jackknife error bars
%                                   [0 p] or 0 - no error bars) - optional. Default 0.
%       segave - (0 for don't average over segments, 1 for average) - optional - default  1
%       fscorr   (finite size corrections, 0 (don't use finite size corrections) or
%                1 (use finite size corrections) - optional
%                (available only for spikes). Defaults 0.
% Output:
%       S       (spectrum in form frequency x segments if segave=0; function of frequency if segave=1)
%       f       (frequencies)
%       R       (spike rate)
%       varS    (variance of the spectrum as a function of frequency)
%       zerosp  (0 for segments in which spikes were found, 1 for segments
%       C       (covariance matrix of the log spectrum - frequency x
%       frequency matrix)
%       Serr    (error bars) - only if err(1)>=1

if nargin < 2; error('Need data and segment information'); end;
if nargin < 3; params=[]; end;
if nargin < 4 || isempty(segave); segave=1; end;
[tapers,pad,Fs,fpass,err,trialave,params]=getparams(params);
clear params trialave
if nargin < 5 || isempty(fscorr); fscorr=0;end;

if nargout > 4 && err(1)==0; error('cannot compute error bars with err(1)=0; change params and run again'); end;

dtmp=change_row_to_column(data);
T=max(dtmp); % total length of data
minT=min(dtmp);
E=minT:win:T-win; % fictitious event triggers
win=[0 win]; % use window length to define left and right limits of windows around triggers
dtmp=createdatamatpt(dtmp,E,win); % create segmented data set
[mintime,maxtime]=minmaxsptimes(dtmp);
dt=1/Fs; % sampling time
t=mintime-dt:dt:maxtime+dt; % time grid for prolates
N=length(t); % number of points in grid for dpss
nfft=max(2^(nextpow2(N)+pad),N); % number of points in fft of prolates
[f,findx]=getfgrid(Fs,nfft,fpass); % get frequency grid for evaluation
tapers=dpsschk(tapers,N,Fs); % check tapers
[J,Msp,Nsp]=mtfftpt(dtmp,tapers,nfft,t,f,findx);% mt fft for point process times
R=Msp*Fs;
S=squeeze(mean(conj(J).*J,2)); % spectra of non-overlapping segments (averaged over tapers)
if segave==1; SS=squeeze(mean(S,2));R=mean(R);else;SS=S;end;% mean of the spectrum averaged across segments
if nargout > 3
    lS=log(SS); % log spectrum for nonoverlapping segments
%     varS=var(lS,1,2); % variance of log spectrum
    varS=var(lS',1)';% variance of the log spectrum R13
    if nargout > 4
       zerosp=zeros(1,size(data,2));
       zerosp(Nsp==0)=1;
       if nargout > 5
          C=cov(lS'); % covariance matrix of the log spectrum
          if nargout==7;
             if fscorr==1;
                Serr=specerr(SS,J,err,segave,Nsp);
             else
                Serr=specerr(SS,J,err,segave);
             end;
          end;
       end;
    end;
end;
S=SS;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% Purpose: Matlab script to compute the advection terms for incompressible Navier-Stokes

% evaluate flux vectors
fxUx = Ux.^2;  fyUx = Ux.*Uy;  fxUy = Ux.*Uy; fyUy = Uy.^2;

% save old nonlinear terms
NUxold = NUx; NUyold = NUy;

% evaluate inner-product of test function gradient and flux functions
NUx = Div2D(fxUx, fyUx); NUy = Div2D(fxUy, fyUy);

% interpolate velocity to face nodes on element faces
UxM = zeros(Nfp*Nfaces, K); UxP = zeros(Nfp*Nfaces, K);
UyM = zeros(Nfp*Nfaces, K); UyP = zeros(Nfp*Nfaces, K);
UxM(:) = Ux(vmapM); UyM(:) = Uy(vmapM); UxP(:) = Ux(vmapP); UyP(:) = Uy(vmapP);

% set '+' trace of velocity at boundary face nodes
UxP(mapI) = bcUx(mapI); UxP(mapW) = bcUx(mapW); UxP(mapC) = bcUx(mapC);
UyP(mapI) = bcUy(mapI); UyP(mapW) = bcUy(mapW); UyP(mapC) = bcUy(mapC);

% evaluate flux vectors at '-' and '+' traces at face nodes
fxUxM = UxM.^2; fyUxM = UxM.*UyM; fxUyM = UxM.*UyM; fyUyM = UyM.^2;
fxUxP = UxP.^2; fyUxP = UxP.*UyP; fxUyP = UxP.*UyP; fyUyP = UyP.^2;

% evaluate dot product of normal and velocity at face nodes
UDotNM = UxM.*nx + UyM.*ny; UDotNP = UxP.*nx + UyP.*ny;
maxvel = max(abs(UDotNM), abs(UDotNP)) ;

% evaluate maximum normal velocity at face face nodes
maxvel = reshape(maxvel, Nfp, Nfaces*K);
maxvel = ones(Nfp, 1)*max(maxvel, [], 1);
maxvel = reshape(maxvel, Nfp*Nfaces, K);

% form local Lax-Friedrichs/Rusonov fluxes
fluxUx = 0.5*( -nx.*(fxUxM-fxUxP) - ny.*(fyUxM-fyUxP)  - maxvel.*(UxP - UxM));
fluxUy = 0.5*( -nx.*(fxUyM-fxUyP) - ny.*(fyUyM-fyUyP)  - maxvel.*(UyP - UyM));

% put volume and surface terms together
NUx = NUx + LIFT*(Fscale.*fluxUx); NUy = NUy + LIFT*(Fscale.*fluxUy);

% compute (U~,V~)
UxT = ((a0*Ux+a1*Uxold) - dt*(b0*NUx + b1*NUxold))/g0;
UyT = ((a0*Uy+a1*Uyold) - dt*(b0*NUy + b1*NUyold))/g0;

repo: Aero-matlab


text:
function x = nco_abscissas ( n )

%*****************************************************************************80
%
%% NCO_ABSCISSAS computes the Newton Cotes Open abscissas.
%
%  Discussion:
%
%    The interval is [ -1, 1 ].
%
%    The abscissas are the equally spaced points between -1 and 1,
%    not including the endpoints.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the rule.
%
%    Output, real X(N), the abscissas.
%
  for i = 1 : n
    x(i) = ( ( n - i + 1 ) * ( -1.0 )   ...
           + (     i     ) * ( +1.0 ) ) ...
           / ( n     + 1 );
  end

  return
end


text:
function test_ft_artifact_threshold

% WALLTIME 00:10:00
% MEM 1gb
% DEPENDENCY

% this script serves to test the asymetric onset/offset and the peak-detection

%%

data = [];
data.label = {'ramp'};
data.fsample = 1000;
data.time{1} = (1:1000)/1000;
data.trial{1} = zeros(1,1000);
data.sampleinfo = [1 1000];

% ramp up and down
data.trial{1}(  1: 500) = 1:500;
data.trial{1}(501:1000) = 500:-1:1;

%%

cfg = [];
cfg.artfctdef.threshold.bpfilter = 'no';
cfg.artfctdef.threshold.max = 400;
[cfg1, artifact1] = ft_artifact_threshold(cfg, data);

cfg = [];
cfg.trl = artifact1;
data1 = ft_redefinetrial(cfg, data);

cfg = [];
timelock1 = ft_timelockanalysis(cfg, data1);

figure
plot(timelock1.time, timelock1.avg, '.-');

% the peak should be at t=0
[val, indx] = max(timelock1.avg);
assert(isalmostequal(timelock1.time(indx), 0));

%%

cfg = [];
cfg.artfctdef.threshold.bpfilter = 'no';
cfg.artfctdef.threshold.onset  = 400;
cfg.artfctdef.threshold.offset = 400;
[cfg2, artifact2] = ft_artifact_threshold(cfg, data);

assert(isequal(artifact2, artifact1));

%%

cfg = [];
cfg.artfctdef.threshold.bpfilter = 'no';
cfg.artfctdef.threshold.onset  = 400;
cfg.artfctdef.threshold.offset = 300;
[cfg3, artifact3] = ft_artifact_threshold(cfg, data);

assert(~isequal(artifact3, artifact1));

cfg = [];
cfg.trl = artifact3;
data3 = ft_redefinetrial(cfg, data);

cfg = [];
timelock3 = ft_timelockanalysis(cfg, data3);

figure
plot(timelock3.time, timelock3.avg, '.-'); % it should be asymetric

% the peak should be at t=0
[val, indx] = max(timelock3.avg);
assert(isalmostequal(timelock3.time(indx), 0, 'abstol', 100*eps));

%%

cfg = [];
cfg.artfctdef.threshold.bpfilter = 'no';
cfg.artfctdef.threshold.onset  = 1;
cfg.artfctdef.threshold.offset = 1;
[cfg4, artifact4] = ft_artifact_threshold(cfg, data);

assert(isequal(artifact4, [1 1000 -499]));

%%
% invert the data and invert the thresholds

data.trial{1} = -data.trial{1};

cfg = [];
cfg.artfctdef.threshold.bpfilter = 'no';
cfg.artfctdef.threshold.onset  = -400;
cfg.artfctdef.threshold.offset = -300;
[cfg5, artifact5] = ft_artifact_threshold(cfg, data);

assert(isequal(artifact5, artifact3));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [tstop,res] = testProblem(optObj,solvers)
%TESTPROBLEM  Test an OPTI Problem against all available solvers
%
%   [times,results] = testProblem(optObj) runs the supplied OPTI problem
%   against all available solvers for the given problem type.
%
%   [times,results] = testProblems(optObj,solvers) runs only the selected
%   solvers across the supplied problem.

%   Copyright (C) 2011 Jonathan Currie (IPL)

%Setup default arguments
if(~exist('solvers','var') || isempty(solvers))
    solvers = optiSolver(['all_' optObj.prob.type]);
end
if(~iscell(solvers))
    solvers = {solvers};
end

%Remove troublesome solvers
ind = strcmpi(solvers,'lipsol');
if(any(ind))
    solvers(ind) = [];
end

%Setup problem generation function (for warm up)
switch(lower(optObj.prob.type))
    case 'lp'
        testSet = @lp_prob;
    case 'milp'
        testSet = @milp_prob;
    case 'bilp'
        testSet = @bilp_prob;
    case 'qp'
        testSet = @qp_prob;
    case 'qcqp'
        testSet = @qcqp_prob;
    case 'sdp'
        testSet = @sdp_prob;
    case 'miqp'
        testSet = @miqp_prob;
    case 'nlp'
        testSet = @nlp_prob;
        if(isempty(optObj.prob.x0))
            error('Please include x0 in optiprob when using this function on NLPs');
        end
    case 'minlp'
        testSet = @minlp_prob;
        if(isempty(optObj.prob.x0))
            error('Please include x0 in optiprob when using this function on MINLPs');
        end
    otherwise
        error('Unknown problem type: %s',optObj.prob.type);
end

%Remove troublesome solvers
if(strcmpi(optObj.prob.type,'lp'))
    ind = strcmpi(solvers,'SEDUMI'); solvers(ind) = []; %no good with neq = nvar
end

%Setup no runs
len = length(solvers);
%confirm solver exists!
for i = 1:len
    rub = optiSolver(solvers{i}); %#ok<NASGU>
end

%Setup wait bar
h = waitbar(0,['Solver: ' upper(solvers{1})],'name','OPTI Problem Test','CreateCancelBtn','setappdata(gcbf,''canceling'',1)');

%Preallocate
tstop = zeros(len,1);
res = struct('sol',{cell(len,1)},'fval',zeros(len,1),'ok',zeros(len,1),'msg',{cell(len,1)});

%Clear display warnings etc
opts = optiset(optObj.prob,'display','off','warnings','off');

try
    for i = 1:len
        %Check Cancel
        if getappdata(h,'canceling')
            break;
        end
        %Assign the new solver
        solver = solvers{i};
        opts = optiset(opts,'solver',solver);
        %Initialize with two solves of problem 1
        prob = testSet(1);
        Opt = opti(prob,opts);
        solve(Opt); solve(Opt);
        %Now run our problem
        Opt = opti(getProb(optObj),opts);
        [x,sval,~,stat] = solve(Opt);
        res.sol{i} = x;
        res.fval(i) = sval;
        res.ok(i) = checkSol(Opt);
        res.msg{i} = stat.Status;
        tstop(i) = stat.Time;
        %Update Waitbar
        waitbar(i/len,h,['Solver: ' upper(solvers{min(len,i+1)})]);
    end
catch ME
    delete(h);
    rethrow(ME);
end

delete(h);

%Draw Table
fprintf('\nTest Problem [%s] Result\n',Opt.prob.type);
fprintf('--------------------------------------------------\n');
fprintf('  Solver    Status        Fval         Time\n');
for i = 1:len
    fprintf('%8s    %4s    %+12.4f     %1.4fs\n',solvers{i},status(res.ok(i)),res.fval(i),tstop(i));
end
fprintf('--------------------------------------------------\n');

function msg = status(ok)

if(ok)
    msg = 'OK';
else
    msg = 'FAIL';
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ctmediumphantom,ctdensityphantom] = dicomrt_ctcreate(inpct_matrix,ct_xmesh,ct_ymesh,ct_zmesh,rtdose_xmesh,rtdose_ymesh,rtdose_zmesh,plotramp,int)
% dicomrt_ctcreate(inpct_matrix,ct_xmesh,ct_ymesh,ct_zmesh,rtdose_xmesh,rtdose_ymesh,rtdose_zmesh,plotramp,int)
%
% Create a ct phantom for BEAM/DOSXYZ simulation. Light Version.
%
% inpct_matrix is the CT dataset
% ct_xmesh, ct_ymesh, ct_zmesh are the coordinates of the voxels for the CT dataset
% rtdose_xmesh, rtdose_ymesh, rtdose_zmesh are the coordinates of the voxels
% of the 3D dose distribution
% plotramp is an option for plotting the CT ramp (~=0 plot, =0[default] no plot)
% Int is a parameter for interactive session.
%     If int ~=1 or if int is not given session is not interactive
%     and default parameters (e.g. directories names) will be used. If int = 1 the user will
%     be asked to set some parameters.
%
% Parameters for this specific m file are:
%
% nmat				    number of material
% materials			    materials name
% medium				materials number
% mat_ct_up_bound		materials upper bound (see dosxyz user's manual)
% density_lo_bound		density lower bound (see dosxyz user's manual)
% density_up_bound		density upper bound (see dosxyz user's manual)
% estepe				estepe (see dosxyz user's manual)
% filename			    ct phantom filename for export
%
% See also dicomrt_loaddose, dicomrt_rotate180x
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check number of argument
error(nargchk(8,9,nargin))

if nargin<=8
    int=0;
end

% Check case and set-up some parameters and variables
[ct_matrix_temp,type_dose,label,PatientPosition]=dicomrt_checkinput(inpct_matrix);
ct_matrix=dicomrt_varfilter(ct_matrix_temp);

[ct_xmesh_2d,ct_ymesh_2d]=dicomrt_build2dgrid(ct_xmesh,ct_ymesh);
[dose_xmesh_2d,dose_ymesh_2d]=dicomrt_build2dgrid(rtdose_xmesh,rtdose_ymesh);

% 1a) resize ct_matrix using dose_xmesh, dose_ymesh, dose_zmesh
% loop over dose_zmesh size since the rt plan could have been performed on a subset of the original ct images
for i=1:length(rtdose_zmesh)
    if i==1 % locate z position of the first slice to start with
      first_zlocation=rtdose_zmesh(i);
      for k=1:length(ct_zmesh)
         temp=num2str(ct_zmesh(k));
         array=char(temp,num2str(first_zlocation));
         if array(1,:)==array(2,:);
            break
         end
      end
      zstart=k; % match with dose matrix
   end
   ctphantom(:,:,i)=interp2(ct_xmesh_2d,ct_ymesh_2d,ct_matrix(:,:,k+i-1), ...
       dose_xmesh_2d,dose_ymesh_2d,'nearest');
end

% 1b) Cut off negative values form outside CT FOV
ctphantom(find((ctphantom)<=0))=10; % a CT value for air

% 2a) Retrieve ctphantom parameters: if int==1 start interactive session
if isnumeric(int)==1 &  int==1
    nmat = input('Input the number of materials [Enter for default]: '); % Number of material
    if isempty(nmat)==1
        nmat=4;
        %disp('The following materials will be used: ');
        materials          = char('AIR521ICRU','LUNG521ICRU','ICRUTISSUE521ICRU','ICRPBONE521ICRU');
        medium             = [1;     2;     3;     4];
        mat_ct_up_bound    = [50;    300;   1125;  5000];
        density_lo_bound   = [0.001; 0.044; 0.302; 1.101];
        density_up_bound   = [0.044; 0.302; 1.101; 3.48];
        estepe             = [0.25;  0.25;  0.25;  0.25];
    else
        medium = [1:nmat]';
        for i=1:nmat
            mat = input('Input the name of material: ','s'); % Name of material
            if i==1
                materials = mat;
            else
                materials = char(materials,mat);
            end
            mat_ct_up_bound(i) = input('Input the material ct upper bound: ');
            density_lo_bound(i) = input('Input the density ct lower bound: ');
            density_up_bound(i) = input('Input the density ct upper bound: ');
            estepe(i) = input('Input estepe for this material: ');
        end
        estepe=estepe';
    end % ctphantom parameters retrieved
else
    nmat=4;
    %disp('The following materials will be used: ');
    materials          = char('AIR521ICRU','LUNG521ICRU','ICRUTISSUE521ICRU','ICRPBONE521ICRU');
    medium             = [1;     2;     3;     4];
    mat_ct_up_bound    = [50;    300;   1125;  5000];
    density_lo_bound   = [0.001; 0.044; 0.302; 1.101];
    density_up_bound   = [0.044; 0.302; 1.101; 3.84];
    estepe             = [0.25;  0.25;  0.25;  0.25];
end

% 2b) Plot CT ramp
if plotramp~=0
    disp('Plotting CT conversion ramp being used: ');
    dicomrt_plotctramp(inputname(1),materials,mat_ct_up_bound,density_lo_bound,density_up_bound);
end

% 3) Build ctphantom

% 3a) CT to medium
ctmediumphantom=ctphantom;
for i=1:nmat
    if i==1
        ctmediumphantom(find((ctmediumphantom<=mat_ct_up_bound(i) & ...
            ctmediumphantom>0)))=medium(i);
    else
        ctmediumphantom(find((ctmediumphantom<=mat_ct_up_bound(i) & ...
            ctmediumphantom>mat_ct_up_bound(i-1))))=medium(i);
    end
end

% 3b) Convert CT data to density
ctdensityphantom=ones(size(ctmediumphantom,1),size(ctmediumphantom,2),size(ctmediumphantom,3));
for i=1:nmat
    if i==1
        a=density_lo_bound(i);
        b=density_up_bound(i)-density_lo_bound(i);
        c=mat_ct_up_bound(i);
        d=ctphantom(find((ctmediumphantom==medium(i))));
        ctdensityphantom(find((ctmediumphantom==medium(i)))) = a+(b/c)*d;
    else
        a=density_lo_bound(i);
        b=density_up_bound(i)-density_lo_bound(i);
        c=mat_ct_up_bound(i)-mat_ct_up_bound(i-1);
        d=ctphantom(find((ctmediumphantom==medium(i))))-mat_ct_up_bound(i-1);
        ctdensityphantom(find((ctmediumphantom==medium(i)))) = a+(b/c)*d;
    end
end

% 5) Export to file.
%
% ImagePositionPatient is the (x,y,z) coordinate of the first pixel (mm).
% Therefore the mesh produced so far represent the coordinates of the center of the voxels.
% CTphantom need the coordinates of the boundaries of the voxels.
%
if isnumeric(int)==1 &  int==1
    export = input('Do you want to export the ctphantom ? Y/N [N]: ','s');
    if export == 'Y' | export == 'y';
        % Prepare to export
        % pixel spacing
        xthick=abs(rtdose_xmesh(1)-rtdose_xmesh(2));
        ythick=abs(rtdose_ymesh(1)-rtdose_ymesh(2));
        zthick=abs(rtdose_zmesh(1)-rtdose_zmesh(2));

        x_bound=dicomrt_createboundgrid(rtdose_xmesh);
        y_bound=dicomrt_createboundgrid(rtdose_ymesh);
        z_bound=dicomrt_createboundgrid(rtdose_zmesh);

        filename = input('Input the egs4phant filename (no ext): ','s'); % Name of the file where to store ctphantom

        if filename == 'N' | filename == 'n';
            warning('dicomrt_ctcreate: No filename was input. Filename "filename" will be used');
            filename=['filename.egs4phant'];
            dicomrt_writeegs4phant(PatientPosition,nmat,materials,estepe,length(rtdose_xmesh),length(rtdose_ymesh), ...
                length(rtdose_zmesh),x_bound,y_bound,z_bound,ctmediumphantom,ctdensityphantom,filename);
        else
            filename=[filename,'.egs4phant'];
            dicomrt_writeegs4phant(PatientPosition,nmat,materials,estepe,length(rtdose_xmesh),length(rtdose_ymesh), ...
                length(rtdose_zmesh),x_bound,y_bound,z_bound,ctmediumphantom,ctdensityphantom,filename);
        end
    else
        disp('No export will be performed. Data will be stored on disk');
        save
    end
else
    % Prepare to export
    % pixel spacing
    xthick=abs(rtdose_xmesh(1)-rtdose_xmesh(2));
    ythick=abs(rtdose_ymesh(1)-rtdose_ymesh(2));
    zthick=abs(rtdose_zmesh(1)-rtdose_zmesh(2));

    x_bound=dicomrt_createboundgrid(rtdose_xmesh);
    y_bound=dicomrt_createboundgrid(rtdose_ymesh);
    z_bound=dicomrt_createboundgrid(rtdose_zmesh);

    filename=inputname(1);
    filename=[filename,'.egs4phant'];

    dicomrt_writeegs4phant(PatientPosition,nmat,materials,estepe,length(rtdose_xmesh),length(rtdose_ymesh), ...
        length(rtdose_zmesh),x_bound,y_bound,z_bound,ctmediumphantom,ctdensityphantom,filename);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% Test file for ADCHEBFUN exp, log related functions

function pass = test_expLog

% List of trigonometric functions to test.
funcList = {@exp, @expm1, @log, @log10, @log1p, @log2};

% Call the ADCHEBFUN testUnary() method to do the tests.
pass = adchebfun.testUnary(funcList);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function  varargout = partitionData(Ndata, varargin)
% PARTITIONDATA Partition a vector of indices into random sets
% [a,b,c,...] = partitionData(N, 0.3, 0.2, 0.5, ...)
%
% Examples:
% [a,b,c]=partitionData(105,0.3,0.2,0.5);
% a= 1:30, b=32:52, c=52:105 (last bin gets all the left over)

Npartitions = length(varargin);
perm = randperm(Ndata);
%perm = 1:Ndata;
ndx = 1;
for i=1:Npartitions
  pc(i) = varargin{i};
  Nbin(i) = fix(Ndata*pc(i));
  low(i) = ndx;
  if i==Npartitions
    high(i) = Ndata;
  else
    high(i) = low(i)+Nbin(i)-1;
  end
  varargout{i} = perm(low(i):high(i));
  ndx = ndx+Nbin(i);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function v=scatterquad2(x,y,z)
% SCATTERQUAD2 - calculates the volume under a surface defined by scattered points
%
% Usage:
%   v = scatterquad2(x,y,z)
%
% Inputs:
%   X,Y,Z are vectors of equal size specifying the coordinates of the points
%   defining the surface.  Z can also be a scalar or a function handle accepting
%   two inputs.
%
% Outputs:
%   V is the volume under the surface defined by linear interpolation of Z on
%   the Delaunay triangulation of the points (X(i),Y(i)), assuming that z=0
%   outside the convex hull of the points (X(i),Y(i)).
%
% Method:
%   The linear interpolation is a linear combination of basis functions of the
%   form z=t if (x,y) is in triangle T such that
%     (x,y)=r*(u1,v1)+s*(u2,v2)+t*(u3,v3)  where r+s+t=1 (all non-negative) and
%   (u,v) are the vertices of T, or z=0 otherwise.  The integral of z is
%   (1/3)*Area(T) where the area of T is determined by the cross product of two
%   edge vectors.
%
% Examples:
%   load seamount
%   scatterquad2(x,y,z-min(z)) % returns 190.7996
%   inR = (x>=211.1 & x<=211.4 & y>=-48.35 & y<=-48);
%   scatterquad2(x,y,(z-min(z)).*inR) % returns 142.3083
%   scatterquad2(x,y,1) % returns 0.2696
%
% For functions that can be evaluated at arbitrary points, use DBLQUAD or
% QUAD2D.  For regular grids, use TRAPZ twice.

if nargin<3
  error('scatterquad2:nargin','Expect 3 input arguments');
elseif ~isnumeric(x) || ~isnumeric(y)
  error('scatterquad2:xytype','Inputs X and Y must be numeric');
elseif ~isequal(numel(x),numel(y))
  error('scatterquad2:xysize','Inputs X and Y must have same number of elements');
elseif numel(x)<3
  error('scatterquad2:xsize','Inputs X and Y must have at least 3 elements');
elseif ~isnumeric(z)
  if isa(z,'function_handle')
    try
      z=z(x,y);
      if ~isnumeric(z)
        error('scatterquad2:zftype','Input Z failed to return a numeric array');
      end;
    catch ME
      error('scatterquad2:zfeval','Input Z failed to evaluate at (X,Y)');
    end;
  else
    error('scatterquad2:ztype','Input Z must be numeric or a function handle');
  end;
end;

if ~isequal(numel(z),numel(x))
  if numel(z)==1
    z=repmat(z,numel(x),1);
  else
    error('scatterquad2:zsize','Input Z must be scalar or the same size as X and Y');
  end;
end;

DT=DelaunayTri(x,y);
i=DT(:,1);
j=DT(:,2);
k=DT(:,3);
A=abs((x(j)-x(i)).*(y(k)-y(i))-(y(j)-y(i)).*(x(k)-x(i)))/2; % area of triangles
v=sum((z(i)+z(j)+z(k)).*A)/3;

end % main function scatterquad2(...)


text:
function [count,cycles] = CountSpikesPerCycle(spikes,phases)

%CountSpikesPerCycle - Count number of spikes per LFP cycle.
%
% Count the number of spikes per cycle in the ongoing oscillatory LFP
% (e.g. during theta).
%
%  USAGE
%
%    [count,cycles] = CountSpikesPerCycle(spikes,phases)
%
%    spikes         list of spike timestamps
%    phases         instantaneous phases in radians (see Phase)
%
%  OUTPUT
%
%    count          spike count in each cycle
%    cycles         list of [start,stop] times for each cycle
%
%  SEE
%
%    See also Phase.

% Copyright (C) 2004-2011 by Michaël Zugaro
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% Check number of parameters
if nargin < 2,
	error('Incorrect number of parameters (type ''help CountSpikesPerCycle'' for details).');
end

% Check parameter sizes
if ~isdvector(spikes),
	error('Parameter ''spikes'' is not a vector (type ''help CountSpikesPerCycle'' for details).');
end
if ~isdmatrix(phases) | size(phases,2) ~= 2,
	error('Parameter ''phases'' is not an Nx2 matrix (type ''help CountSpikesPerCycle'' for details).');
end
isradians(phases(:,2));

% Find theta peaks
p = phases(:,2);
p(p>pi) = p(p>pi) - 2*pi;
[up,unused] = ZeroCrossings([phases(:,1) p]);
start = find(up);
cycles = [phases(start(1:end-1),1) phases(start(2:end),1)];

% Intervals between successive theta peaks
[unused,interval] = InIntervals(spikes,cycles);

% Count
used = interval ~= 0;
count = Accumulate(interval(used));
%  count = count(interval(used));


text:
function [val,retData]=unbiasedMomentCumulant(x,r,isMoment,algorithm,retData)
%%UNBIASEDMOMENTCUMULANT Given a set of samples of a scalar distribution,
%           return an unbiased estimate of the rth order central moment (h
%           statistic), which is E{(x-mu)^r}, where E is the expected value
%           operation, when given only samples of x and not knowing the
%           mean mu. Alternatively, one could compute an unbiased estimate
%           of the rth order cumulant (k statistic). Note that for r=0, if
%           a moment is desired, then the mean is returned, rather than
%           returning 0. If r is a vector, one can estimate generalized h
%           statistics, which are an estimate of
%           E{(x-mu)^r(1)}*E{(x-mu)^r(2)}*...*E{(x-mu)^r(v)} However,
%           algorithm 2 is used for generalized h statistics and its
%           complexity increases rapidly with the number of samples.
%
%INPUTS: x A 1XN or NX1 array of scalar samples of the distribution.
%        r The scalar order of the cumulant or of the central moment.
%          Alternatively, if this is a vector, then one will get a
%          generalized h statistic, and in that instance, algorithm 2 is
%          required.
% isMoment A parameter selecting whether a moment or a cumulant is desired.
%          The default if omitted or an empty matrix is passed is true,
%          indicating that a central moment is desired.
% algorithm An optional parameter selecting which algorithm to use.
%          possible values are:
%          0 Use an explicit solution with precomputed coefficients. These
%            are only available for scalar r=1,...,5. The explicit
%            solutions use the parameters in Table I of [1]. This is the
%            default if r<=5 is scalar.
%          1 Use Dwyer's method in [1] to determine the coefficients of the
%            expression of the moment in terms of the sample power sums
%            sum(x^k) for k=1 to r.
%          2 Use the symmetric function based algorithm of [2]. This
%            algorithm can only estimate central moments, not cumulants.
%            The complexity of this algorithm increases rapidly with N.
%            This can be used with a vector r and is thus the default when
%            r is a vector.
% retData If algorithm=1 is used, then the retData value is returned. If
%         one calls this function passing back retData and keeps all other
%         inputs the same except replacing x with a different set of
%         samples, then the function will execute faster. One should not
%         pass back retData if any inputs other than x have changed.
%
%OUTPUTS: val The scalar value of the moment.
%     retData If algorithm=1, then this is a structure that can be passed
%     back so that if the function is called a second time with the same
%     inputs, changing only x to different samples, then the algorithm will
%     be faster.
%
%When mu is not available and must be found from the samples, as is the
%case here, then the average mean((x-mean(x)).^p) to approximate
%E{(x-mu)^p} is a biased estimator. Hence, h statistics are more
%complicated.
%
%EXAMPLE 1:
%Here, we estimate the sixth central moment of the scalar Laplace
%distribution from 1000 samples. Since the distribution is symmetric and
%unbounded, we can assume that if the estimator is unbiased, the averaging
%a large number of these estimates will approach the true moment value.
%That is done here and compared to the true moment value. One can also see
%how retData is used to speed up the algorithm.
% n=1000;
% lambda=0.1;
% mu=2;
% Gamma=1;
% momentNum=6;
% numRuns=1e4;
%
% hAvg=0;
% retData=[];
% for curRun=1:numRuns
%     x=LaplaceD.rand(n,lambda,mu,Gamma);
%     [val,retData]=unbiasedMomentCumulant(x,momentNum,true,[],retData);
%     hAvg=hAvg+val;
% end
% avgMomentEst=hAvg/numRuns
% muList=zeros(momentNum,1);
% for k=1:momentNum
%     muList(k)=LaplaceD.momentGenFun(lambda,mu,Gamma,k);
% end
% centralMoment=raw2CentralMoments(muList);
% trueMoment=centralMoment(end)
%One should see that the average moment and the true moment should agree to
%a few decimal places.
%
%EXAMPLE 2:
%Now, we estimate the product of the second central moment and the fourth
%central moment. This necessitates the use of algorith 2, which can be
%rather slow. This example will typically take a few seconds to run.
% rng(0);%To make the results reproducible
% n=10;
% lambda=0.25;
% mu=2;
% Gamma=1;
% p=[2,4];
% numRuns=1000;
%
% hAvg=0;
% for curRun=1:numRuns
%     x=LaplaceD.rand(n,lambda,mu,Gamma);
%     hAvg=hAvg+unbiasedMomentCumulant(x,p);
% end
% avgMomentEst=hAvg/numRuns
% muList=zeros(max(p),1);
% for k=1:max(p)
%     muList(k)=LaplaceD.momentGenFun(lambda,mu,Gamma,k);
% end
% centralMoment=raw2CentralMoments(muList);
% trueMoment=prod(centralMoment(p))
%The average moment estimate will be about 0.0998 and the true moment will
%be about 0.0938.
%
%REFERENCES:
%[1] P. S. Dwyer, "Moments of any rational integral isobaric sample moment
%    function," The Annals of Mathematical Statistics, vol. 8, no. 1, pp.
%    21-65, Mar. 1937.
%[2] D. S. Tracy and B. C. Gupta, "Generalized h-statistics and other
%    symmetric functions," The Annals of Statistics, vol. 2, no. 4, pp.
%    837-844, 1974.
%
%November 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release

if(nargin<3||isempty(isMoment))
    isMoment=true;
end

if(nargin<4||isempty(algorithm))
    if(isscalar(r))
        if(r<=5)
            algorithm=0;%Use the explicit solution.
        else
            algorithm=1;%Use the power sum solution of Dwyer.
        end
    else
        algorithm=2;

        if(isMoment==false)
            error('Vector values of r are not allowed for algorithm~=2.')
        end
    end
end

if(nargin<5)
    retData=[];
end

if(isMoment)%h-statistics
    switch(algorithm)
        case 0
            val=hStatExplicit(x,r);
            retData=[];
        case 1
            [val,retData]=momCumDwyer(x,r,false,retData);
        case 2
            val=genHStat(x,r);

            retData=[];
        otherwise
            error('Unknown Algorithm specified.')
    end
else%k-statistics
    switch(algorithm)
        case 0
            val=kStatExplicit(x,r);
            retData=[];
        case 1
            [val,retData]=momCumDwyer(x,r,true,retData);
        otherwise
            error('Unknown Algorithm specified.')
    end
end
end

function [momentVal,retData]=momCumDwyer(x,r1,isCumulant,retData)
%%MOMCUMDWYER Compute unbiased central moments (h statistics) or cumulants
%             (k statistics) using Dwyer's algorithm, described in [1]. The
%             definition of the symmetric powers in the algorithm can be
%             bettwer understood when considering [2].
%
%The algorithm uses a weird definition of a product, which involves
%appending the indices of various "b" variables on each other. Thus, an
%initial phase of the algorithm involves just determinign the indices of
%the necessary coefficients. Sections 17 and 18 of [1] give the b variables
%for k statistics and h statistics. Section 16 discusses how to turn the b
%variables into a variables. The expression for turning the a variables
%into moments, F_r, is in Section 1.
%
%REFERENCES:
%[1] P. S. Dwyer, "Moments of any rational integral isobaric sample moment
%    function," The Annals of Mathematical Statistics, vol. 8, no. 1, pp.
%    21-65, Mar. 1937.
%[2] P. S. Dwyer, "Combined expansions of products of symmetric power sums
%    and of sums of symmetric power products with application to sampling,"
%    The Annals of Mathematical Statistics, vol. 9, no. 1, pp. 1-47, Mar.
%    1938.
%
%November 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

n=length(x);

numParts=numberOfPartitions(r1,true);

if(nargin<4||isempty(retData))
    %If we have to compute the coefficients.
    %AParts holds the indices of the b that go into the sum to form the a as
    %well as their coefficients.
    aParts=cell(r1,2);
    for r=1:r1
        %Allocate space for the indices of the b values needed and their
        %coefficients.
        bIdx=cell(numParts(r+1),1);
        bCoeff=zeros(numParts(r+1),1);
        [thePartition,theData]=getNextPartition(r);
        curB=1;
        while(~isempty(thePartition))
            s=theData.d;
            p=theData.r(1:s);
            piVals=theData.m(1:s);
            rho=sum(piVals);

            [p,idx]=sort(p,'descend');
            piVals=piVals(idx);

            bCoeff(curB)=partComboNum(r,p,piVals)*(-1)^(rho-1)*factorial(rho-1);
            idx=zeros(rho,1);
            startIdx=1;
            for k=1:s
                numRep=piVals(k);
                sel=startIdx:(startIdx+numRep-1);
                idx(sel)=p(k);
                startIdx=startIdx+numRep;
            end
            bIdx{curB}=idx;

            [thePartition,theData]=getNextPartition(r,theData);
            curB=curB+1;
        end
        aParts{r,1}=bIdx;
        aParts{r,2}=bCoeff;
    end

    %Now, determine the b indices needed to construct the sum in terms of the
    %as. (The Fr equation).
    aPartNeeded=cell(r1,3);
    %This stores the a indices, the coefficient of the a term and the powers of
    %the P sums for each a term.
    [thePartition,theData]=getNextPartition(r1);
    curA=1;
    while(~isempty(thePartition))
        s=theData.d;
        p=theData.r(1:s);
        piVals=theData.m(1:s);
        rho=sum(piVals);

        [p,idx]=sort(p,'descend');
        piVals=piVals(idx);
        coeffVal=partComboNum(r1,p,piVals);

        idx=zeros(rho,1);
        startIdx=1;
        for k=1:s
            numRep=piVals(k);
            sel=startIdx:(startIdx+numRep-1);
            idx(sel)=p(k);
            startIdx=startIdx+numRep;
        end

        aPartNeeded{curA,1}=idx;
        aPartNeeded{curA,2}=coeffVal;
        aPartNeeded{curA,3}=piVals(:);
        aPartNeeded{curA,4}=p(:);

        [thePartition,theData]=getNextPartition(r1,theData);
        curA=curA+1;
    end

    %To build each of the required "a" terms with the parts needed from the
    %single index a terms.
    numCoeffs=curA-1;
    %Allocate space for the a coefficients. These must be built from the suffix
    %multiplication rule of Section 16.
    aCoeffs=zeros(numCoeffs,1);
    for curPart=1:numCoeffs
        idx=aPartNeeded{curPart,1};
        r=length(idx);

        bIdxList=aParts{idx(1),1};
        bCoeffList=aParts{idx(1),2};
        for k=2:r
            bIdxCur=aParts{idx(k),1};
            bCoeffCur=aParts{idx(k),2};
            numB=length(bIdxList);
            numCur=length(bIdxCur);

            bIdxListNew=cell(numB*numCur,1);
            bCoeffListNew=zeros(numB*numCur,1);
            newCur=1;
            for i1=1:numB
                for i2=1:numCur
                    bIdxListNew{newCur}=[bIdxList{i1};bIdxCur{i2}];
                    bCoeffListNew(newCur)=bCoeffList(i1)*bCoeffCur(i2);

                    newCur=newCur+1;
                end
            end
            bIdxList=bIdxListNew;
            bCoeffList=bCoeffListNew;
        end

        %We now have the b indices and coefficients to construct the a
        %coefficients once we compute the b values.
        numB=length(bIdxList);
        aSum=0;
        for curB=1:numB
            bIdx=bIdxList{curB};
            bCoeff=bCoeffList(curB);

            numBIdx=length(bIdx);
            s=numBIdx;

            bIdx=sort(bIdx,'descend');
            p=zeros(s,1);
            piVals=zeros(s,1);

            p(1)=bIdx(1);
            numUnique=1;
            piVals(1)=1;
            for curIdx=2:numBIdx
                if(bIdx(curIdx)==bIdx(curIdx-1))
                    s=s-1;
                    piVals(numUnique)=piVals(numUnique)+1;
                else
                    numUnique=numUnique+1;
                    p(numUnique)=bIdx(curIdx);
                    piVals(numUnique)=1;
                end
            end
            %Note that s==numUnique.
            piVals=piVals(1:s);
            p=p(1:s);
            rho=sum(piVals);

            if(isCumulant)%For Fisher k statistics
                b=(-1)^(rho-1)*factorial(rho-1)/fallingFactorial(n,rho);
            else%For h statistics
                if(s==1&&piVals(1)==1&&p(1)==r1)
                    AVal=1;
                elseif(p(1)>1&&piVals(1)==1&&s==2&&p(2)==1)
                    AVal=(-1)^piVals(2);
                elseif(p(1)==1&&s==1&&piVals(1)==r1)
                    AVal=(-1)^(r1-1)*(r1-1);
                else
                    AVal=0;
                end
                b=AVal/fallingFactorial(n,rho);
            end

            aSum=aSum+bCoeff*b;
        end

        aCoeffs(curPart)=aSum;
    end

    retData.aPartNeeded=aPartNeeded;
    retData.aCoeffs=aCoeffs;
else
    aPartNeeded=retData.aPartNeeded;
    aCoeffs=retData.aCoeffs;
    numCoeffs=length(aCoeffs);
end

%We have all of the coefficients. To get this moment with a generic length
%n set of data, we just need the power sums up to degree r.
P=zeros(r1,1);
for k=1:r1
    P(k)=sum(x.^k);
end

momentVal=0;
for curA=1:numCoeffs
    aCur=aCoeffs(curA);
    coeff=aPartNeeded{curA,2};
    piVals=aPartNeeded{curA,3};
    idx=aPartNeeded{curA,4};

    momentVal=momentVal+coeff*aCur*prod(P(idx).^piVals);
end

end

function val=partComboNum(r,p,piVals)
%%PARTCOMBONUM The expression for this multinomial value is given on page
%              23 of [1].
%
%REFERENCES:
%[1] P. S. Dwyer, "Moments of any rational integral isobaric sample moment
%    function," The Annals of Mathematical Statistics, vol. 8, no. 1, pp.
%    21-65, Mar. 1937.
%
%November 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

    %For large r, use this form to help reduce overflow.
    val=exp(gammaln(r+1)-sum(piVals.*gammaln(p+1))-sum(gammaln(piVals+1)));
end

function val=hStatExplicit(x,r)
%%HSTATEXPLICIT This implements the explicit h statistics of orders 1 to 5
%               that are given in the table I of [1] using the formula for
%               Fr given in Section 1 of [1]. h statistics are unbiased
%               estimates of the moment about the mean.
%
%INPUTS: x A 1XN or NX1 array of scalar samples of the distribution.
%        r The order of the moment. 1<=r<=5.
%
%OUTPUT: val The rth order h statistic. This is a scalar value.
%
%REFERENCES:
%[1] P. S. Dwyer, "Moments of any rational integral isobaric sample moment
%    function," The Annals of Mathematical Statistics, vol. 8, no. 1, pp.
%    21-65, Mar. 1937.
%
%November 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

n=length(x);
switch(r)
    case 1
        a1=1/n;

        P1=sum(x);
        val=a1*P1;
    case 2
        a2=1/(n-1);
        a11=-1/(n*(n-1));

        P1=sum(x);
        P2=sum(x.^2);

        val=a2*P2+a11*P1^2;
    case 3
        a3=n/((n-1)*(n-2));
        a21=-1/((n-1)*(n-2));
        c21=3;
        a111=2/(n*(n-1)*(n-2));

        P1=sum(x);
        P2=sum(x.^2);
        P3=sum(x.^3);

        val=a3*P3+c21*a21*P1*P2+a111*P1^3;
    case 4
        n13=(n-1)*(n-2)*(n-3);
        n4=n*n13;

        a4=(n^2-2*n+3)/n13;
        a31=-(n^2-2*n+3)/n4;
        c31=4;
        a22=-(2*n-3)/n4;
        c22=3;
        a211=1/n13;
        c211=6;
        a1111=-3/n4;

        P1=sum(x);
        P2=sum(x.^2);
        P3=sum(x.^3);
        P4=sum(x.^4);

        val=a4*P4+c31*a31*P3*P1+c22*a22*P2^2+c211*a211*P2*P1^2+a1111*P1^4;
    case 5
        n14=(n-1)*(n-2)*(n-3)*(n-4);
        n5=n*n14;

        a5=n*(n^2-5*n+10)/n14;
        a41=-(n^2-5*n+10)/n14;
        c41=5;
        a32=-(n-2)/n14;
        c32=10;
        a311=(n^2-4*n+8)/n5;
        c311=10;
        a221=(2*n-4)/n5;
        c221=15;
        a2111=-1/n14;
        c2111=10;
        a11111=4/n5;

        P1=sum(x);
        P2=sum(x.^2);
        P3=sum(x.^3);
        P4=sum(x.^4);
        P5=sum(x.^5);
        val=P5*a5+c41*a41*P4*P1+c32*a32*P3*P2+c311*a311*P3*P1^2+c221*a221*P2^2*P1+c2111*a2111*P2*P1^3+a11111*P1^5;
    otherwise
        error('The explicit solutions are only available up to order 5.');
end
end

function val=kStatExplicit(x,r)
%%KSTATEXPLICIT This implements the explicit Fisher k statistics of orders
%               1 to 5 that are given in the table I of [1] using the
%               formula for Fr given in Section 1 of [1]. h statistics are
%               unbiased estimates of the cumulants of the distribution. Up
%               to order 3, they are the same as the h statistics.
%
%INPUTS: x A 1XN or NX1 array of scalar samples of the distribution.
%        r The order of the moment. 1<=r<=5.
%
%OUTPUT: val The rth order h statistic. This is a scalar value.
%
%REFERENCES:
%[1] P. S. Dwyer, "Moments of any rational integral isobaric sample moment
%    function," The Annals of Mathematical Statistics, vol. 8, no. 1, pp.
%    21-65, Mar. 1937.
%
%November 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

n=length(x);
switch(r)
    case 1
        a1=1/n;

        P1=sum(x);
        val=a1*P1;
    case 2
        a2=1/(n-1);
        a11=-1/(n*(n-1));

        P1=sum(x);
        P2=sum(x.^2);

        val=a2*P2+a11*P1^2;
    case 3
        a3=n/((n-1)*(n-2));
        a21=-1/((n-1)*(n-2));
        c21=3;
        a111=2/(n*(n-1)*(n-2));

        P1=sum(x);
        P2=sum(x.^2);
        P3=sum(x.^3);

        val=a3*P3+c21*a21*P1*P2+a111*P1^3;
    case 4
        n22=(n-2)*(n-3);
        n13=(n-1)*n22;
        n4=n*n13;

        a4=n*(n+1)/n13;
        a31=-(n+1)/n13;
        c31=4;
        a22=-1/n22;
        c22=3;
        a211=2/n13;
        c211=6;
        a1111=-6/n4;

        P1=sum(x);
        P2=sum(x.^2);
        P3=sum(x.^3);
        P4=sum(x.^4);

        val=a4*P4+c31*a31*P3*P1+c22*a22*P2^2+c211*a211*P2*P1^2+a1111*P1^4;
    case 5
        n14=(n-1)*(n-2)*(n-3)*(n-4);
        n5=n*n14;

        a5=n^2*(n+5)/n14;
        a41=-n*(n+5)/n14;
        c41=5;
        a32=-n*(n-1)/n14;
        c32=10;
        a311=2*(n+2)/n14;
        c311=10;
        a221=2*(n-1)/n14;
        c221=15;
        a2111=-6/n14;
        c2111=10;
        a11111=24/n5;

        P1=sum(x);
        P2=sum(x.^2);
        P3=sum(x.^3);
        P4=sum(x.^4);
        P5=sum(x.^5);
        val=P5*a5+c41*a41*P4*P1+c32*a32*P3*P2+c311*a311*P3*P1^2+c221*a221*P2^2*P1+c2111*a2111*P2*P1^3+a11111*P1^5;
    otherwise
        error('The explicit solutions are only available up to order 5.');
end
end

function h=genHStat(x,p)
%%GENHSTAT Given a set of scalar samples, obtain an unbiased estimate of
%          the pth central sample moment. That is, estimate E{(x-mu)^p},
%          where E is the expected value operation and mu is the mean of
%          the distribution. Alternatively, p can be a length-v vector, and
%          an unbiased estimate of
%          E{(x-mu)^p(1)}*E{(x-mu)^p(2)}*...*E{(x-mu)^p(v)} is obtained.
%          When mu is not available and must be found from the samples, as
%          is the case here, then the average  mean((x-mean(x)).^p) is a
%          biased estimator. Additionally multiplying unbiased estimators
%          can lead to a biased estimator. The unbiased estimator for this
%          type of problem when p is a scalar is called an h-statistic and
%          when p is a vector is a generalized statistic.
%
%INPUTS: x An nX1 or 1Xn vector of real samples. Note that this function
%          can be slow for values of p>4 when a large number of samples is
%          used.
%        p The integer moment desired. p>=1. p can be a vector of all
%          unique elements if the product of central moments should be
%          estimated.
%
%OUTPUTS: h The unbiased estimate of the pth central moment given the
%           samples.
%
%The estimator of [1], which is an h-statistic, is rather complicated. An
%explicit expression for h-statistics in terms of symmetric means is given
%in [2], and an expression for generalized h statististics in terms of
%symmetric means is also provided. These are how this function is
%implemented. Equation 2.2 is used when p is a scalar (a central moment)
%and Equation 3.3 is used when p is a vector (a product of central
%moments). Note that the evaluation of symmetric means involves the
%evaluation of matrix permanents and thus the complexity scales
%exponentialy with the number of samples in x. Only small numbers of
%samples can be used, because the computational complexity increases too
%rapidly.
%
%REFERENCES:
%[1] P. R. Halmos, "The theory of unbiased estimation," The Annals of
%    Mathematical Statistics, vol. 17, no. 1, pp. 34-43, Mar. 1946.
%[2] D. S. Tracy and B. C. Gupta, "Generalized h-statistics and other
%    symmetric functions," The Annals of Statistics, vol. 2, no. 4, pp.
%    837-844, 1974.
%
%September 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

n=length(x);

if(any(p>=n))
    error('p must be >=n for an unbiased estimator to exist.')
end

if(isscalar(p))
    if(p==1)%Return the mean instead of 0.
        h=mean(x);
        return
    end

    h=0;
    for r=0:p
        h=h+(-1)^r*binomial(p,r)*symmetricMean(x,[p-r,ones(1,r)]);
    end
else
    v=length(p);

    p=p(:)';
    h=0;

    r=zeros(1,v);
    while(~isempty(r))
        sumR=sum(r);

        prodVal=1;
        for k=1:v
            prodVal=prodVal*binomial(p(k),r(k));
        end

        %We want p-r without the zero elements.
        pV=p-r;
        pV(pV==0)=[];

        h=h+(-1)^sumR*prodVal*symmetricMean(x,[pV,ones(1,sumR)]);

        r=getNextTuple(r,p);
    end
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 plotFit(min_x, max_x, mu, sigma, theta, p)
%PLOTFIT Plots a learned polynomial regression fit over an existing figure.
%Also works with linear regression.
%   PLOTFIT(min_x, max_x, mu, sigma, theta, p) plots the learned polynomial
%   fit with power p and feature normalization (mu, sigma).

% Hold on to the current figure
hold on;

% We plot a range slightly bigger than the min and max values to get
% an idea of how the fit will vary outside the range of the data points
x = (min_x - 15: 0.05 : max_x + 25)';

% Map the X values
X_poly = polyFeatures(x, p);
X_poly = bsxfun(@minus, X_poly, mu);
X_poly = bsxfun(@rdivide, X_poly, sigma);

% Add ones
X_poly = [ones(size(x, 1), 1) X_poly];

% Plot
plot(x, X_poly * theta, '--', 'LineWidth', 2)

% Hold off to the current figure
hold off

end

author: Ayatans

repo: Machine-Learning-homework