text:
% [INPUT]
% data = A float t-by-n matrix containing the time series to be sanitized.
% x = A vector of length t containing the numeric observation dates of the time series (if empty, observations are assumed to be linearly spaced between 1 and t).
% w = An integer [5,21] representing the length of the moving window used to detect outliers (if empty, no outliers replacement is performed).
% m = A vector of 2 floats (-Inf,Inf) containing minimum and maximum clamping values (if empty, no clamping is performed).
%
% [OUTPUT]
% data = A float t-by-n matrix containing the sanitized time series.

function data = sanitize_data(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('data',@(x)validateattributes(x,{'double'},{'real' '2d' 'nonempty'}));
        ip.addRequired('x',@(x)validateattributes(x,{'double'},{'real'}));
        ip.addRequired('w',@(x)validateattributes(x,{'double'},{'real'}));
        ip.addRequired('m',@(x)validateattributes(x,{'double'},{'real'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    data = ipr.data;
    [x,w,m] = validate_input(data,ipr.x,ipr.w,ipr.m);

    nargoutchk(1,1);

    data = sanitize_data_internal(data,x,w,m);

end

function data = sanitize_data_internal(data,x,w,m)

    for i = 1:size(data,2)
        y = data(:,i);

        nan_indices = isnan(y);

        if (all(nan_indices))
            continue;
        end

        if (any(nan_indices))
            y = fill_missing_values(y,x,nan_indices);
        end

        if (~isempty(w))
            y = replace_outliers(y,x,w);
        end

        if (~isempty(m))
            y = min(max(y,m(1)),m(2));
        end

        data(:,i) = y;
    end

end

function y = fill_missing_values(y,x,nan_indices)

    d = diff(nan_indices);

    if (nan_indices(1))
        z = find(d == -1,1,'first');
        y(1:z) = y(z+1);
        nan_indices(1:z) = 0;
    end

    if (nan_indices(end))
        z = find(d == 1,1,'last');
        y(z+1:end) = y(z);
        nan_indices(z+1:end) = 0;
    end

    y(nan_indices) = spline(x(~nan_indices),y(~nan_indices),x(nan_indices));

end

function y = replace_outliers(y,x,w)

    f = 3 * (-1 / (sqrt(2) * erfcinv(1.5)));
    k = floor(w * 0.5);
    n = numel(y);
    p = nan(k,1);

    xp = [p; y; p];

    m_med = zeros(n,1);
    m_mad = zeros(n,1);

    for i = 1:k
        x_i = y(1:k+i);
        m = median(x_i);

        m_med(i) = m;
        m_mad(i) = median(abs(x_i - m));
    end

    for i = k+1:n-k-1
        x_i = xp(i:i+w,:);
        m = median(x_i);

        m_med(i) = m;
        m_mad(i) = median(abs(x_i - m));
    end

    for i = n-k:n
        x_i = y(i-k:end);
        m = median(x_i);

        m_med(i) = m;
        m_mad(i) = median(abs(x_i - m));
    end

    b = m_mad .* f;
    lb = m_med - b;
    ub = m_med + b;

    is_outlier = (y > ub) | (y < lb);
    d = diff(is_outlier);

    if (is_outlier(1) == 1)
        z = find(d == -1,1,'first');
        y(1:z) = y(z+1);
        is_outlier(1:z) = 0;
    end

    if (is_outlier(end))
        z = find(d == 1,1,'last');
        y(z+1:end) = y(z);
        is_outlier(z+1:end) = 0;
    end

    y(is_outlier) = spline(x(~is_outlier),y(~is_outlier),x(is_outlier));

end

function [x,w,m] = validate_input(data,x,w,m)

    t = size(data,1);

    if (isempty(x))
        x = 1:t;
    else
        if (~isvector(x))
            error('The value of ''x'' is invalid. Expected input to be a vector.');
        end

        if (numel(x) ~= t)
            error(['The value of ''x'' is invalid. Expected input to contain ' num2str(t) ' elements.']);
        end

        if (~all(isfinite(x)))
            error('The value of ''x'' is invalid. Expected input to contain finite elements.');
        end

        if (~all(diff(x) > 0))
            error('The value of ''x'' is invalid. Expected input to contain increasing elements.');
        end
    end

    if (~isempty(w))
        if (~isscalar(w))
            error('The value of ''w'' is invalid. Expected input to be a scalar.');
        end

        if (~isfinite(w))
            error('The value of ''w'' is invalid. Expected input to be finite.');
        end

        if (floor(w) ~= w)
            error('The value of ''w'' is invalid. Expected input to be an integer.');
        end

        if ((w < 5) || (w > 21))
            error('The value of ''w'' is invalid. Expected input to have a value >= 5 and <= 21.');
        end

        w = w - 1;
    end

    if (~isempty(m))
        if (~isvector(m))
            error('The value of ''m'' is invalid. Expected input to be a vector.');
        end

        if (numel(m) ~= 2)
            error('The value of ''m'' is invalid. Expected input to contain 2 elements.');
        end

        if (~all(isfinite(m)))
            error('The value of ''m'' is invalid. Expected input to contain finite elements.');
        end

        if (any(floor(m) ~= m))
            error('The value of ''m'' is invalid. Expected input to contain integer elements.');
        end

        if (m(1) >= m(2))
            error('The value of ''m'' is invalid. Expected input first element to be less than the input second element.');
        end
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function varargout = quiver( F, varargin )
%QUIVER   Quiver plot of DISKFUNV.
%   QUIVER(F) plots the vector velocity field of F. QUIVER automatically
%   attempts to scale the arrows to fit within the grid.
%
%   QUIVER(F,S) automatically scales the arrows to fit within the grid and then
%   stretches them by S.  Use S=0 to plot the arrows without the automatic
%   scaling.
%
%   QUIVER(X,Y,F,...) is the same as QUIVER(F,...) except the arrows are on the
%   grid given in X and Y.
%
%   QUIVER(...,LINESPEC) uses the plot linestyle specified for the velocity
%   vectors.  Any marker in LINESPEC is drawn at the base instead of an arrow on
%   the tip.  Use a marker of '.' to specify no marker at all.  See PLOT for
%   other possibilities.
%
%   QUIVER(...,'numpts',N) plots arrows on a N by N grid.
%
%   H = QUIVER(...) returns a quivergroup handle.

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

numpts = 30;

% Empty check:
if ( isempty( F ) )
    quiver( [] )
    return
end

if ( isempty( varargin ) )
    varargin = {};
end

holdState = ishold;
if ( ~holdState )
    hold on;
end

if ( ~holdState )
    % Generate a unit disk
    N = 200;
    th = linspace(-pi, pi, N)';
    r = exp( 1i*th );
    plot(real(r), imag(r), 'k--');
end

% Number of points to plot
j = 1;
argin = {};
while ( ~isempty( varargin ) )
    if strcmpi( varargin{1}, 'numpts' )
        numpts = varargin{2};
        varargin(1:2) = [];
    else
        argin{j} = varargin{1};
        varargin(1) = [];
        j = j+1;
    end
end
varargin = argin;

if ( isa(F, 'diskfunv') )             % quiver(F,...)

        % Plot quiver with arrows at equally spaced points:
        [xx, yy] = diskpts(numpts);
        F1 = F.components{1};
        F2 = F.components{2};
        vals1 = feval(F1, xx, yy, 'cart');
        vals2 = feval(F2, xx, yy, 'cart');
        h = quiver(xx, yy, vals1,vals2, varargin{:});
        if ( ~holdState )
            axis tight;
            axis( max([max(abs(xlim)) abs(ylim)])*[-1 1 -1 1] );
            axis equal;
        end

elseif ( nargin >= 3 )                 % quiver(x,y,F,...)

    % First two arguments contain arrow locations: we assume these are
    % Cartesian coords.
    xx = F;
    yy = varargin{1};

    if ( isa(varargin{2}, 'diskfunv') )

        F = varargin{2};
            F1 = F.components{1};
            F2 = F.components{2};
            vals1 = feval(F1, xx, yy, 'cart');
            vals2 = feval(F2, xx, yy, 'cart');
            h = quiver( xx, yy, vals1, vals2, varargin{3:end} );
            if ( ~holdState )
                axis tight;
                axis(max([max(abs(xlim)) abs(ylim)])*[-1 1 -1 1]);
                axis equal;
            end

    else

        error('DISKFUN:DISKFUNV:quiver:inputs', ...
                                  'Third argument should be a diskfunv.');

    end

end

if ( ~holdState )
    hold off;
end

if ( nargout > 0 )
    varargout = {h};
end

end

% Generates a nice set of points on the unit disk that are roughly equally
% spaced.
function [xx,yy] = diskpts(numpts)
% Idea is to use a polar grid (r,theta), but instead of using the same
% number of points in theta for every r, we make it a function of r.  We
% start with the origin, then as r increases from the origin there will be
% 6, 9, 15, 21, ... until the outer r=1 radius is reached.  This appears to
% give a nice sampling of the disk.

% The number of radii to use to get approximately numpts is
n = floor( numpts/sqrt(3) );

% Increment for the radii
dr = 1/n;

% Add the origin
xx = 0;
yy = 0;

% Handle the second ring from the origin a bit differently (i.e. do 6
% points instead of 3).
th = trigpts(6,[-pi pi]);
xx = [xx; dr*cos(th)];
yy = [yy; dr*sin(th)];

% Add points 3*(2*k-1) points for radii k.
for k = 2:n
    th = trigpts(3*(2*k-1), [-pi pi]);
    xx = [xx ; dr*k*cos(th)];
    yy = [yy ; dr*k*sin(th)];
end
end

% ANOTHER OPTION, BUT I DON'T THINK IT'S AS NICE.
% % Generates a nice set of points on the unit disk using the technique
% % described in Section 3.3 of
% % D. Calhoun, C. Helzel, R. J. LeVeque. Logically rectangular grids and
% % finite volume methods for PDEs in circular and spherical domains. SIAM
% % Review Vol. 50, Issue 4, pp. 723-752. (2008)
% function [xx,yy] = diskpts(numpts)
%
% % Generate equally spaced points over [-1,1]x[-1,1].  These will be mapped
% % to the unit disk according the the algorithm given in Section 3.3. of the
% % paper referenced above.
% [xc, yc] = meshgrid(linspace(-1,1,numpts));
%
% r1 = 1;   % map [-1,1] x [-1,1] to circle of radius r1
% d = max(abs(xc),abs(yc));
% r = sqrt(xc.^2 + yc.^2);
% r = max(r,1e-10);
% xx = r1 * d .* xc./r;
% yy = r1 * d .* yc./r;
% w = d.^2;
% xx = w.*xx + (1-w).*xc/sqrt(2);
% yy = w.*yy + (1-w).*yc/sqrt(2);
%
% end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function  a = i4vec_indicator0 ( n )

%*****************************************************************************80
%
%% I4VEC_INDICATOR0 sets an I4VEC to the indicator vector (0,1,2,...).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Output, integer A(N), the vector.
%
  a = ( 0 : n - 1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function aux_plotvolcano(params, hParentFigure)
% function aux_FMD(params, hParentFigure);
%-------------------------------------------
%
% Incoming variables:
% params        : all variables
% hParentFigure : Handle of the parent figure
%
% Thomas van Stiphout, thomas@sed.ethz.ch
% last update: 7.9.2005

report_this_filefun(mfilename('fullpath'));
hold on;plot3(-150,62,0,'r^','MarkerSize',15)
hold on;plot3([-150 -150],[62 62],[0 -200],'k','Linewidth',4)

author: CelsoReyes


text:
%% Check rate of convergence for 3D Hcurl interface problem
%     curl(A curl u) + B u  = f,    x\in \Omega
%      where A and B are piecewise constants on Omega^+ and Omega^-.
%
% Domain: Rectangular domain: [xmin,xmax] X [ymin,ymax] X [zmin,zmax]
% Mesh: Cartesian triangular mesh.
% Method: FE
% for this method, need to mannualy import a mesh data

%% Geometry and Boundary Conditions
%clear
%close all
%clc

%path(pathdef)
addpath(genpath(pwd),'-begin');
rmpath(genpath('./.git'));
rmpath(genpath('./docs'));
savepath;

domain = [-1,1,-1,1,-1,1];
bc = [1,1,1,1,1,1]; % Dirichelet BC

%% Finite Element Type
femtype = '1st kind Nedelec';
disp(['FEM Type =  Conforming ', femtype]);

%% Initial Partition
nx0 = 10;
ny0 = nx0;
nz0 = nx0;

%% Task
showErr = 0;
showMesh = 0;
computErr = 1;

%% PDE
test = 3;
switch test
    case 0
        pde = poissonNonPoly3D;
    case 1 % circular interface
        r = pi/5; bm = 1; bp = 1; am = 1; ap = 1;
        x0 = 0; y0 = 0; z0 = 0; a11 = 1; a12 = 1; a = 1;
        pde = elli3DcircIntf2(am,ap,bm,bp,r,x0,y0,z0,a11,a12,a);
    case 2
        r = pi/5; bm = 1; bp = 1;
        x0 = 0; y0 = 0; z0 = 0; coef = 1;
        pde = constant_fun(bm,bp,r,x0,y0,z0,coef);
    case 3
        r = pi/5; bm = 1; bp = 1; am = 1; ap = 1;
        x0 = 0; y0 = 0; z0 = 0; coef = 1;
        pde = linear_fun(am,ap,bm,bp,r,x0,y0,z0,coef);
end

%% Max Iteration

%% 1. Generate Mesh
disp('*******************************************************************************')
disp('Import mesh');
load('BoxTorusMesh.mat')
mesh = enrichMesh3D(mesh,0); % Mesh detail level = 1 (for IFE).

%% 2. Generate FEM DoF
fem = genNedFEM3D(mesh,bc);
disp(['number of DoF =  ', int2str(length(fem.p))]);

%% 3. Assemble Matrix
disp(' '); disp('Start Assembling Matrix');
S = globMatrixNedFit3D(am,ap,1,mesh,fem,fem);
M = globMatrixNedFit3D(bm,bp,0,mesh,fem,fem);
rhsF1 = globNedFitRHS3D(pde.fm1,pde.fp1, mesh, fem, 0, 1);
rhsF2 = globNedFitRHS3D(pde.fm2,pde.fp2, mesh, fem, 0, 2);
rhsF3 = globNedFitRHS3D(pde.fm3,pde.fp3, mesh, fem, 0, 3);
rhsF = rhsF1 + rhsF2 + rhsF3;
Atotal = S + M;
tu1 = sum(feval(pde.exactu1,fem.gex,fem.gey,fem.gez).*fem.gew,2);
tu2 = sum(feval(pde.exactu2,fem.gex,fem.gey,fem.gez).*fem.gew,2);
tu3 = sum(feval(pde.exactu3,fem.gex,fem.gey,fem.gez).*fem.gew,2);
tgt = mesh.p(mesh.e(:,2),:) - mesh.p(mesh.e(:,1),:);
tgt = tgt./sum(tgt.^2,2).^(1/2);
tu = tu1.*tgt(:,1) + tu2.*tgt(:,2) + tu3.*tgt(:,3);

Ndof = size(mesh.e,1);
bdidx = zeros(Ndof,1);
isBdEdge = true(Ndof,1);
isBdEdge(fem.mapper) = false;
bdidx(isBdEdge) = 1;
Tbd = spdiags(bdidx,0,Ndof,Ndof);
T = spdiags(1-bdidx,0,Ndof,Ndof);
A = T*Atotal*T + Tbd;
ub = tu;
ub(fem.mapper) = 0;
rhsB = Atotal*ub;
f = rhsF - rhsB;
f(isBdEdge) = tu(isBdEdge);

%% 3. Solve the linear system Au = f
option.outsolver = 'cg';
tID1 = (mesh.tLoc == 1);
tID2 = (mesh.tLoc == 2);
e1tmp = unique(reshape(mesh.t_e(tID1,:),[],1));
e2tmp = unique(reshape(mesh.t_e(tID2,:),[],1));
eInt = intersect(e1tmp,e2tmp);
%eIntp1 = mesh.p(mesh.e(eInt,1),:);
%eIntp2 = mesh.p(mesh.e(eInt,2),:);
eid1 = setdiff(e1tmp,eInt);
eid2 = setdiff(e2tmp,eInt);
alpha = am*ones(size(mesh.e,1),1);
alpha(eid2) = ap;
alpha(eInt) = (am+ap)/2;
beta = bm*ones(size(mesh.e,1),1);
beta(eid2) = bp;
beta(eInt) = (bm+bp)/2;
option.alpha = alpha;
option.beta = beta;
option.solver = 'amg';
edge = mesh.e;
option.isBdEdge = isBdEdge;
option.smoother = 'BD';
NEdof = size(A,1);
option.blklevel = 0;
option.blkId = NEdof;
option.fact = 'chol';
[x,info] = amgMaxwellinterface(A,f,mesh.p,edge,option);
uh = x;

%% 4. Postprocess: Calculating Errors
tic
errND = max(abs(uh - tu)); % Error on nodes
err.nd = errND; err.inf = 0; err.l2 = 0; err.h1 = 0;
if computErr == 1
    disp(' ');  disp('Start computing error in L2 norm');
    eNorm = 'L2'; disp(['Start computing error in ',eNorm,'  norm']);
    [errL2,errL2K1,errL2K2,errL2K3] = getCurlErr3D(uh, pde, fem, eNorm);

    eNorm = 'Curl'; disp(['Start computing error in ',eNorm,' norm']);
    [errCurl,errCurl1,errCurl2,errCurl3] = getCurlErr3D(uh, pde, fem, eNorm);
    err.l2 = errL2;
    err.curl = errCurl;
end
time(6) = toc;
time(7) = 1e6*sum(time(2:6))/length(mesh.t);

%% 5: Output

disp(' ')
disp('Errors')
disp('Node   L2 norm     H1 norm')
formatSpec = '%6.4e %6.4e  %6.4e\n';
fprintf(formatSpec, err.nd, err.l2, err.curl)

err0 = err; h0 = h;

disp(' '); disp('CPU Time')
disp('   N     Mesh     FEM      Matrix   Solve    Error    Time/1M cell')
formatSpec = '%4i  %7.2f  %7.2f  %7.2f  %7.2f  %7.2f   %7.2f\n';
fprintf(formatSpec, time)

%% 6. Plot Solution and Error
if showErr == 1
end

repo: ifem


text:
function [targetMissDist, dVNorm, hOrbit, dVVect, dVVectNTW, eRVect, xferOrbit, departUT, timeSOI] = departureTargetMissDistance(x, dVVectNTW, eOrbit, departBodyInfo, arriveBodyInfo, parentBodyInfo, departUT, arrivalUT)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
    eTru = x(1);
    progradeDV = x(2);
    dVVectNTW(1) = progradeDV;

    eSMA = eOrbit(1);
    eEcc = eOrbit(2);
    eInc = eOrbit(3);
    eRAAN = eOrbit(4);
    eArg = eOrbit(5);
%     eTru = eOrbit(6);

    eTruAtDepartUT = eOrbit(6);
    eMeanMotion = computeMeanMotion(eSMA, departBodyInfo.gm);

    eMean = computeMeanFromTrueAnom(eTru, eEcc);
    eMeanAtDepartUT = computeMeanFromTrueAnom(eTruAtDepartUT, eEcc);
    departDeltaTBurnPosAdjust = (eMean-eMeanAtDepartUT)/eMeanMotion;
    departUT = departUT + departDeltaTBurnPosAdjust;

    [eRVect,eVvect] = getStatefromKepler(eSMA, eEcc, eInc, eRAAN, eArg, eTru, departBodyInfo.gm);
    tHat = eVvect/norm(eVvect);
    wHat = cross(eRVect,eVvect)/norm(cross(eRVect,eVvect));
    nHat = cross(tHat,wHat)/norm(cross(tHat,wHat));
    ECI2TWNRotMat = [tHat,wHat,nHat];
    dVVect = ECI2TWNRotMat * dVVectNTW;
    dVNorm = norm(dVVect);

    hRVect = eRVect;
    hVVect = eVvect + dVVect;
    [hSMA, hECC, hINC, hRAAN, hARG, hTRU] = getKeplerFromState(hRVect,hVVect,departBodyInfo.gm);
    hOrbit = [hSMA, hECC, hINC, hRAAN, hARG, hTRU];

    rSOI = getSOIRadius(departBodyInfo, parentBodyInfo);

    hTruSOI = computeTrueAFromRadiusEcc(rSOI, hOrbit(1), hOrbit(2));
    [rVectSOI,vVectSOI]=getStatefromKepler(hOrbit(1), hOrbit(2), hOrbit(3), hOrbit(4), hOrbit(5), hTruSOI, departBodyInfo.gm);
    hOrbit(7) = hTruSOI;

    meanDepartTime = computeMeanFromTrueAnom(hOrbit(6), hOrbit(2));
    meanSOITime = computeMeanFromTrueAnom(hTruSOI, hOrbit(2));
    hMeanMotion = computeMeanMotion(hOrbit(1), departBodyInfo.gm);
    deltaT2SOI = (meanSOITime-meanDepartTime)/hMeanMotion;

    timeSOI = departUT - deltaT2SOI; %may be a + sign here
    [departBodyRVect, departBodyVVect] = getStateAtTime(departBodyInfo, timeSOI, parentBodyInfo.gm);
    scInParentOrbitRVect = departBodyRVect + rVectSOI;
    scInParentOrbitVVect = departBodyVVect + vVectSOI;
    [smaXAct, eccXAct, incXAct, raanXAct, argXAct, truDXAct] = getKeplerFromState(scInParentOrbitRVect,scInParentOrbitVVect,parentBodyInfo.gm);
    xferOrbit = [smaXAct, eccXAct, incXAct, raanXAct, argXAct, truDXAct];

    bodyInfoXAct.sma = smaXAct;
    bodyInfoXAct.ecc = eccXAct;
    bodyInfoXAct.inc = rad2deg(incXAct);
    bodyInfoXAct.raan = rad2deg(raanXAct);
    bodyInfoXAct.arg = rad2deg(argXAct);
    bodyInfoXAct.mean = computeMeanFromTrueAnom(truDXAct, eccXAct);
    bodyInfoXAct.epoch = timeSOI;

    [rVectXAct, vVectXAct] = getStateAtTime(bodyInfoXAct, arrivalUT, parentBodyInfo.gm);
    [~, ~, ~, ~, ~, truAXAct] = getKeplerFromState(rVectXAct,vVectXAct,parentBodyInfo.gm);
    xferOrbit(7) = truAXAct;
    [rVectArriveBody, ~] = getStateAtTime(arriveBodyInfo, arrivalUT, parentBodyInfo.gm);
    targetMissDist = norm(rVectArriveBody-rVectXAct);
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function g=aicplot(alpha,varargin)
%
% Computes and plots the -2*AIC
% for local fits with different smoothing parameters.
%
% The first argument to aicplot(), alpha, should be a matrix with one
% or two columns (first column = nearest neighbor component, second
% column = constant component). Each row of this matrix is, in turn,
% passed as the 'alpha' argument to aic() (and locfit()). The results
% are stored in a matrix, and aic score ploted against the degrees of
% freedom.

k = size(alpha,1);
z = zeros(k,4);

for i=1:k
  z(i,:) = aic(varargin{:},'alpha',alpha(i,:));
end;

plot(z(:,3),z(:,4));
xlabel('Fitted DF');
ylabel('AIC');

g = [alpha z];
return;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function out = ir_dct8(x, adj)
%| 2D DCT of each 8x8 block
%| set adj=1 for adjoint (transpose)

if nargin < 1, ir_usage, end
if nargin < 2, adj = false; end

if streq(x, 'test'), ir_dct8_test, return, end

dc = ir_dctmtx(8);
if adj
	dc = dc';
end

out = zeros(size(x));
[nx ny] = size(x);
for by=1:floor(ny/8)
	for bx=1:floor(nx/8)
		ix = (bx-1)*8+(1:8);
		iy = (by-1)*8+(1:8);
		block = x(ix,iy);
		tmp = (dc * block) * dc';
		out(ix,iy) = tmp;
	end
end

end

function ir_dct8_test
	nx = 32;
	ny = 16;
	X = rand(nx,ny);
	d1 = ir_dct8(X);
%	im(d1)

	b8 = dctmtx(8);
	bx = kron(eye(nx/8), b8);
	by = kron(eye(ny/8), b8);
	d2 = bx * X * by';
%	im(d2)
	assert(max_percent_diff(d1, d2) < 1e-12)

	d3 = kron(by, bx) * X(:);
	d3 = reshape(d3, nx, ny);
%	im(d3)
	assert(max_percent_diff(d1, d3) < 1e-12)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function ll = gpDynamicsPointLogLikelihood(model, x, y)

% GPDYNAMICSPOINTLOGLIKELIHOOD Compute the log likelihood of a point under the GP dynamics model.
% DESC computes the log likelihood of a given point under the GP
% dynamics model.
% ARG model : the model for which the log likelihood is to be
% computed.
% ARG x : the input location for the model.
% ARG y : the target location for the model.
% RETURN ll : the log likelihood of the given point.
%
% COPYRIGHT : Neil D. Lawrence, 2006
%
% SEEALSO : gpDynamicsCreate, gpDynamicsLogLikelihood, gpPointLogLikelihood
%

% MODIFICATION: Carl Henrik Ek, 2009

% FGPLVM

if(isfield(model,'indexIn'))
  x = x(:,model.indexIn);
  y = y(:,model.indexIn);
end
ll = gpPointLogLikelihood(model, x, y);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%%Ex. 2 The meaning of "a = b"

%In Matlab and in any programming language, the statement "a = b" does not mean
%"a equals b". Instead, it prompts the action of replacing the content of a by the
%content of b.
a = 3;
b = a;
b

%Output:3

author: TheAlgorithms

repo: MATLAB-Octave


text:
function [DVal,totalError,exitCode]=KLDistApproxGaussMix(w1,mu1,P1,w2,mu2,P2,algorithm,param8,w)
%%KLDISTAPPROXGAUSSMIX Approximate the Kullback-Leibler (KL) divergence
%          (also known as the KL distance or relative entropy) between two
%          multivariate Gaussian mixture distributions. Whereas an explicit
%          solution is available for the KL distance between two Gaussian
%          PDFs, no explicit solution exists between a Gaussian PDF and a
%          Gaussian mixture. Thus, various approximations are necessary.
%          As defined in Chapter 8.5 of [4], the relative entropy between
%          two continuous PDFs, f and g, is
%          D(f||g)=integral_{x} f(x) log(f(x)/g(x)) dx
%          The integral is over an appropriate domain of x, in this case
%          the d-dimensional real space for an n-dimensional state. Note
%          that changing the order of f and g changes the result (the
%          operation is not commutative).
%
%INPUTS: w1 The 1Xn1 or n1X1 set of weights of the first Gaussian mixture's
%           components such that all w>=0 and sum(w)=1.
%       mu1 The dXn1 set of mean values of the first Gaussian mixture's
%           components.
%        P1 The dXxdXn1 set of covariance matrices of the first Gaussian
%           mixture's components.
% w2,mu2,P2 The equivalent of w1, mu1 and P1 for the second mixtures
%           components. There length-n2, dXn2, and dXdXn2 in size.
% algorithm An optional parameter specifying the algorithm used to obtain
%           the approximation of the KL divergence. Possible values are
%           for:
%           0 Use adaptive numeric integration to directly evaluate the
%             Kullback-Leiberler divergence integral over a rectangular
%             region that is standard deviations from the means of the
%             distributions. The default number of standard deviations in
%             each direction is 4.5 and can be adjusted by passing a member
%             of param8 named 'numStdDev'. In 1D, the function
%             integral1DAdaptiveCC is used for the integration and the
%             inputs to that function 'RelTol', 'AbsTol', 'NPowMin', and
%             'NPowMax' can be passed as members of param8 if one wishes to
%             not use the defaults. If higher dimensions, the function
%             integrateUnifCub57Adaptive is used for the numeric
%             integration and the inputs to that function 'maxSearchReg',
%             'AbsTol', and 'RelTol' can be passed if one wishes to not use
%              the defaults. This approximation can produce negative
%              values, though it is often good.
%           1 Monte Carlo sampling from Equation 4 of Section 2 of [3].
%             This uses param8 as the number of samples (if given). The
%             default if param8 is not given is 3e3. As the number of
%             samples increases, this method should become increasingly
%             accurate. This approximation can produce negative values.
%           2 (The default if omitted or an empty matrix is passed) The
%             cubature integration approximation from Equation 8 in Section
%             3 of [3]. This uses param8 as the cubature points and w as
%             the weights, if provided. If not provided, then
%             fifthOrderCubPoints is used to obtain fifth-order xDim-
%             dimensional cubature points. This approximation can produce
%             negative values.
%           3 Use the variational upper bound from Section 8 of [3]. This
%             approximation should always be positive. 10 iterations are
%             performed. This bound is typically not zero when the two
%             input PDFs are equal.
%           4 The variational approximation from Section 7 of [3] (also
%             given in [1]). This approximation can produce negative
%             values.
%           5 The product of Gaussians approximation in [1] and Section 5
%             of [3]. This approximation can produce negative values.
%           6 Take the average of algorithms 5 and 4, which tend to be
%             upper and lower bounds.
%           7 Goldberger's approximation from Equation 2 of [2]. This is
%             also given in Section 6 of [3] as the matched bound
%             approximation. This approximation can produce negative
%             values and is often quite bad with certain classes of
%             distributions.
%    param8 An optional parameter. If algorithm=0, then this is a structure
%           that can take members as described above for algorithm 0. If
%           algorithm=1, then this is the number of samples to use. If
%           algorithm=2, then this is a set of cubature points as a
%           dXnumPoints matrix. For other values of algorithm, this input
%           is not used.
%         w This input is only used if algorithm=2, in which case it is the
%           set of cubature weights (sum(w)=1). that are associated with
%           the cubature points in param8. If param8 is omitted or an empty
%           matrix is passed, then w will be the weights associated with
%           the points from the fifthOrderCubPoints function.
%
%%OUTPUTS: val The approximate value of the Kullback-Leibler divergence
%              between the first and the second Gaussian mixture
%              distributions.
%
%Some approximations can produce negative values (the true KL divergence
%can never be negative).
%
%EXAMPLE 1:
%Though many of the algorithms are often cited in the literature, many of
%the approximations can be quite bad. Here, we consider a 1D example and
%then the mixture missing a few components and also the moment-matched
%Gaussian approximation.
% w1=[0.03,0.18,0.12,0.19,0.02,0.16,0.06,0.1,0.08,0.06];
% n1=length(w1);
% mu1=[1.45,2.20,0.67,0.48,1.49,0.91,1.01,1.42,2.77,0.89];
% P1=[0.0487,0.0305,0.1171,0.0174,0.0295,0.0102, 0.0323, 0.0380, 0.0115, 0.0679];
% P1=reshape(P1,[1,1,n1]);
%
% %The second PDF is the first with the five least-weight components deleted.
% w2=[0.18,0.12,0.19,0.16,0.1,0.08];
% w2=w2/sum(w2);
% n2=length(w2);
% mu2=[2.20,0.67,0.48,0.91,1.42,2.77];
% P2=[0.0305,0.1171,0.0174,0.0102,0.0380,0.0115];
% P2=reshape(P2,[1,1,n2]);
%
% %The third PDF is the second PDF with two more components deleted.
% w3=[0.18,0.12,0.19,0.16];
% w3=w3/sum(w3);
% n3=length(w3);
% mu3=[2.20,0.67,0.48,0.917];
% P3=[0.0305,0.1171,0.0174,0.0102];
% P3=reshape(P3,[1,1,n3]);
%
% %The fourth PDF is just the Gaussian matching the first two moments of the
% %original PDF.
% w4=1;
% [mu4,P4]=calcMixtureMoments(mu1,w1,P1);
%
% DVal=zeros(8,3);
% for curAlg=0:7
%     DVal(curAlg+1,1)=KLDistApproxGaussMix(w1,mu1,P1,w2,mu2,P2,curAlg);
%     DVal(curAlg+1,2)=KLDistApproxGaussMix(w1,mu1,P1,w3,mu3,P3,curAlg);
%     DVal(curAlg+1,3)=KLDistApproxGaussMix(w1,mu1,P1,w4,mu4,P4,curAlg);
% end
% RelErr=abs(bsxfun(@rdivide,bsxfun(@minus,DVal,DVal(1,:)),DVal(1,:)))
% %The rows of RelErr are the algorithms and the PDF chosen for comparison
% %against the original PDF is given by the column. The value obtained
% %through the numeric integration (algorithm 0) is nearly exact, so it is
% %taken as a basis of comparison for the relative error of the
% %approximations. It can be seen that some approximations are much worse
% %than others.
% %We can visualize the above PDFs.
% numPoints=500;
% xVals=linspace(0,3,numPoints);
% PDFVals1=GaussianMixtureD.PDF(xVals,w1,mu1,P1);
% PDFVals2=GaussianMixtureD.PDF(xVals,w2,mu2,P2);
% PDFVals3=GaussianMixtureD.PDF(xVals,w3,mu3,P3);
% PDFVals4=GaussianMixtureD.PDF(xVals,w4,mu4,P4);
% figure(1)
% clf
% hold on
% plot(xVals,PDFVals1,'-k','linewidth',4)
% plot(xVals,PDFVals2,'-r','linewidth',2)
% plot(xVals,PDFVals3,'--g','linewidth',2)
% plot(xVals,PDFVals4,'-.b','linewidth',2)
% legend('10 Elements', '5 Elements', '3 Elements','Matched Gaussian')
%
%EXAMPLE 2:
%Here is an example with bivariare Gaussian mixtures.
% w1=[0.25;0.5;0.25];
% mu1=zeros(2,2);
% mu1(:,1)=[1;-1];
% mu1(:,2)=[-1;1];
% mu1(:,3)=[0;0];
% P1=zeros(2,2,2);
% P1(:,:,1)=[4/9,  14/45;
%            14/45,4/9];
% P1(:,:,2)=[4/9, 0;
%            0, 4/9];
% P1(:,:,3)=[2/9, -1/9;
%           -1/9, 3/9];
%
% %The second distribution just throws out the first component.
% w2=w1(2:3);
% w2=w2/sum(w2);
% mu2=mu1(:,2:3);
% P2=P1(:,:,2:3);
%
% %The third distribution replaces the original distribution with a single
% %moment-matched Gaussian.
% w3=1;
% [mu3,P3]=calcMixtureMoments(mu1,w1,P1);
%
% DVal=zeros(8,2);
% for curAlg=0:7
%     DVal(curAlg+1,1)=KLDistApproxGaussMix(w1,mu1,P1,w2,mu2,P2,curAlg);
%     DVal(curAlg+1,2)=KLDistApproxGaussMix(w1,mu1,P1,w3,mu3,P3,curAlg);
% end
% RelErr=abs(bsxfun(@rdivide,bsxfun(@minus,DVal,DVal(1,:)),DVal(1,:)))
% %The rows of RelErr are the algorithms and the PDF chosen for comparison
% %against the original PDF is given by the column. The value obtained
% %through the numeric integration (algorithm 0) has a few digits of
% %accuracy, so it is taken as a basis of comparison for the relative
% %error of the approximations. It can be seen that some approximations are
% %much worse than others.
% %We can visualize the above PDFs.
% figure(1)
% clf
% subplot(3,1,1);
% numPoints=250;
% vals=linspace(-3,3,numPoints);
% [X,Y]=meshgrid(vals,vals);
% points=[X(:)';Y(:)'];
% PDFVals=GaussianMixtureD.PDF(points,w1,mu1,P1);
% PDFVals=reshape(PDFVals,numPoints,numPoints);
% surf(X,Y,PDFVals,'EdgeColor','none')
% title('Full PDF')
% %Plot the reduced PDF
% subplot(3,1,2);
% numPoints=250;
% vals=linspace(-3,3,numPoints);
% [X,Y]=meshgrid(vals,vals);
% points=[X(:)';Y(:)'];
% PDFVals=GaussianMixtureD.PDF(points,w2,mu2,P2);
% PDFVals=reshape(PDFVals,numPoints,numPoints);
% surf(X,Y,PDFVals,'EdgeColor','none')
% title('Missing 1 Component')
% %Plot the Gaussian approximation.
% subplot(3,1,3);
% numPoints=250;
% vals=linspace(-3,3,numPoints);
% [X,Y]=meshgrid(vals,vals);
% points=[X(:)';Y(:)'];
% PDFVals=GaussianMixtureD.PDF(points,w3,mu3,P3);
% PDFVals=reshape(PDFVals,numPoints,numPoints);
% surf(X,Y,PDFVals,'EdgeColor','none')
% title('Gaussian Approximation')
%
%REFERENCES:
%[1] J. L. Durrieu, J. P. Thiran, and F. Kelly, "Lower and upper bounds for
%    approximation of Kullback-Leibler divergence between Gaussian mixture
%    models," in Proceedings of the IEEE International Conference on
%    Acoustics, Speech and Signal Processing, Kyoto, Japan, 25-30 Mar.
%    2012, pp. 4833-4836.
%[2] J. Goldberger, S. Gordon, and H. Greenspan, "An efficient image
%    similarity measure based on approximations of KL-divergence between
%    two Gaussian mixtures," in Proceedings of the Ninth IEEE International
%    Conference on Computer Vision, Nice, France, 2003.
%[3] J. R. Hershey and P. A. Olden, "Approximating the Kullback Leibler
%    divergence between Gaussian mixture models," in Proceedings of the
%    IEEE International Conference on Acoustics, Speech and Signal
%    Processing, Honolulu, HI, 15-20 Apr. 2007.
%[4] T. M. Cover and J. A. Thomas, Elements of Information Theory, 2nd ed.
%    Hoboken, NJ: Wiley-Interscience, 2006.
%
%July 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<7||isempty(algorithm))
    algorithm=2;
end

d=size(mu1,1);
n1=length(w1);
n2=length(w2);

totalError=[];
exitCode=0;
switch(algorithm)
    case 0%Numeric integration over a bounded region.
        if(nargin>7&&~isempty(param8))
            if(isfield(param8,'numStdDev'))
                numStdDev=param8.numStdDev;
            end
        else
            numStdDev=4.5;
        end

        [mu1Merged,P1Merged]=calcMixtureMoments(mu1,w1,P1);
        [mu2Merged,P2Merged]=calcMixtureMoments(mu2,w2,P2);

        stdDevStep1=numStdDev*sqrt(diag(P1Merged));
        stdDevStep2=numStdDev*sqrt(diag(P2Merged));

        P1Inv=zeros(d,d,n1);
        P1InvDet=zeros(n1,1);
        P2Inv=zeros(d,d,n2);
        P2InvDet=zeros(n2,1);
        for k=1:n1
            P1Inv(:,:,k)=inv(P1(:,:,k));
            P1InvDet(k)=det(P1Inv(:,:,k));
        end
        for k=1:n2
            P2Inv(:,:,k)=inv(P2(:,:,k));
            P2InvDet(k)=det(P2Inv(:,:,k));
        end

        upperBounds=max(mu1Merged+stdDevStep1,mu2Merged+stdDevStep2);
        lowerBounds=min(mu1Merged-stdDevStep1,mu2Merged-stdDevStep2);

        f=@(x)KLArgGaussMix(x,w1,mu1,P1Inv,P1InvDet,w2,mu2,P2Inv,P2InvDet);

        if(d==1)
            %A different integration routine is used for 1D veruss
            %arbitrary dimensional integrals.
            bounds=[lowerBounds;upperBounds];

            RelTol=[];
            AbsTol=[];
            NPowMin=[];
            NPowMax=[];
            if(nargin>7&&~isempty(param8))
                if(isfield(param8,'RelTol'))
                    RelTol=param8.RelTol;
                end

                if(isfield(param8,'AbsTol'))
                    AbsTol=param8.AbsTol;
                end

                if(isfield(param8,'NPowMin'))
                    NPowMin=param8.NPowMin;
                end

                if(isfield(param8,'NPowMax'))
                    NPowMax=param8.NPowMax;
                end
            end

            [DVal,totalError,exitCode]=integral1DAdaptiveCC(f,bounds,RelTol,AbsTol,NPowMin,NPowMax);
        else
            maxSearchReg=[];
            AbsTol=[];
            RelTol=[];
            if(nargin>7&&~isempty(param8))
                if(isfield(param8,'RelTol'))
                    RelTol=param8.RelTol;
                end

                if(isfield(param8,'AbsTol'))
                    AbsTol=param8.AbsTol;
                end

                if(isfield(param8,'maxSearchReg'))
                    maxSearchReg=param8.maxSearchReg;
                end
            end

            [DVal,totalError,exitCode]=integrateUnifCub57Adaptive(f,lowerBounds,upperBounds,maxSearchReg,AbsTol,RelTol);
        end
    case 1%Monte Carlo sampling from Equation 4 of Section 2 of [3].
        if((nargin<8||isempty(param8)))
            numSamp=2e3;
        else
            numSamp=param8;
        end

        S1=zeros(d,d,n1);
        S2=zeros(d,d,n2);
        for k=1:n1
            S1(:,:,k)=chol(P1(:,:,k),'lower');
        end
        for k=1:n2
            S2(:,:,k)=chol(P2(:,:,k),'lower');
        end

        x1Samp=GaussianMixtureD.randS(numSamp,w1,mu1,S1);
        PDFVals1=GaussianMixtureD.PDFS(x1Samp,w1,mu1,S1);

        logVals=log(PDFVals1)-log(GaussianMixtureD.PDFS(x1Samp,w2,mu2,S2));
        %Get rid of points where PDFVals1==0 due to finite precision
        %errors. This probably won't occur. However, eliminating them also
        %gets rid of any unlikely points where both PDFs evaluate to 0 and
        %thus avoids returning NaN values.
        sel=(PDFVals1==0);
        logVals(sel)=0;

        DVal=sum(logVals)/numSamp;
    case 2%The cubature integration approximation from Equation 8 in
          %Section 3 of [3].
        if(nargin<8||isempty(param8))
            [xi,w]=fifthOrderCubPoints(d);
        else
            xi=param8;
        end

        S1=zeros(d,d,n1);
        S2=zeros(d,d,n2);
        for k=1:n1
            S1(:,:,k)=chol(P1(:,:,k),'lower');
        end
        for k=1:n2
            S2(:,:,k)=chol(P2(:,:,k),'lower');
        end

        w=w(:).';
        DVal=0;
        for a=1:n1
            %The component over which the expected value is taken.
            xiCur=transformCubPoints(xi,mu1(:,a),S1(:,:,a));

            logVals=log(GaussianMixtureD.PDFS(xiCur,w1,mu1,S1))-log(GaussianMixtureD.PDFS(xiCur,w2,mu2,S2));
            %Get rid of NaN values.
            %logVals(isnan(logVals))=0;

            DVal=DVal+w1(a)*sum(logVals.*w);
        end
    case 3%The variational upper bound from Section 8 of [3].
        psi=zeros(n1,n2);
        phi=zeros(n2,n1);

        %We start with psi and phi being the same and equal to what would
        %produce the convexity bound:
        for a=1:n1
            for b=1:n2
                psi(a,b)=w1(a)*w2(b);
                phi(b,a)=psi(a,b);
            end
        end

        %Get all pairwise KL-divergences.
        DKLg=zeros(n1,n2);
        for a=1:n1
            for b=1:n2
                DKLg(a,b)=KLDistGauss(mu1(:,a),P1(:,:,a),mu2(:,b),P2(:,:,b));
            end
        end
        expnDKLg=exp(-DKLg);

        numIter=10;
        for curIter=1:numIter
            %Equation 24 for phi.
            for a=1:n1
                denom=sum(psi(a,:).*expnDKLg(a,:));
                for b=1:n2
                    phi(b,a)=w1(a)*psi(a,b)*expnDKLg(a,b)/denom;
                end
            end

            %Equation 23 for psi
            for b=1:n2
                denom=sum(phi(b,:));
                for a=1:n1
                    psi(a,b)=w2(b)*phi(b,a)/denom;
                end
            end
        end

        %The first term in Equation 22
        term1=0;
        term2=0;
        for a=1:n1
            for b=1:n2
                val=phi(b,a)*log(phi(b,a)/psi(a,b));
                if(isnan(val))
                    val=0;
                end
                term1=term1+val;

                term2=term2+phi(b,a)*DKLg(a,b);
            end
        end

        DVal=term1+term2;
    case 4%The variational approximation from Section 2.3 of [1] and
          %Section 7 of [3].

        %Compute the necessary KL divergences.
        DKLf=zeros(n1,n1);
        DKLg=zeros(n1,n2);
        for a=1:n1
            DKLf(a,a)=KLDistGauss(mu1(:,a),P1(:,:,a),mu1(:,a),P1(:,:,a));
            for ap=(a+1):n1
                DKLf(a,ap)=KLDistGauss(mu1(:,a),P1(:,:,a),mu1(:,ap),P1(:,:,ap));
                DKLf(ap,a)=KLDistGauss(mu1(:,ap),P1(:,:,ap),mu1(:,a),P1(:,:,a));
            end

            for b=1:n2
                DKLg(a,b)=KLDistGauss(mu1(:,a),P1(:,:,a),mu2(:,b),P2(:,:,b));
            end
        end

        %Evaluate the sum in Equation 18 of [1], which is the same as
        %Equation 20 of [3].
        DVal=0;
        w1=w1(:).';
        w2=w2(:).';
        for a=1:n1
            DVal=DVal+w1(a)*(log(sum(w1.*exp(-DKLf(a,:))))-log(sum(w2.*exp(-DKLg(a,:)))));
        end
    case 5%The product of Gaussians approximation in [1] and Section 5 of
          %[3].

        P1Inv=zeros(d,d,n1);
        P2Inv=zeros(d,d,n2);
        for k=1:n1
            P1Inv(:,:,k)=inv(P1(:,:,k));
        end
        for k=1:n2
            P2Inv(:,:,k)=inv(P2(:,:,k));
        end

        %We have to first compute the z and t terms. The logarithms of the z
        %and t terms are given in Equation 22 in [1]. However, that
        %expression is incorrect; the corrected expression is used here

        zaap=zeros(n1,n1);
        tab=zeros(n1,n2);
        for a=1:n1
            for ap=a:n1
                diff=mu1(:,a)-mu1(:,ap);
                Sigmad=inv(P1Inv(:,:,a)+P1Inv(:,:,ap));
                const=sqrt(det(2*pi*Sigmad))/(sqrt(det(2*pi*P1(:,:,a))*det(2*pi*P1(:,:,ap))));

                B=P1Inv(:,:,a)*Sigmad*P1Inv(:,:,ap);
                zaap(a,ap)=const*exp(-(1/2)*diff'*B*diff);
                zaap(ap,a)=zaap(a,ap);
            end

            for b=1:n2
                diff=mu1(:,a)-mu2(:,b);
                Sigmad=inv(P1Inv(:,:,a)+P2Inv(:,:,b));
                const=sqrt(det(2*pi*Sigmad))/(sqrt(det(2*pi*P1(:,:,a))*det(2*pi*P2(:,:,b))));

                B=P1Inv(:,:,a)*Sigmad*P2Inv(:,:,b);
                tab(a,b)=const*exp(-(1/2)*diff'*B*diff);
            end
        end

        %Evaluate the sum in Equation 12 of [1], which is the same as
        %Equation 12 of [3].
        DVal=0;
        w1=w1(:).';
        w2=w2(:).';
        for a=1:n1
            DVal=DVal+w1(a)*(log(sum(w1.*zaap(a,:)))-log(sum(w2.*tab(a,:))));
        end
    case 6
        DProd=KLDistApproxGaussMix(w1,mu1,P1,w2,mu2,P2,5);
        DVar=KLDistApproxGaussMix(w1,mu1,P1,w2,mu2,P2,4);
        DVal=(DProd+DVar)/2;
    case 7%Goldberger's approximation from Equation 2 of [2].
        DKLg=zeros(n1,n2);
        for a=1:n1
            for b=1:n2
                DKLg(a,b)=KLDistGauss(mu1(:,a),P1(:,:,a),mu2(:,b),P2(:,:,b));
            end
        end

        w2=w2(:).';

        %Perform the minimization in Equation 1 of [2], which is also given
        %in Equation 14 of [3] and simultaneously evaluate Equation 2 of
        %[2], which is the same as Equation 15 in [3].
        DVal=0;
        logw2=log(w2);
        for a=1:n1
            minVal=min(DKLg(a,:)-logw2);
            DVal=DVal+w1(a)*(log(w1(a))+minVal);
        end
    otherwise
        error('Unknown algorithm specified.')
end
end

function val=KLArgGaussMix(x,w1,mu1,P1Inv,P1InvDet,w2,mu2,P2Inv,P2InvDet)
%%KLSAMPGAUSSMIX The KL divergence is an integral. This evaluates the
%                argument of the integral for the KL divergences between
%                two Gaussian mixtures.
%
%July 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.

PDF1=GaussianMixtureD.PDFI(x,w1,mu1,P1Inv,P1InvDet);
PDF2=GaussianMixtureD.PDFI(x,w2,mu2,P2Inv,P2InvDet);

val=PDF1.*(log(PDF1)-log(PDF2));
val(isnan(val))=0;

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 varargout = rotation3dToEulerAngles(mat, varargin)
%ROTATION3DTOEULERANGLES Extract Euler angles from a rotation matrix.
%
%   [PHI, THETA, PSI] = rotation3dToEulerAngles(MAT)
%   Computes Euler angles PHI, THETA and PSI (in degrees) from a 3D 4-by-4
%   or 3-by-3 rotation matrix.
%
%   ANGLES = rotation3dToEulerAngles(MAT)
%   Concatenates results in a single 1-by-3 row vector. This format is used
%   for representing some 3D shapes like ellipsoids.
%
%   ... = rotation3dToEulerAngles(MAT, CONVENTION)
%   CONVENTION specifies the axis rotation sequence. Default is 'ZYX'.
%   Supported conventions are:
%       'ZYX','ZXY','YXZ','YZX','XYZ','XZY'
%       'ZYZ','ZXZ','YZY','YXY','XZX','XYX'
%
%   Example
%   rotation3dToEulerAngles
%
%   References
%   Code from '1994 - Shoemake - Graphics Gems IV: Euler Angle Conversion:
%   http://webdocs.cs.ualberta.ca/~graphics/books/GraphicsGems/gemsiv/euler_angle/EulerAngles.c
%   (see rotm2eul, that is part of MATLAB's Robotics System Toolbox)
%   Modified using explanations in:
%   http://www.gregslabaugh.net/publications/euler.pdf
%   https://www.geometrictools.com/Documentation/EulerAngles.pdf
%
%   See also
%   transforms3d, rotation3dAxisAndAngle, createRotation3dLineAngle,
%   eulerAnglesToRotation3d
%

% ------
% Authors: David Legland, oqilipo
% E-mail: david.legland@grignon.inra.fr
% Created: 2010-08-11, using Matlab 7.9.0.529 (R2009b)
% Copyright 2010-2022 INRA - Cepia Software Platform

p = inputParser;
validStrings = {...
    'ZYX','ZXY','YXZ','YZX','XYZ','XZY',...
    'ZYZ','ZXZ','YZY','YXY','XZX','XYX'};
addOptional(p,'convention','ZYX',@(x) any(validatestring(x,validStrings)));
logParValidFunc = @(x) (islogical(x) || isequal(x,1) || isequal(x,0));
addParameter(p,'IsRotation', 1, logParValidFunc);
valTol = @(x) validateattributes(x,{'numeric'},{'scalar', '>=',eps(class(mat)), '<=',1});
addParameter(p,'tolerance', 1e-8, valTol);
parse(p,varargin{:});
convention=p.Results.convention;
isRotation = p.Results.IsRotation;
tolerance = p.Results.tolerance;

if isRotation
    if ~isTransform3d(mat(1:3,1:3), 'rotation', 1, 'tolerance', tolerance)
        warning(['Rotation matrix contains reflection or scaling ' ...
            'tested with a tolerance of ' num2str(tolerance) '.' newline ...
            'Calculation of euler angles might be incorrect.'])
    end
end

switch convention
    case 'ZYX'
        % extract |cos(theta)|
        cy = hypot(mat(1,1), mat(2,1));
        % avoid dividing by 0
        if cy > 16*eps
            % normal case: theta <> 0
            phi   = atan2( mat(2,1), mat(1,1));
            theta = atan2(-mat(3,1), cy);
            psi   = atan2( mat(3,2), mat(3,3));
        else
            phi   = 0;
            theta = atan2(-mat(3,1), cy);
            psi   = atan2(-mat(2,3), mat(2,2));
        end
    case 'ZXY'
        cy = hypot(mat(2,2), mat(1,2));
        if cy > 16*eps
            phi   = -atan2( mat(1,2), mat(2,2));
            theta = -atan2(-mat(3,2), cy);
            psi   = -atan2( mat(3,1), mat(3,3));
        else
            phi   = 0;
            theta = -atan2(-mat(3,2), cy);
            psi   = -atan2(-mat(1,3), mat(1,1));
        end
    case 'YXZ'
        cy = hypot(mat(3,3), mat(1,3));
        if cy > 16*eps
            phi   = atan2( mat(1,3), mat(3,3));
            theta = atan2(-mat(2,3), cy);
            psi   = atan2( mat(2,1), mat(2,2));
        else
            phi   = 0;
            theta = atan2(-mat(2,3), cy);
            psi   = atan2(-mat(1,2), mat(1,1));
        end
    case 'YZX'
        cy = hypot(mat(1,1), mat(3,1));
        if cy > 16*eps
            phi   = -atan2( mat(3,1), mat(1,1));
            theta = -atan2(-mat(2,1), cy);
            psi   = -atan2( mat(2,3), mat(2,2));
        else
            phi   = 0;
            theta = -atan2(-mat(2,1), cy);
            psi   = -atan2(-mat(3,2), mat(3,3));
        end
    case 'XYZ'
        cy = hypot(mat(3,3), mat(2,3));
        if cy > 16*eps
            phi   = -atan2( mat(2,3), mat(3,3));
            theta = -atan2(-mat(1,3), cy);
            psi   = -atan2( mat(1,2), mat(1,1));
        else
            phi   = 0;
            theta = -atan2(-mat(1,3), cy);
            psi   = -atan2(-mat(2,1), mat(2,2));
        end
    case 'XZY'
        cy = hypot(mat(2,2), mat(3,2));
        if cy > 16*eps
            phi   = atan2( mat(3,2), mat(2,2));
            theta = atan2(-mat(1,2), cy);
            psi   = atan2( mat(1,3), mat(1,1));
        else
            phi   = 0;
            theta = atan2(-mat(1,2), cy);
            psi   = atan2(-mat(3,1), mat(3,3));
        end

    case 'ZYZ'
        cy = hypot(mat(3,2), mat(3,1));
        if cy > 16*eps
            phi   = -atan2(mat(2,3), -mat(1,3));
            theta = -atan2(cy, mat(3,3));
            psi   = -atan2(mat(3,2), mat(3,1));
        else
            phi   = 0;
            theta = -atan2(cy, mat(3,3));
            psi   = -atan2(-mat(2,1), mat(2,2));
        end
    case 'ZXZ'
        cy = hypot(mat(3,2), mat(3,1));
        if cy > 16*eps
            phi   = atan2(mat(1,3), -mat(2,3));
            theta = atan2(cy, mat(3,3));
            psi   = atan2(mat(3,1), mat(3,2));
        else
            phi   = 0;
            theta = atan2(cy, mat(3,3));
            psi   = atan2(-mat(1,2), mat(1,1));
        end
    case 'YZY'
        cy = hypot(mat(2,3), mat(2,1));
        if cy > 16*eps
            phi   = atan2(mat(3,2), -mat(1,2));
            theta = atan2(cy, mat(2,2));
            psi   = atan2(mat(2,3), mat(2,1));
        else
            phi   = 0;
            theta = atan2(cy, mat(2,2));
            psi   = atan2(-mat(3,1), mat(3,3));
        end
    case 'YXY'
        cy = hypot(mat(2,3), mat(2,1));
        if cy > 16*eps
            phi   = -atan2(mat(1,2), -mat(3,2));
            theta = -atan2(cy, mat(2,2));
            psi   = -atan2(mat(2,1), mat(2,3));
        else
            phi   = 0;
            theta = -atan2(cy, mat(2,2));
            psi   = -atan2(-mat(1,3), mat(1,1));
        end
    case 'XZX'
        cy = hypot(mat(1,3), mat(1,2));
        if cy > 16*eps
            phi   = -atan2(mat(3,1), -mat(2,1));
            theta = -atan2(cy, mat(1,1));
            psi   = -atan2(mat(1,3), mat(1,2));
        else
            phi   = 0;
            theta = -atan2(cy, mat(1,1));
            psi   = -atan2(-mat(3,2), mat(3,3));
        end
    case 'XYX'
        cy = hypot(mat(1,2), mat(1,3));
        if cy > 16*eps
            phi   = atan2(mat(2,1), -mat(3,1));
            theta = atan2(cy, mat(1,1));
            psi   = atan2(mat(1,2), mat(1,3));
        else
            phi   = 0;
            theta = atan2(cy, mat(1,1));
            psi   = atan2(-mat(2,3), mat(2,2));
        end
end

% format output arguments
if nargout <= 1
    % one array
    varargout{1} = rad2deg([phi theta psi]);
else
    % three separate arrays
    varargout = cellfun(@rad2deg, {phi theta psi},'uni',0);
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function findelem(node,elem,range,varargin)
%% FINDELEM highlights some elements
%
%    FINDELEM(node,elem,range) finds all elements whose indices are in the
%    range array by displaying these elements in yellow. The size of the marker
%    is proportional to the number of digits of the indices. The font will
%    magnify as one resizes the figure.
%
%    FINDELEM(node,elem) finds all elements.
%
%    FINDELEM(node,elem,range,'noindex') skip the display of indices.
%
%    FINDELEM(node,elem,range,'param','value','param','value'...) allows
%    additional patch param/value pairs to highlight the elements.
%
% Example:
%     node = [1,0; 1,1; 0,1; -1,1; -1,0; -1,-1; 0,-1; 0,0];
%     elem = [1,2,8; 3,8,2; 8,3,5; 4,5,3; 7,8,6; 5,6,8];
%     subplot(1,2,1);
%     showmesh(node,elem);
%     findelem(node,elem,1,'index','FaceColor','r');
%     subplot(1,2,2);
%     showmesh(node,elem);
%     findelem(node,elem);
%
%   See also findelem3, findnode3, findedge.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

hold on

if (nargin==2) || isempty(range) || (ischar(range) && strcmp(range,'all'))
    range = (1:size(elem,1))';
end
if islogical(range)
    range = find(range);
end
if size(range,2)>size(range,1)
    range = range';
end
nV = size(elem,2);
if nV == 3
    center=(node(elem(range,1),:)+node(elem(range,2),:)+node(elem(range,3),:))/3;
elseif nV == 4
    center=(node(elem(range,1),:)+node(elem(range,2),:)+node(elem(range,3),:)+node(elem(range,4),:))/4;
end
if length(range) < size(elem,1)
    x = reshape(node(elem(range,:),1),size(range,1), size(elem,2))';
    y = reshape(node(elem(range,:),2),size(range,1), size(elem,2))';
    h = patch(x,y,'y');
    if nargin > 3
        if strcmp(varargin{1},'noindex')
           if size(varargin,2)>=2
                set(h,varargin{2:end});
           end
        else
            set(h,varargin{:});
        end
    end
end
nDigit =ceil(log10(range+1));
oSize = 400*nDigit;
f = gcf;
set(f, 'Units', 'pixel');
fHeight = f.Position(4);
if (nargin <=3) || ~(strcmp(varargin{1},'noindex'))
    if size(node,2) == 2
        s = scatter(center(:,1),center(:,2),oSize,'o','LineWidth',1,'MarkerEdgeColor','k',...
            'MarkerFaceColor','y');
        if max(f.Position(3:4)) <= 1
            fSz = 20*f.Position(4); % font size is proportional to height
        else
            fSz = log(f.Position(3)*f.Position(4));
        end
        tShift = 0.02*(max(nDigit) + 1 - nDigit);
        tShift(nDigit==1) = 0.5*tShift(nDigit==1);
        t = text(center(:,1)-tShift,center(:,2),int2str(range),...
            'Fontsize',fSz,'FontWeight','bold','Color','k');
        nElemShown = length(t);
        set(f,'SizeChangedFcn',@resizeCallback);
        set(f,'Visible','on');
        drawnow;
    elseif size(node,2) == 3 % surface mesh
        scatter3(center(:,1),center(:,2),center(:,3),oSize,'o','LineWidth',1,...
            'MarkerEdgeColor','k','MarkerFaceColor','y');
        text(center(:,1)-0.02*nDigit,center(:,2),center(:,3),int2str(range),'FontSize',12,...
        'FontWeight','bold','Color','k');
    end
end
hold off
%% relative position of numbers
% add by Shuhao Cao
    function resizeCallback(f, ~)
        if max(f.Position(3:4)) <= 1 % relative
            newFSz = 80*f.Position(4); % font size is proportional to height
            newOSize = 10*oSize*(f.Position(4));
            newTshift = 4*tShift*f.Position(4);
        else % pixel
            newFSz = log(f.Position(3)*f.Position(4));
            newOSize = oSize*(f.Position(4)/fHeight);
            newTshift = tShift/(f.Position(4)/fHeight);
        end
        % change font size accordingly
        set(t,'FontSize',newFSz);
        set(s,'SizeData',newOSize);
        for i = 1:nElemShown
            set(t(i), 'Position', [center(i,1)-newTshift(i), center(i,2), 0]);
        end
    end

end


text:
%PRGLOBAL Define or reset PRTools globals
%
%   PRGLOBAL
%   PRGLOBAL RESET
%   PRGLOBAL(PAR1,VAL1,PAR2,VAL2, ....)
%
% INPUT
%   PAR1,PAR2  String with variable name of global variable
%   VAL1,VAL2  Desired value of global variable
%
% DESCRIPTION
% The functioning of PRTools depends on the values of a number of global
% variables. Many of them have their own routine to change them. This
% routine is a central routine to list them (no arguments), reset them
% to their initial values or set them to some user preference.
%
% The following global variables can be changed:
% PRMEMORY,  see corresponding routine
% DEFAULTBATCHSIZE,  see SETBATCH
% GRIDSIZE,  see corresponding routine and PLOTC
% CHECK_SIZES,  TRUE or FALSE, disables or enables some size checking
% STAMP_MAP,  TRUE or FALSE, see corresponding routine
% PRWAITBAR,  ON or OFF, see corresponding routine
% PRWARNING,  set level, see corresponding routine
% REGOPT_NFOLDS,  see REGOPTC
% REGOPT_ITERMAX,  see REGOPTC
% REGOPT_REPS,  see REGOPTC
%
% For more information just type PRGLOBAL
%
% SEE ALSO (PRTools Guide)
% PRMEMORY, SETBATCH, GRIDSIZE, STAMP_MAP, PRWAITBAR, PRWARNING, REGOPTC

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com

function prglobal(varargin)

global CHECK_SIZES
global REGOPT_NFOLDS REGOPT_REPS REGOPT_ITERMAX REGOPT_OPTCRIT
global DEFAULTBATCHSIZE
if isempty(DEFAULTBATCHSIZE)
  DEFAULTBATCHSIZE = 1000;
end

if nargin == 0
  s = '\n          ';
  fprintf(['\n%5i  PRMEMORY/1e6, the maximum size of matrices intermediately created by' ...
    s 'PRTools. If set too high Matlab may crash. If set low PRTools may' s ...
    'generate many loops. \n'],prmemory/1e6);
  fprintf(['\n%5i  DEFAULTBATCHSIZE, the batch size set as default by setbatch. If set ' ...
    s 'too high Matlab may crash. If set low PRTools may generate many loops. \n'], ...
    DEFAULTBATCHSIZE);
  fprintf(['\n%5i  GRIDSIZE, to be used in scatterplots. Its value determines the' ...
    s 'accuracy by which classifiers are drawn using plotc. The default is' s ...
    '30 but sometimes 300 or more is needed causing long computing times.\n'],gridsize);
  fprintf(['\n%5i  CHECK_SIZES, true (1) or false(0), determines whether in sequential' ...
    s 'combining of datassets and mappings output and input sizes should be' s ...
    'checked on matching.\n'],CHECK_SIZES);
  fprintf(['\n%5i  STAMP_MAP, disabled (0), retrieval only (1), storage and retrieval (2)' ...
    s 'of trained mappings to avoid recomputation of precomputed results.\n'],stamp_map);
  fprintf(['\n%5s  PRWAITBAR, the status of prwaitbar, ''on'' or ''off''\n'],prwaitbar);
  fprintf(['\n%5i  PRWARNING, the level of prwarning.\n'],prwarning);
  fprintf(['\n%5i  REGOPT_NFOLDS, the number of folds used in the crossvalidation' ...
    s 'procedure by which the classifier performance is estimated for a' s ...
    'single setting of the parameters during their optimization by REGOPTC.\n'],REGOPT_NFOLDS);
  fprintf(['\n%5i  REGOPT_ITERMAX, the maximum number of iterations in the optimization' ...
    s 'procedure by which the classifier performance is estimated for a' s ...
    'single setting of the parameters during their optimization by REGOPTC.\n'],REGOPT_ITERMAX);
  fprintf(['\n%5i  REGOPT_REPS, the number of repetitions used in the crossvalidation' ...
    s 'procedure by which the classifier performance is estimated for a' s ...
    'single setting of the parameters during their optimization by REGOPTC.\n'],REGOPT_REPS);
  fprintf(['\n          Note that automatic parameter optimization for classifiers may slow' ...
    s 'down training by a factor REGOPT_NFOLDS*REGOPT_ITERMAX*REGOPT_REPS.\n']);

elseif (nargin == 1) & strcmpi(varargin{1},'reset')
  prmemory(50000000);
  gridsize(30);
  prwarning(1);
  prwaitbar;
  CHECK_SIZES = true;
  stamp_map(0);
  REGOPT_NFOLDS = 5;
  REGOPT_REPS = 1;
  REGOPT_ITERMAX = 20;
  DEFAULTBATCHSIZE = 1000;
elseif nargin ~= 2*floor(nargin/2)
  error('Wrong number of input arguments')
else
  for j=1:2:nargin
    switch lower(varargin{j})
      case 'globalprmemory'
        prmemory(varargin{j+1});
      case 'defaultbatchsize'
        DEFAULTBATCHSIZE = varargin{j+1};
      case 'current_gridsize'
        gridsize(varargin{j+1});
      case 'prwaitbar'
        prwaitbar(varargin{j+1});
      case 'prwarning'
        prwarning(varargin{j+1});
      case 'check_sizes'
        if all(varargin{j+1} ~= [true,false])
          error('CHECK_SIZES should either be ''true'' or ''false''')
        end
        CHECK_SIZES = varargin{j+1};
      case 'stamp_map'
        stamp_map(varargin{j+1});
      case 'regopt_nfolds'
        REGOPT_NFOLDS = varargin{j+1};
      case 'regopt_reps'
        REGOPT_REPS = varargin{j+1};
      case 'regopt_itermax'
        REGOPT_ITERMAX = varargin{j+1};
      otherwise
        error('Unknown global parameter')
    end
  end
end

author: marianux

repo: ecg-kit

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


text:
function [nodes2, edges2] = grMergeNodeClusters(nodes, edges)
%GRMERGENODECLUSTERS Merge cluster of connected nodes in a graph
%
%   grMergeNodeClusters(nodes, edges)
%   Detects groups of nodes that belongs to the same global node, and
%   replace them by a unique node. Coordinates of reference node is given
%   by the median coordinates of cluster nodes.
%
%   This function is intended to be used as filter after a binary image
%   skeletonization and vectorization.
%
%
%   See Also
%   grMergeNodesMedian
%

%   -----
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 13/08/2003.
%

%   HISTORY

%% Initialization

% intialize result
nodes2 = nodes;
edges2 = edges;

% compute degree of each node
degrees = grNodeDegree(1:size(nodes, 1), edges)';

% find index of multiple nodes
indMul = find(degrees > 2);

% indices of edges that link several multiple nodes
indEdges = sum(ismember(edges, indMul), 2) == 2;

% associate a label to each cluster
labels = grLabel(nodes, edges(indEdges, :));
clusterLabels = unique(labels(indMul));

%% Replace each cluster by median point

% iterate on clusters
for i = 1:length(clusterLabels)
    % indices of nodes of the current cluster
    inds = find(labels == clusterLabels(i));

    % coordinates of new reference node
    clusterNodes = nodes(inds, :);
    medianNode = median(clusterNodes, 1);

    % replace coordinates of reference node
    refNode = min(inds);
    nodes2(refNode, :) = medianNode;

    % replace node indices in edge array
    edges2(ismember(edges2, inds)) = refNode;
end

%% Clean up

% keep only relevant nodes
inds = unique(edges2(:));
nodes2 = nodes2(inds, :);

% relabeling of edges
for i = 1:length(inds)
    edges2(edges2 == inds(i)) = i;
end

% remove double edges
edges2 = unique(sort(edges2, 2), 'rows');

% remove 'loops'
edges2(edges2(:,1) == edges2(:,2), :) = [];

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function Offspring = Exploration(Problem,PC,NPC,nND,N)
% Individual exploration in BCE

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

    PCObj  = PC.objs;
    NPCObj = NPC.objs;

    %% Normalization
    fmax   = max(PCObj,[],1);
    fmin   = min(PCObj,[],1);
    PCObj  = (PCObj-repmat(fmin,size(PCObj,1),1))./repmat(fmax-fmin,size(PCObj,1),1);
    NPCObj = (NPCObj-repmat(fmin,size(NPCObj,1),1))./repmat(fmax-fmin,size(NPCObj,1),1);

    %% Determine the size of the niche
    d  = pdist2(PCObj,PCObj);
    d(logical(eye(length(d)))) = inf;
    d  = sort(d,2);
    r0 = mean(d(:,min(3,size(d,2))));
    r  = nND/N*r0;

    %% Detect the solutions in PC to be explored
    d = pdist2(PCObj,NPCObj);
    S = find(sum(d<=r,2)<=1);

    %% Generate new solutions
    if ~isempty(S)
        MatingPool = randi(length(PC),1,length(S));
        Offspring  = OperatorGAhalf(Problem,PC([S',MatingPool]));
    else
        Offspring = [];
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function ft_plot_cloud(pos, val, varargin)

% FT_PLOT_CLOUD visualizes spatially sparse scalar data as spheres or
% spherical clouds of points and optionally 2D slices through those clouds
%
% Use as
%   ft_plot_cloud(pos, val, ...)
% where the first argument are the positions and the second argument are the values
% for each location.
%
% Optional input arguments should come in key-value pairs and can include
%   'cloudtype'          = 'cloud' (default) plots a group of spherically arranged points at each sensor position
%                          'surf' plots a single spherical surface mesh at each sensor position
%   'scalerad'           = scale radius with val, can be 'yes' or 'no' (default = 'yes')
%   'radius'             = scalar, maximum radius of cloud (default = 4 mm)
%   'clim'               = 1x2 vector specifying the min and max for the colorscale
%   'unit'               = string, convert the sensor array to the specified geometrical units (default = [])
%   'mesh'               = string or Nx1 cell-array, triangulated mesh(es), see FT_PREPARE_MESH
%   'slice'              = requires 'mesh' as input (default = 'none')
%                          '2d', plots 2D slices through the cloud with an outline of the mesh
%                          '3d', draws an outline around the mesh at a particular slice
%
% The following inputs apply when 'cloudtype' = 'cloud'
%   'rmin'               = scalar >= 1, minimum radius of cloud if scalerad = 'yes' (default = 1 mm)
%   'colormap'           = colormap for functional data, see COLORMAP
%   'colorgrad'          = 'white' or a scalar (e.g. 1), degree to which the saturatoin of points in cloud changes from its center
%   'ptsize'             = scalar, size of points in cloud (default = 1 mm)
%   'ptdensity'          = scalar, density of points in cloud (default = 20 per mm^3)
%   'ptgradient'         = scalar, degree to which density of points in cloud changes from its center (default = 0.5, i.e. uniform density)
%
% The following inputs apply when 'slice' = '2d' or '3d'
%   'ori'                = 'x', 'y', or 'z', specifies the orthogonal plane which will be plotted (default = 'y')
%   'slicepos'           = 'auto' or Nx1 vector specifying the position of the
%                          slice plane along the orientation axis (default = 'auto': chooses slice(s) with
%                          the most data)
%   'nslices'            = scalar, number of slices to plot if 'slicepos' = 'auto (default = 1)
%   'minspace'           = scalar, minimum spacing between slices if nslices>1
%                          (default = 1)
%   'intersectcolor'     = string, Nx1 cell-array, or Nx3 vector specifying line color (default = 'k')
%   'intersectlinestyle' = string or Nx1 cell-array, line style specification (default = '-')
%   'intersectlinewidth' = scalar or Nx1 vector, line width specification (default = 2)
%
% The following inputs apply when 'cloudtype' = 'surf' and 'slice' = '2d'
%   'ncirc'           = scalar, number of concentric circles to plot for each
%                       cloud slice (default = 15) make this hidden or scale
%   'scalealpha'      = 'yes' or 'no', scale the maximum alpha value of the center circle
%                       with distance from center of cloud
%
% See also FT_ELECTRODEPLACEMENT, FT_PLOT_SENS, FT_PLOT_TOPO, FT_PLOT_TOPO3D

% Copyright (C) 2017-2018, Arjen Stolk, Sandon Griffin
%
% 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$

% run some checks
if size(pos,2)~=3
  ft_error('positions shoudl be specified as an Nx3 array')
end

if isempty(val)
  val = ones(size(pos,1),1); % vector of ones
end

assert(isrow(val) || iscolumn(val), 'values should be represented as a vector')
val = val(:); % ensure it is a column

unit = ft_getopt(varargin, 'unit');
if isempty(unit)
  % estimate the unit of geometry of the positions, this is used for the other defaults
  unit = ft_estimate_units(norm(range(pos)));
end

% get the generic input arguments
cloudtype          = ft_getopt(varargin, 'cloudtype', 'cloud');
radius             = ft_getopt(varargin, 'radius', 4 * ft_scalingfactor('mm', unit));
rmin               = ft_getopt(varargin, 'rmin',   1 * ft_scalingfactor('mm', unit));
scalerad           = ft_getopt(varargin, 'scalerad', 'yes');
colorgrad          = ft_getopt(varargin, 'colorgrad', 'white');
clim               = ft_getopt(varargin, 'clim');
meshplot           = ft_getopt(varargin, 'mesh');

if ft_platform_supports('parula')
  cmap             = ft_getopt(varargin, 'colormap', 'parula');
else
  cmap             = ft_getopt(varargin, 'colormap', 'jet');
end

% cloud related inputs
ptsize             = ft_getopt(varargin, 'ptsize', 1 * ft_scalingfactor('mm', unit));
ptdensity          = ft_getopt(varargin, 'ptdensityity', 20 / ft_scalingfactor('mm', unit)^3);  % points per unit of volume
ptgradient         = ft_getopt(varargin, 'ptgradient', .5);

% slice related inputs
sli                = ft_getopt(varargin, 'slice', 'none');
ori                = ft_getopt(varargin, 'ori', 'y');
slicepos           = ft_getopt(varargin, 'slicepos', 'auto');
nslices            = ft_getopt(varargin, 'nslices', 1);
minspace           = ft_getopt(varargin, 'minspace', 1);
intersectcolor     = ft_getopt(varargin, 'intersectcolor', {'k'});
intersectlinestyle = ft_getopt(varargin, 'intersectlinestyle', {'-'});
intersectlinewidth = ft_getopt(varargin, 'intersectlinewidth', 2);
ncirc              = ft_getopt(varargin, 'ncirc', 15);
scalealpha         = ft_getopt(varargin, 'scalealpha', 'no');

% mesh related inputs
facecolor          = ft_getopt(varargin, 'facecolor', [0.781 0.762 0.664]);
edgecolor          = ft_getopt(varargin, 'edgecolor', 'none');
facealpha          = ft_getopt(varargin, 'facealpha', 1);
edgealpha          = ft_getopt(varargin, 'edgealpha', 0);
vertexcolor        = ft_getopt(varargin, 'vertexcolor', 'curv');

% give some warnings for possibly inappropriate input
if radius < 1 * ft_scalingfactor('mm', unit)
  ft_warning('radius is smaller than 1 mm');
end
if rmin < 1 * ft_scalingfactor('mm', unit)
  ft_warning('rmin is smaller than 1 mm');
end

if ~isempty(meshplot)
  % mesh should be a cell-array
  if isstruct(meshplot)
    tmp = meshplot;
    meshplot = cell(size(tmp));
    for i=1:numel(tmp)
      meshplot{i} = tmp(i);
    end
  elseif iscell(meshplot)
    % do nothing
  else
    meshplot = {};
  end

  % replace pnt by pos
  for k = 1:numel(meshplot)
    meshplot{k} = fixpos(meshplot{k});
  end

  for k = 1:numel(meshplot)
    if ~isfield(meshplot{k}, 'pos') || ~isfield(meshplot{k}, 'tri')
      % ft_error('the mesh should be a structure with pos and tri');
      meshplot{k}.pos = [];
      meshplot{k}.tri = [];
    end
  end

  % facecolor, edgecolor, and vertexcolor should be cell-array
  if ~iscell(facecolor)
    tmp = facecolor;
    if ischar(tmp)
      facecolor = {tmp};
    elseif ismatrix(tmp) && size(tmp, 2) == 3
      facecolor = cell(size(tmp,1), 1);
      for i=1:size(tmp,1)
        facecolor{i} = tmp(i, 1:3);
      end
    else
      facecolor = {};
    end
  end
  if ~iscell(edgecolor)
    tmp = edgecolor;
    if ischar(tmp)
      edgecolor = {tmp};
    elseif ismatrix(tmp) && size(tmp, 2) == 3
      edgecolor = cell(size(tmp,1), 1);
      for i=1:size(tmp,1)
        edgecolor{i} = tmp(i, 1:3);
      end
    else
      edgecolor = {};
    end
  end
  if ~iscell(vertexcolor)
    tmp = vertexcolor;
    if ischar(tmp)
      vertexcolor = {tmp};
    elseif ismatrix(tmp) && size(tmp, 2) == 3
      vertexcolor = cell(size(tmp,1), 1);
      for i=1:size(tmp,1)
        vertexcolor{i} = tmp(i, 1:3);
      end
    else
      vertexcolor = {};
    end
  end

  % make sure each mesh has plotting options specified
  if numel(meshplot) > 1
    nmesh = numel(meshplot);
    if numel(facecolor) < numel(meshplot)
      for m = numel(facecolor)+1:nmesh
        facecolor{m} = facecolor{1};
      end
    end
    if numel(edgecolor) < numel(meshplot)
      for m = numel(edgecolor)+1:nmesh
        edgecolor{m} = edgecolor{1};
      end
    end
    if numel(facealpha) < numel(meshplot)
      for m = numel(facealpha)+1:nmesh
        facealpha(m) = facealpha(1);
      end
    end
    if numel(edgealpha) < numel(meshplot)
      for m = numel(edgealpha)+1:nmesh
        edgealpha(m) = edgealpha(1);
      end
    end
    if numel(vertexcolor) < numel(meshplot)
      for m = numel(vertexcolor)+1:nmesh
        vertexcolor(m) = vertexcolor(1);
      end
    end
  end
end

if strcmp(sli, '2d') || strcmp(sli, '3d')
  if isempty(meshplot)
    ft_error('plotting a slice requires a mesh as input')
  else
    dointersect = 1;
  end
else
  dointersect = 0;
end

if dointersect % check intersection inputs
  % color and linestyle should be cell-array
  if ~iscell(intersectcolor)
    tmp = intersectcolor;
    if ischar(tmp)
      intersectcolor = {tmp};
    elseif ismatrix(tmp) && size(tmp, 2) == 3
      intersectcolor = cell(size(tmp,1), 1);
      for i=1:size(tmp,1)
        intersectcolor{i} = tmp(i, 1:3);
      end
    else
      intersectcolor = {};
    end
  end
  if ischar(intersectlinestyle)
    intersectlinestyle = {intersectlinestyle};
  elseif iscell(intersectlinestyle)
    % do nothing
  else
    intersectlinestyle = {};
  end

  % make sure each intersection has plotting options specified
  if numel(meshplot) > 1
    nmesh = numel(meshplot);
    if numel(intersectcolor) < numel(meshplot)
      for m = numel(intersectcolor)+1:nmesh
        intersectcolor{m} = intersectcolor{1};
      end
    end
    if numel(intersectlinewidth) < numel(meshplot)
      for m = numel(intersectlinewidth)+1:nmesh
        intersectlinewidth(m) = intersectlinewidth(1);
      end
    end
    if numel(intersectlinestyle) < numel(meshplot)
      for m = numel(intersectlinestyle)+1:nmesh
        intersectlinestyle{m} = intersectlinestyle{1};
      end
    end
  end % end intersection plotting checks
end % end dointersect checks

if dointersect
  % set the orientation of the slice plane
  if strcmp(ori, 'x')
    oriX = 1; oriY = 0; oriZ = 0;
  elseif strcmp(ori, 'y')
    oriX = 0; oriY = 1; oriZ = 0;
  elseif strcmp(ori, 'z')
    oriX = 0; oriY = 0; oriZ = 1;
  else
    ft_error('ori must be "x", "y" or "z"')
  end
end

if isempty(clim)
  clim = [min(val) max(val)]; % use the data
end

% functional data scaling factors
if ischar(cmap)
  if strcmp(cmap, 'default')
    cmapsc = get(0, 'DefaultFigureColormap');
  else
    cmapsc = feval(cmap, 201); % an odd number
  end
else
  cmapsc = cmap;
end

cmid   = size(cmapsc,1)/2;                          % colorbar middle
colscf = 2*( (val-clim(1)) / (clim(2)-clim(1)) )-1; % color between -1 and 1, bottom vs. top colorbar
colscf(colscf>1)=1; colscf(colscf<-1)=-1;           % clip values outside the [-1 1] range

radscf = abs(val-min(abs(val)) * sign(max(val)));   % 'vdown' representation
radscf = 2*( radscf / (clim(2)-clim(1)) );          % radius between 0 and 1, small vs. large pos/neg effect
radscf(radscf>1)=1; radscf(radscf<0)=0;             % clip values outside the [0 1] range

if strcmp(scalerad, 'yes')
  rmax = rmin+(radius-rmin)*radscf; % maximum radius of the clouds
else
  rmax = ones(length(pos), 1)*radius; % each cloud has the same radius
end

if dointersect
  % generate circle points
  angles = linspace(0,2*pi,50);
  x = cos(angles)';
  y = sin(angles)';
  slicedim = zeros(length(angles),1);

  if strcmp(slicepos, 'auto') % search each slice for largest area of data
    % find the potential limits of the interpolation
    intxmax = max(pos(:,1))+radius; intxmin = min(pos(:,1))-radius;
    intymax = max(pos(:,2))+radius; intymin = min(pos(:,2))-radius;
    intzmax = max(pos(:,3))+radius; intzmin = min(pos(:,3))-radius;

    % define potential slices with data
    if oriX; potent_slices = round(intxmin):round(intxmax); end
    if oriY; potent_slices = round(intymin):round(intymax); end
    if oriZ; potent_slices = round(intzmin):round(intzmax); end

    area = NaN(length(pos),length(potent_slices));  % preallocate matrix of electrode interpolation area for each slice
    for s = 1:length(potent_slices)                 % only search slices that could potentially contain data
      distance = NaN(length(pos),1);                % preallocate vector for each electrodes distance from the slice
      for c = 1:length(pos)
        indpos = pos(c, :);
        if oriX; distance(c) = abs(indpos(1)-potent_slices(s)); end
        if oriY; distance(c) = abs(indpos(2)-potent_slices(s)); end
        if oriZ; distance(c) = abs(indpos(3)-potent_slices(s)); end

        if distance(c) < rmax(c) % if there is any data from this electrode in this slice
          % find the circle points for the interpolation of this electrode
          xmax = rmax(c)*x;
          ymax = rmax(c)*y;

          % find the maximum radius of a cross section of the virtual sphere in the given slice
          xmaxdif = abs(xmax-distance(c));
          imindif = find(xmaxdif == min(xmaxdif), 1); % index of x value closest to distance(e)
          rcmax = abs(ymax(imindif));

          area(c, s) = 0.5*pi*rcmax^2;
        else % area must be zero
          area(c, s) = 0;
        end
      end % for each loop
    end % for each slice
    totalarea = sum(area);

    slicepos = zeros(nslices,1);
    for n = 1:nslices
      imaxslice = find(totalarea == max(totalarea), 1);       % index of the slice with the maximum area
      slicepos(n) = potent_slices(imaxslice);                 % position of the yet unlisted slice with the maximum area
      totalarea(imaxslice-minspace:imaxslice+minspace) = 0;   % change the totalarea of the chosen slice and those within minspace to 0 so that it is not chosen again
    end
  end

  % pre-allocate logical array specifying whether intersection with a given mesh (k) actually exists within a given slice (s)
  intersect_exists = zeros(numel(slicepos), numel(meshplot));
end

% draw figure
if strcmp(sli, '2d')
  % pre-allocate interpolation limits of each slice to facilitate inding overall limits of all slices after plotting
  xsmax = NaN(numel(slicepos),1); xsmin = NaN(numel(slicepos),1);
  ysmax = NaN(numel(slicepos),1); ysmin = NaN(numel(slicepos),1);
  zsmax = NaN(numel(slicepos),1); zsmin = NaN(numel(slicepos),1);

  for s = 1:numel(slicepos) % slice loop
    subplot(numel(slicepos),1,s); hold on;

    % pre-allocate interpolation limits of each cloud to facilitate finding slice limits after plotting
    xcmax = NaN(length(pos(:,1)),1); xcmin = NaN(length(pos(:,1)),1);
    ycmax = NaN(length(pos(:,1)),1); ycmin = NaN(length(pos(:,1)),1);
    zcmax = NaN(length(pos(:,1)),1); zcmin = NaN(length(pos(:,1)),1);

    for c = 1:length(pos(:,1)) % cloud loop
      indpos = pos(c, :);
      % calculate distance from slice
      if oriX; distance = abs(indpos(1)-slicepos(s)); end
      if oriY; distance = abs(indpos(2)-slicepos(s)); end
      if oriZ; distance = abs(indpos(3)-slicepos(s)); end

      if distance < rmax(c)
        switch (cloudtype)
          case 'surf'
            xmax = rmax(c)*x;
            ymax = rmax(c)*y;

            if strcmp(scalealpha, 'yes')
              maxalpha = (rmax(c)-distance)/rmax(c);
            else
              maxalpha = 1;
            end

            % find the maximum radius of a cross section of the virtual sphere in the given slice
            xmaxdif = abs(xmax-distance);
            imindif = find(xmaxdif == min(xmaxdif), 1); % index of x value closest to distance(e)
            rcmax = abs(ymax(imindif));

            % determine points along outermost circle
            xe = rcmax*x;
            ye = rcmax*y;

            % jitter values of points in the slice plane so no surfaces overlap
            slicedime = slicedim+(0.01*rand*ones(length(x), 1));

            % plot concentric circles
            for n = 0:ncirc-1 % circle loop
              xo = xe*((ncirc-n)/ncirc);    % outer x points
              yo = ye*((ncirc-n)/ncirc);    % outer z points
              xi = xe*((ncirc-1-n)/ncirc);  % inner x points
              yi = ye*((ncirc-1-n)/ncirc);  % inner z points
              if n == ncirc-1 % for the last concentric circle
                if oriX; hs = fill3(slicepos(s)+slicedime, indpos(2)+xo, indpos(3)+yo, val(c)); end
                if oriY; hs = fill3(indpos(1)+xo, slicepos(s)+slicedime, indpos(3)+yo, val(c)); end
                if oriZ; hs = fill3(indpos(1)+xo, indpos(2)+yo, slicepos(s)+slicedime, val(c)); end
              else
                if oriX; hs = fill3([slicepos(s)+slicedime; slicepos(s)+slicedim], [indpos(2)+xo; indpos(2)+xi], [indpos(3)+yo; indpos(3)+yi], val(c)); end
                if oriY; hs = fill3([indpos(1)+xo; indpos(1)+xi], [slicepos(s)+slicedime; slicepos(s)+slicedim], [indpos(3)+yo; indpos(3)+yi], val(c)); end
                if oriZ; hs = fill3([indpos(1)+xo; indpos(1)+xi], [indpos(2)+yo; indpos(2)+yi], [slicepos(s)+slicedime; slicepos(s)+slicedim], val(c)); end
              end
              set(hs, 'EdgeColor', 'none', 'FaceAlpha', maxalpha*n/ncirc)
            end % end circle loop

            % find the limits of the plotted surfaces for this electrode
            if oriX
              xcmax(c) = max(slicedime+slicepos(s)); xcmin(c) = min(slicedime+slicepos(s));
              ycmax(c) = max(xe+indpos(2)); ycmin(c) = min(xe+indpos(2));
              zcmax(c) = max(ye+indpos(3)); zcmin(c) = min(ye+indpos(3));
            elseif oriY
              xcmax(c) = max(xe+indpos(1)); xcmin(c) = min(xe+indpos(1));
              ycmax(c) = max(slicedime+slicepos(s)); ycmin(c) = min(slicedime+slicepos(s));
              zcmax(c) = max(ye+indpos(3)); zcmin(c) = min(ye+indpos(3));
            elseif oriZ
              xcmax(c) = max(xe+indpos(1)); xcmin(c) = min(xe+indpos(1));
              ycmax(c) = max(ye+indpos(2)); ycmin(c) = min(ye+indpos(2));
              zcmax(c) = max(slicedime+slicepos(s)); zcmin(c) = min(slicedime+indpos(s));
            end

          case 'cloud'
            rng(0, 'twister');                          % random number generator
            npoints = round(ptdensity*pi*rmax(c)^2);       % number of points based on area of cloud cross section
            azimuth = 2*pi*rand(npoints,1);             % azimuthal angle for each point
            radii = rmax(c)*(rand(npoints,1).^ptgradient);  % radius value for each point
            radii = sort(radii);                        % sort radii in ascending order so they are plotted from inside out
            % convert to Carthesian; note that second input controls third output
            if oriX; [y,z,x] = sph2cart(azimuth, zeros(npoints,1)+0.01*rand(npoints,1), radii); end
            if oriY; [x,z,y] = sph2cart(azimuth, zeros(npoints,1)+0.01*rand(npoints,1), radii); end
            if oriZ; [x,y,z] = sph2cart(azimuth, zeros(npoints,1)+0.01*rand(npoints,1), radii); end

            % color axis with radius scaling
            if strcmp(colorgrad, 'white') % color runs up to white
              fcolidx = ceil(cmid) + sign(colscf(c))*floor(abs(colscf(c)*cmid));
              if fcolidx == 0; fcolidx = 1; end
              fcol = cmapsc(fcolidx,:); % color [Nx3]
              ptcol = [linspace(fcol(1), 1, npoints)' linspace(fcol(2), 1, npoints)' linspace(fcol(3), 1, npoints)'];
            elseif isscalar(colorgrad) % color runs down towards colorbar middle
              rnorm = radii/rmax(c); % normalized radius
              if radscf(c)>=.5 % extreme values
                ptcol = val(c) - (flip(1-rnorm).^inv(colorgrad))*val(c); % scaled fun [Nx1]
              elseif radscf(c)<.5 % values closest to zero
                ptcol = val(c) + (flip(1-rnorm).^inv(colorgrad))*abs(val(c)); % scaled fun [Nx1]
              end
            else
              ft_error('color gradient should be either ''white'' or a scalar determining the fall-off')
            end

            % draw the points
            if oriX; scatter3(x+slicepos(s), y+indpos(2), z+indpos(3), ptsize, ptcol, '.'); end
            if oriY; scatter3(x+indpos(1), y+slicepos(s), z+indpos(3), ptsize, ptcol, '.'); end
            if oriZ; scatter3(x+indpos(1), y+indpos(2), z+slicepos(s), ptsize, ptcol, '.'); end

            % find the limits of the plotted points for this electrode
            if oriX
              xcmax(c) = max(x+slicepos(s)); xcmin(c) = min(x+slicepos(s));
              ycmax(c) = max(y+indpos(2)); ycmin(c) = min(y+indpos(2));
              zcmax(c) = max(z+indpos(3)); zcmin(c) = min(z+indpos(3));
            elseif oriY
              xcmax(c) = max(x+indpos(1)); xcmin(c) = min(x+indpos(1));
              ycmax(c) = max(y+slicepos(s)); ycmin(c) = min(y+slicepos(s));
              zcmax(c) = max(z+indpos(3)); zcmin(c) = min(z+indpos(3));
            elseif oriZ
              xcmax(c) = max(x+indpos(1)); xcmin(c) = min(x+indpos(1));
              ycmax(c) = max(y+indpos(2)); ycmin(c) = min(y+indpos(2));
              zcmax(c) = max(z+slicepos(s)); zcmin(c) = min(z+slicepos(s));
            end

          otherwise
            ft_error('unsupported cloudtype "%s"', cloudtype);

        end % switch cloudtype
      end % if distance < rmax(c)
    end % for each position

    if dointersect
      if oriX; ori = [1 0 0]; loc = [slicepos(s) 0 0]; end
      if oriY; ori = [0 1 0]; loc = [0 slicepos(s) 0]; end
      if oriZ; ori = [0 0 1]; loc = [0 0 slicepos(s)]; end

      % normalise the orientation vector to one
      ori = ori./sqrt(sum(ori.^2));

      % shift the location to be along the orientation vector
      loc = ori*dot(loc,ori);

      % determine three points on the plane
      inplane = eye(3) - (eye(3) * ori') * ori;
      v1 = loc + inplane(1,:);
      v2 = loc + inplane(2,:);
      v3 = loc + inplane(3,:);

      for k = 1:numel(meshplot)

        % only plot if the mesh actually intersects the plane
        xmmax = max(meshplot{k}.pos(:,1)); xmmin = min(meshplot{k}.pos(:,1));
        ymmax = max(meshplot{k}.pos(:,2)); ymmin = min(meshplot{k}.pos(:,2));
        zmmax = max(meshplot{k}.pos(:,3)); zmmin = min(meshplot{k}.pos(:,3));

        if oriX
          if slicepos(s) < xmmax && slicepos(s) > xmmin
            intersect_exists(s,k) = 1;
          else
            intersect_exists(s,k) = 0;
          end
        elseif oriY
          if slicepos(s) < ymmax && slicepos(s) > ymmin
            intersect_exists(s,k) = 1;
          else
            intersect_exists(s,k) = 0;
          end
        elseif oriZ
          if slicepos(s) < zmmax && slicepos(s) > zmmin
            intersect_exists(s,k) = 1;
          else
            intersect_exists(s,k) = 0;
          end
        end

        if intersect_exists(s,k)
          [xmesh, ymesh, zmesh] = intersect_plane(meshplot{k}.pos, meshplot{k}.tri, v1, v2, v3);

          % draw each individual line segment of the intersection
          if ~isempty(xmesh)
            p = patch(xmesh', ymesh', zmesh', nan(1, size(xmesh,1)));
            if ~isempty(intersectcolor),     set(p, 'EdgeColor', intersectcolor{k}); end
            if ~isempty(intersectlinewidth), set(p, 'LineWidth', intersectlinewidth(k)); end
            if ~isempty(intersectlinestyle), set(p, 'LineStyle', intersectlinestyle{k}); end
          end

          % find the limits of the lines and add them to the limits of the
          % interpolation to facilitate finding the limits of the slice
          xcmax(end+1) = max(xmesh(:)); xcmin(end+1) = min(xmesh(:));
          ycmax(end+1) = max(ymesh(:)); ycmin(end+1) = min(ymesh(:));
          zcmax(end+1) = max(zmesh(:)); zcmin(end+1) = min(zmesh(:));
        end
      end % for each mesh
    end % if dointersect

    % find limits of this particular slice
    xsmax(s) = max(xcmax); xsmin(s) = min(xcmin);
    ysmax(s) = max(ycmax); ysmin(s) = min(ycmin);
    zsmax(s) = max(zcmax); zsmin(s) = min(zcmin);

    % color settings
    ft_colormap(cmap);
    if ~isempty(clim) && clim(2)>clim(1)
      caxis(gca, clim);
    end

    % axis and view settings
    set(gca, 'DataAspectRatio', [1 1 1])
    if oriX
      view([90 0]);
    elseif oriY
      view([180 0]);
    elseif oriZ
      view([90 90]);
    end

    % add title to differentiate slices
    if oriX; title(['slicepos = [' num2str(slicepos(s)) ' 0 0]']); end
    if oriY; title(['slicepos = [0 ' num2str(slicepos(s)) ' 0]']); end
    if oriZ; title(['slicepos = [0 0 ' num2str(slicepos(s)) ']']); end
  end

  % set matching limits in the non-slice dimensions for each slice
  for s = 1:numel(slicepos) % slice loop
    subplot(numel(slicepos),1,s);
    if oriX
      xlim([xsmin(s)-2 xsmax(s)+2]);
      ylim([min(ysmin)-2 max(ysmax)+2]);
      zlim([min(zsmin)-2 max(zsmax)+2]);
    elseif oriY
      xlim([min(xsmin)-2 max(xsmax)+2]);
      ylim([ysmin(s)-2 ysmax(s)+2]);
      zlim([min(zsmin)-2 max(zsmax)+2]);
    elseif oriZ
      xlim([min(xsmin)-2 max(xsmax)+2]);
      ylim([min(ysmin)-2 max(ysmax)+2]);
      zlim([zsmin(s)-2 zsmax(s)+2]);
    end
  end

else % plot 3d cloud
  hold on;
  for n = 1:size(pos,1) % cloud loop
    switch (cloudtype)
      case 'cloud'
        % point cloud with radius scaling
        rng(0, 'twister'); % random number generator
        npoints   = round(ptdensity*(4/3)*pi*rmax(n)^3);      % number of points based on cloud volume
        elevation = asin(2*rand(npoints,1)-1);                % elevation angle for each point
        azimuth   = 2*pi*rand(npoints,1);                     % azimuth angle for each point
        radii     = rmax(n)*(rand(npoints,1).^ptgradient);    % radius value for each point
        radii     = sort(radii);                              % sort radii in ascending order so they are plotted from inside out
        [x,y,z]   = sph2cart(azimuth, elevation, radii);      % convert to Carthesian

        % color axis with radius scaling
        if strcmp(colorgrad, 'white')                   % color runs up to white
          indx  = ceil(cmid) + sign(colscf(n))*floor(abs(colscf(n)*cmid));
          indx  = max(min(indx,size(cmapsc,1)),1);      % index should fall within the colormap
          fcol  = cmapsc(indx,:);                       % color [Nx3]
          ptcol = [linspace(fcol(1), 1, npoints)' linspace(fcol(2), 1, npoints)' linspace(fcol(3), 1, npoints)'];
        elseif isscalar(colorgrad)                      % color runs down towards colorbar middle
          rnorm = radii/rmax(n);                        % normalized radius
          if radscf(n)>=.5                              % extreme values
            ptcol = val(n) - (flip(1-rnorm).^inv(colorgrad))*val(n); % scaled fun [Nx1]
          elseif radscf(n)<.5                           % values closest to zero
            ptcol = val(n) + (flip(1-rnorm).^inv(colorgrad))*abs(val(n)); % scaled fun [Nx1]
          end
        else
          ft_error('color gradient should be either ''white'' or a scalar determining color falloff')
        end

        % draw the points
        scatter3(x+pos(n,1), y+pos(n,2), z+pos(n,3), ptsize, ptcol, '.');

      case 'surf'
        indx  = ceil(cmid) + sign(colscf(n))*floor(abs(colscf(n)*cmid));
        indx  = max(min(indx,size(cmapsc,1)),1);  % index should fall within the colormap
        fcol  = cmapsc(indx,:);                   % color [Nx3]
        [xsp, ysp, zsp] = sphere(100);
        hs = surf(rmax(n)*xsp+pos(n,1), rmax(n)*ysp+pos(n,2), rmax(n)*zsp+pos(n,3));
        set(hs, 'EdgeColor', 'none', 'FaceColor', fcol, 'FaceAlpha', 1);

      otherwise
        ft_error('unsupported cloudtype "%s"', cloudtype);

    end % switch cloudtype
  end % end cloud loop

  if ~isempty(meshplot)
    for k = 1:numel(meshplot) % mesh loop
      ft_plot_mesh(meshplot{k}, 'facecolor', facecolor{k}, 'EdgeColor', edgecolor{k}, ...
        'facealpha', facealpha(k), 'edgealpha', edgealpha(k), 'vertexcolor', vertexcolor{k});
      material dull
    end % end mesh loop

    if dointersect % plot the outlines on the mesh
      for s = 1:numel(slicepos) % slice loop
        if oriX; ori = [1 0 0]; loc = [slicepos(s) 0 0]; end
        if oriY; ori = [0 1 0]; loc = [0 slicepos(s) 0]; end
        if oriZ; ori = [0 0 1]; loc = [0 0 slicepos(s)]; end

        % normalise the orientation vector to one
        ori = ori./sqrt(sum(ori.^2));

        % shift the location to be along the orientation vector
        loc = ori*dot(loc,ori);

        % determine three points on the plane
        inplane = eye(3) - (eye(3) * ori') * ori;
        v1 = loc + inplane(1,:);
        v2 = loc + inplane(2,:);
        v3 = loc + inplane(3,:);

        for k = 1:numel(meshplot)

          % only plot if the mesh actually intersects the plane
          xmmax = max(meshplot{k}.pos(:,1)); xmmin = min(meshplot{k}.pos(:,1));
          ymmax = max(meshplot{k}.pos(:,2)); ymmin = min(meshplot{k}.pos(:,2));
          zmmax = max(meshplot{k}.pos(:,3)); zmmin = min(meshplot{k}.pos(:,3));

          if oriX
            if slicepos(s) < xmmax && slicepos(s) > xmmin
              intersect_exists(s,k) = 1;
            else
              intersect_exists(s,k) = 0;
            end
          elseif oriY
            if slicepos(s) < ymmax && slicepos(s) > ymmin
              intersect_exists(s,k) = 1;
            else
              intersect_exists(s,k) = 0;
            end
          elseif oriZ
            if slicepos(s) < zmmax && slicepos(s) > zmmin
              intersect_exists(s,k) = 1;
            else
              intersect_exists(s,k) = 0;
            end
          end

          if intersect_exists(s,k)
            [xmesh, ymesh, zmesh] = intersect_plane(meshplot{k}.pos, meshplot{k}.tri, v1, v2, v3);

            % draw each individual line segment of the intersection
            if ~isempty(xmesh)
              p = patch(xmesh', ymesh', zmesh', nan(1, size(xmesh,1)));
              if ~isempty(intersectcolor),     set(p, 'EdgeColor', intersectcolor{k}); end
              if ~isempty(intersectlinewidth), set(p, 'LineWidth', intersectlinewidth(k)); end
              if ~isempty(intersectlinestyle), set(p, 'LineStyle', intersectlinestyle{k}); end
            end
          end
        end % for each mesh
      end % for each slice
    end % if dointersect
  end % if plotting mesh

  % axis settings
  axis off
  axis vis3d
  axis equal

  % color settings
  ft_colormap(cmap);
  if ~isempty(clim) && clim(2)>clim(1)
    caxis(gca, clim);
  end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% For each layer, determine whether we need to pass gradient back
function layer = DetermineGradientPass(layer)

% find the parent of every layer
parent = DetermineLayerParent(layer);

skipBP = zeros(length(layer),1);
for i=1:length(layer)
    if isfield(layer{i}, 'skipBP')
        skipBP(i) = layer{i}.skipBP;
    end
end

update = zeros(length(layer),1);
for i=1:length(layer)
    if isfield(layer{i}, 'update')
        update(i) = layer{i}.update;
    end
end

for i=1:length(layer)
    if IsUpdatableNode(layer{i}.name)==0; continue; end
    if skipBP(i); continue; end

    layer{i}.passGradBack = sum(update(parent{i}))>0;
end

end

author: singaxiong

repo: SignalGraph


text:
function cout = intermutual(cin1, cin2, n)

%tstoolbox/@core/intermutual
%   Syntax:
%     * intermutual(cin1,cin2,n)
%
%   Input Arguments:
%     * cin1,cin2 - core objects
%
%   Calculates the mutual information of cin1 and cin2.
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

N = dlens(cin1,1);
cout = zeros(n,1);
epsilon = 1e-9;
points1 = data(cin1);
points1 = points1 - min(points1);
points2 = data(cin2);
points2 = points2 - min(points2);
ma1 = max(points1);
ma2 = max(points2);
increment = 1/N;

for i=1:n
	partitionen = 2^i;
	pointsA = 1+floor(points1 / (ma1/(partitionen-epsilon)));
	pointsB = 1+floor(points2 / (ma2/(partitionen-epsilon)));
	% Aus den Werten sind jetzt fertige Indizes fuer die Partititonen geworden

	histA = sparse(pointsA, ones(N,1), increment);
	histB = sparse(ones(N,1), pointsB, increment);
	histAB = sparse(pointsA, pointsB, increment);
	[ind1, ind2, value] = find(histAB);
	amf = sum(value .* log2(value ./ (histA(ind1) .* histB(ind2)')));	% '

	cout(i) = amf;
end

%cout = core(amf);

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% %
% Setup default values and run PK computation. Note that the kernel will
% perform BETTER if you learn the parameters via cross validation instead of
% using these defaults.
%
% Marion Neumann (m.neumann@wustl.edu)
% %

% propagation kernel parameter
num_iterations = 3;     % number of iterations (sth small 2 or 3)

% hashing parameters
w              = 1e-5;  % bin width
distance       = 'tv';  % distance to approximately preserve

% load you dataset HERE
load('mutag_mat');

num_nodes   = size(A, 1);
num_classes = max(labels);

initial_label_distributions = accumarray([(1:num_nodes)', labels], 1, [num_nodes, num_classes]);

% create a function handle to a feature transformation. Here we will
% use label diffustion as we have fully labeled graphs.
transformation = @(features) label_diffusion(features, A);

% calculate the graph kernel using the default (linear) base kernel
K = propagation_kernel(initial_label_distributions, graph_ind, transformation, ...
                       num_iterations, ...
                       'distance', distance, ...
                       'w',        w);

% % If you want to set the BASE KERNEL to an RBF kernel instead of the default
% % linear base kernel, use the following.
% length_scale = 3;
% base_kernel = @(counts) ...
%               exp(-(squareform(pdist(counts)).^2 / (2 * length_scale^2)));
%
%           % calculate the graph kernel again using the new parameters
% K = propagation_kernel(initial_label_distributions, graph_ind, transformation, ...
%                        num_iterations, ...
%                        'distance',    distance, ...
%                        'w',           w, ...
%                        'base_kernel', base_kernel);

author: muhanzhang

repo: DGCNN

save_path: github-repos/MATLAB/muhanzhang-DGCNN