text:
function y = vl_nnsoftmaxloss(x,c,dzdy)
%VL_NNSOFTMAXLOSS CNN combined softmax and logistic loss.
%   **Deprecated: use `vl_nnloss` instead**
%
%   Y = VL_NNSOFTMAX(X, C) applies the softmax operator followed by
%   the logistic loss the data X. X has dimension H x W x D x N,
%   packing N arrays of W x H D-dimensional vectors.
%
%   C contains the class labels, which should be integers in the range
%   1 to D. C can be an array with either N elements or with dimensions
%   H x W x 1 x N dimensions. In the fist case, a given class label is
%   applied at all spatial locations; in the second case, different
%   class labels can be specified for different locations.
%
%   DZDX = VL_NNSOFTMAXLOSS(X, C, DZDY) computes the derivative of the
%   block projected onto DZDY. DZDX and DZDY have the same dimensions
%   as X and Y respectively.

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

% work around a bug in MATLAB, where native cast() would slow
% progressively
if isa(x, 'gpuArray')
  switch classUnderlying(x) ;
    case 'single', cast = @(z) single(z) ;
    case 'double', cast = @(z) double(z) ;
  end
else
  switch class(x)
    case 'single', cast = @(z) single(z) ;
    case 'double', cast = @(z) double(z) ;
  end
end

%X = X + 1e-6 ;
sz = [size(x,1) size(x,2) size(x,3) size(x,4)] ;

if numel(c) == sz(4)
  % one label per image
  c = reshape(c, [1 1 1 sz(4)]) ;
end
if size(c,1) == 1 & size(c,2) == 1
  c = repmat(c, [sz(1) sz(2)]) ;
end

% one label per spatial location
sz_ = [size(c,1) size(c,2) size(c,3) size(c,4)] ;
assert(isequal(sz_, [sz(1) sz(2) sz_(3) sz(4)])) ;
assert(sz_(3)==1 | sz_(3)==2) ;

% class c = 0 skips a spatial location
mass = cast(c(:,:,1,:) > 0) ;
if sz_(3) == 2
  % the second channel of c (if present) is used as weights
  mass = mass .* c(:,:,2,:) ;
  c(:,:,2,:) = [] ;
end

% convert to indexes
c = c - 1 ;
c_ = 0:numel(c)-1 ;
c_ = 1 + ...
  mod(c_, sz(1)*sz(2)) + ...
  (sz(1)*sz(2)) * max(c(:), 0)' + ...
  (sz(1)*sz(2)*sz(3)) * floor(c_/(sz(1)*sz(2))) ;

% compute softmaxloss
xmax = max(x,[],3) ;
ex = exp(bsxfun(@minus, x, xmax)) ;

%n = sz(1)*sz(2) ;
if nargin <= 2
  t = xmax + log(sum(ex,3)) - reshape(x(c_), [sz(1:2) 1 sz(4)]) ;
  y = sum(sum(sum(mass .* t,1),2),4) ;
else
  y = bsxfun(@rdivide, ex, sum(ex,3)) ;
  y(c_) = y(c_) - 1;
  y = bsxfun(@times, y, bsxfun(@times, mass, dzdy)) ;
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function printpdf(fname)
%PRINTPDF Prints the current figure into a pdf document

% set(gca, 'LooseInset', get(gca, 'TightInset'));
fname = [regexprep(fname, '^(.*)\.pdf$', '$1'), '.eps'];
print('-depsc', fname);
if ~system(['epstopdf ', fname])
  system(['rm ', fname]);
end

author: hosang

repo: detection-proposals


text:
clear all

dt = 0.05;
t  = 0:dt:20;

Nsamples = length(t);

Xsaved = zeros(Nsamples, 3);
Zsaved = zeros(Nsamples, 1);

for k=1:Nsamples
  r = GetRadar(dt);

  [pos vel alt] = RadarUKF(r, dt);

  Xsaved(k, :) = [pos vel alt];
  Zsaved(k)    = r;
end

PosSaved = Xsaved(:, 1);
VelSaved = Xsaved(:, 2);
AltSaved = Xsaved(:, 3);

t = 0:dt:Nsamples*dt-dt;

figure
plot(t, PosSaved)

figure
plot(t, VelSaved)

figure
plot(t, AltSaved)

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
function [patchSum] = evaluateIntegralImage(ii, row, col, delta)
% This function should calculate the sum over the patch centred at row, col
% of size patchSize of the integral image ii

 % NOTE : the integral image now has an extra row and an extra column
 % so we have to add a 1 to all our indices.
 % We can notice that when we subtract delta we should've had a minus 1 as
 % well but because we're also adding a 1 is just row-delta.
 %

 row_plus  = min(row+delta+1, size(ii,1));
 row_minus = max(row-delta, 1);
 col_plus  = min(col+delta+1, size(ii,2));
 col_minus = max(col-delta, 1);

% SOME DEBUG CODE
% disp(row_plus);
% disp(row_minus);
% disp(clu_plus);
% disp(col_minus);

 patchSum = ii(row_plus,  col_plus)     ...
         + ii(row_minus, col_minus)    ...
         - ii(row_minus, col_plus)     ...
         - ii(row_plus,  col_minus);

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function vectorRgb = dtiTrackAllFibers(samplingRateMm, dataDirectory, faThresh, faTrackThresh, exclusionRoi)
%
% Track all fibers for a particular subject. These are stored in one
% massive data file in the specified path.
%
%   dtiTrackAllFibers([samplingRateMm=2], [dataDirectory],
%   [faMaskThresh=0.25], [faTrackThresh=0.15], [exclusionRoi=[]], [hemisphere='both'])
%
%  samplingRateMm: the grid resolution for seed points, in millimeters.
%  Default is 2 mm.
%
%  dataDirectory: the path to the directory containing the tensors.nii.gz
%  file. If unspecified, a directory selection dialog will pop up.
%
%
% Author: DA

if(~exist('samplingRateMm','var')||isempty(samplingRateMm))
    samplingRateMm = 2;
end;
if(~exist('dataDirectory','var') || isempty(dataDirectory))
    dataDirectory = uigetdir('', 'Set subject bin directory');
    if(isnumeric(dataDirectory)), disp('dtiTrackAllFibers canceled.'); return; end
end;
if(~exist('faThresh','var')||isempty(faThresh))
    faThresh = 0.25;
end
if(~exist('faTrackThresh','var')||isempty(faTrackThresh))
    faTrackThresh = 0.15;
end
if(~exist('exclusionRoi','var'))
    exclusionRoi = [];
end
if(~isempty(exclusionRoi))
  if(~isstruct(exclusionRoi))
	exclusionRoi = dtiReadRoi(exclusionRoi);
  end
end

disp ('Loading tensors...');

opts.stepSizeMm = 1;
opts.faThresh = faTrackThresh;
opts.lengthThreshMm = 30;
opts.angleThresh = 50;
opts.wPuncture = 0.2;
opts.whichAlgorithm = 1;
opts.whichInterp = 1;
opts.seedVoxelOffsets = [0];

dt6 = niftiRead(fullfile(dataDirectory,'tensors.nii.gz'));
% check for the new NIFTI-compliant 5d format
if(dt6.ndim>4)
  % Convert from the 5d, lower-tri row order NIFTI tensor format
  % (Dxx Dxy Dyy Dxz Dyz Dzz) to our 4d tensor format (Dxx Dyy Dzz Dxy Dxz Dyz).
  dt6.data = double(squeeze(dt6.data(:,:,:,1,[1 3 6 2 4 5])));
end
if(exist(fullfile(dataDirectory,'brainMask.nii.gz')))
  mask = niftiRead(fullfile(dataDirectory,'brainMask.nii.gz'));
  for(ii=1:6)
	tmp = dt6.data(:,:,:,ii);
	tmp(~mask.data) = 0;
	dt6.data(:,:,:,ii) = tmp;
  end
  clear tmp mask;
end

% Create WM ROIs
[vec,val] = dtiEig(dt6.data);
fa = dtiComputeFA(val);
fa(fa>1) = 1; fa(fa<0) = 0;
if(nargout>0)
  img = abs(squeeze(vec(:,:,:,[1 2 3],1)));
  img(isnan(img)|img<0) = 0;
  img(img>1) = 1;
  for(ii=1:3) img(:,:,:,ii) = img(:,:,:,ii).*fa; end
  vectorRgb = niftiGetStruct(single(img),dt6.qto_xyz);
  clear img;
end
clear vec val;

samplingRateVox = samplingRateMm./dt6.pixdim;

% first create an enormous array of seeds and figure out which ones fall
% within the masked region:

xSeed = 1:samplingRateVox(1):size(dt6.data,1);
ySeed = 1:samplingRateVox(2):size(dt6.data,2);
zSeed = 1:samplingRateVox(3):size(dt6.data,3);

mask = fa>=faThresh;
midline = mrAnatXformCoords(dt6.qto_ijk,[0 0 0]);

[xIndex, yIndex, zIndex] = ind2sub ([size(xSeed,2), size(ySeed,2), size(zSeed,2)], ...
                                    find (mask (round(xSeed), round(ySeed), round(zSeed))));

pathsFilename = fullfile(dataDirectory, sprintf('all_paths_%dmm.pdb',samplingRateMm));

dtiWritePDBHeader (dt6.qto_xyz, pathsFilename);

fileOffsets = [];

% track in index blocks of width XX:
xStepSize = 10;
nFibers = 0;
for xIter = 1:xStepSize:size(xIndex,1)
    xRange = find (xIndex >= xIter & xIndex < xIter+xStepSize);
    seeds = [xSeed(xIndex(xRange)); ySeed(yIndex(xRange)); zSeed(zIndex(xRange))]';
    acpcSeeds = mrAnatXformCoords(dt6.qto_xyz, seeds);
    seedIndices = [xIndex(xRange), yIndex(xRange), zIndex(xRange)];
    if (isempty(seeds))
        continue;
    end;
    [fg, opts, kept] = dtiFiberTrackWrapper(dt6.data, seeds, dt6.pixdim(1:3), dt6.qto_xyz, ['FG'], opts);
    if (size(fg.fibers,1) == 0)
        continue;
    end;
    if(~isempty(exclusionRoi))
       fg = dtiIntersectFibersWithRoi([], {'not'}, [], exclusionRoi, fg);
    end
    nFibers = nFibers + length(fg.fibers);
    newOffsets = dtiAppendPathwaysToPDB (fg, pathsFilename);
    fileOffsets = [fileOffsets newOffsets];
    clear fg;
end

dtiAppendFileOffsetsToPDB (fileOffsets, pathsFilename);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function mr = mrReadDicomDir(ifileDir)
% Read a directory of GE I-files into an mr struct.
%
% mr = mrReadDicomDir(ifileDir);
%
% ras, 07/05.

% initialize mr struct
[p f ext] = fileparts(ifileDir);
if isequal(ext, '.dcm')
    % oops -- accidentally included an I-file in the path: chop this off:
    ifileDir = p;
    [p f] = fileparts(ifileDir); %#ok<ASGLU>
end

mr = mrCreateEmpty;
mr.pth = ifileDir;
mr.name = f;
mr.format = 'dicom';

[mr.data hdr1 hdr2] = loopOverMRFiles(ifileDir,'*.dcm','mrReadIfile');
mr.hdr = hdr1;

mr.voxelSize = double([mr.hdr.image.pixsize_Y mr.hdr.image.pixsize_X ...
                       mr.hdr.image.slthick]);
mr.dims = size(mr.data);
mr.spaces = mrStandardSpaces(mr);
mr.dimUnits = {'mm' 'mm' 'mm' 'sec'};
mr.dataUnits = 'T1-Weighted Intensity';

mr.info.scanner = char(mr.hdr.exam.hospname');
mr.info.subject = char(mr.hdr.exam.patname');
sex = {'female' 'male'};
mr.info.subjectSex = sex{mr.hdr.exam.patsex+1};
mr.info.subjectAge = mr.hdr.exam.patage;
mr.info.subjectSpecies = 'human'; % may not always be the case, find field
mr.info.date = datestr(parseGEDate(mr.hdr.series.se_datetime));
mr.info.scanStart = datestr(parseGEDate(mr.hdr.image.im_actual_dt));
mr.info.examNum = mr.hdr.exam.ex_no;
mr.info.protocol = char(mr.hdr.series.prtcl');
if checkfields(mr, 'hdr', 'image', 'cname')
    mr.info.coil = mr.hdr.image.cname;
else
    mr.info.coil = 'unknown coil';
end
mr.info.magnetStrength = mr.hdr.exam.magstrength/10000;

% get coordinates of 3 corners of the last image and
% express in R|A|S coords:
tmp = mr.hdr.image;
mr.info.cornerCoords = [tmp.tlhc_R tmp.trhc_R tmp.brhc_R; ...
                        tmp.tlhc_A tmp.trhc_A tmp.brhc_A; ...
                        tmp.tlhc_S tmp.trhc_S tmp.brhc_S];
tmp = hdr2.image;
mr.info.cornerCoords = [mr.info.cornerCoords; ...
                        tmp.tlhc_R tmp.trhc_R tmp.brhc_R; ...
                        tmp.tlhc_A tmp.trhc_A tmp.brhc_A; ...
                        tmp.tlhc_S tmp.trhc_S tmp.brhc_S];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use the corner coords to estimate transform into  %
% scanner space:                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % first, get the coordinates of the corners of the last slice
% % in pixel coords: (Using code Bob wrote, below)
% tlhc = [1 1 mr.dims(3)]'; % top left-hand-corner
% trhc = [1 mr.dims(2:3)]'; % top right-hand-corner
% brhc = [mr.dims]';        % bottom right-hand-corner
% pixCorners = [tlhc trhc brhc];

% now build a space defining the scanner coords, where the
% xform maps from the pixCorners to the R|A|S coords of the
% the three corners from the header:
w = dir(fullfile(ifileDir,'*.dcm'));
firstIfile = fullfile(ifileDir,w(1).name);
mr.spaces(end+1).name = 'Scanner';
mr.spaces(end).xform = inv(affineScanner2Pixels(firstIfile));
mr.spaces(end).dirLabels = {'L <--> R' 'P <--> A'  'I <--> S'};
mr.spaces(end).sliceLabels =  {'Sagittal' 'Coronal' 'Axial'};
mr.spaces(end).units = 'mm';
mr.spaces(end).coords = [];
mr.spaces(end).indices = [];

% also update the direction labels on the pixel and L/R flipped spaces,
% using the header info:
dirs = mrIfileDirections(ifileDir);
for i=1:3, mr.spaces(i).dirLabels = dirs; end
mr.spaces(3).dirLabels{2} = dimFlip(mr.spaces(3).dirLabels{2});

return
% /-----------------------------------------------------------------/ %

% /-----------------------------------------------------------------/ %
function dirText = dimFlip(dirText)
% flip direction 'a <--> b' to read 'b <--> a'.
I = strfind(dirText, ' <--> ');
if isempty(I), return; end
lhs = dirText(1:I);
rhs = dirText(I+5:end);
dirText = [rhs '<-->' lhs];
return


text:
function [M, S] = instantiate(disc)
%INSTANTIATE   Convert a COEFFSDISCRETIZATION to discrete form.
%   M = INSTANTIATE(DISC) converts each item DISC.SOURCE to discrete form
%   using the information in discretization DISC. The result M is return a cell
%   array if DISC.SOURCE has more than one component.
%
%   [M, S] = INSTANTIATE(DISC) retusn a second output, S, which is a cell array
%   containing the dscrete form of the conversion operator for each block
%   of DISC.SOURCE. Only used by ULTRAS class.
%
%   DISC.SOURCE may be one or a cell array of:
%      linBlock (becomes a matrix)
%      chebfun (becomes a vector)
%      numeric (not changed)

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

data = disc.source;
if ( isa(data, 'chebmatrix') )
    data = data.blocks;
end

if ( iscell(data) )
    M = cell(size(data));
    S = cell(size(data));
    for j = 1:size(data, 1)
        for k = 1:size(data, 2)
            discJK = extractBlock(disc, j, k);
            [M{j,k}, S{j,k}] = instantiate(discJK);
        end
    end
    return
else
    [M, S] = instantiateOne(disc, data);
end

end

function [M, S] = instantiateOne(disc, item)
% Instantiate one block of data.

if ( isa(item, 'operatorBlock') )
    % Convert a square block.

    if ( ~isempty(disc.coeffs) )
        % Coefficients of the block are available, convert to a diffmat.
        [M, S] = quasi2diffmat(disc);
    else
        [M, S] = convertOperator(disc, item);
    end

elseif ( isa(item, 'functionalBlock') )
    % Convert a row block.

    [M, S] = convertFunctional(disc, item);

elseif ( isa(item, 'chebfun') )
    % Block is a CHEBFUN. Convert to value space.

    M = toValues(disc, item);
    if ( item.isTransposed )
        M = M.';
    end
    S = zeros(size(M));

elseif ( isnumeric(item) )
    % Block is numeric, don't need to do much.

    M = item;
    S = 1;

else

    error('COEFFSDISCRETIZATION:instantiate:inputType', ...
        'Unrecognized item type.')

end

end

function [M, S] = convertFunctional(disc, item)
    % Developer note: In general we can't represent functional
    % blocks via coeffs. To get around this we instantiate a
    % VALSDISCRETIZAION and convert it to coefficient space
    % using COEFFS2VALS(). (Note it's COEFFS2VALS() rather than
    % VALS2COEFFS() because it's a right-multiply (I think..).)

    % For convenience:
    dim = disc.dimension;
    dom = disc.domain;

    % Create a TECH and a VALSDISCRETIZATION:
    tech = disc.returnTech;
    tech = tech();
    valsDisc = tech.returnValsDisc;
    valsDisc = valsDisc(item, dim, dom);
    M = matrix(valsDisc);

    % Convert from value-space to coeff-space using COEFFS2VALS.
    cumsumDim = [0, cumsum(dim)];
    tmp = cell(1, numel(dom)-1);
    for l = 1:numel(tmp)
        Ml = M(cumsumDim(l) + (1:dim(l)));
        Ml = rot90(Ml);
        tmp{l} = rot90(tech.coeffs2vals(Ml), -1);
    end

    M = cell2mat(tmp);
    S = zeros(size(M));
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [varargout]=rigidTransformationMatrixDirect(V1,V2)

%%
%Force input to 3D
if size(V1,2)==2
    V1(:,3)=0;
end

if size(V2,2)==2
    V2(:,3)=0;
end

[Q]=kabschRotationMatrix(V1,V2);

V1_m=mean(V1,1);
V2_m=mean(V2,1);

T1=eye(4,4);
T1(1:3,end)=V2_m(:);

R=eye(4,4);
R(1:3,1:3)=Q;

T2=eye(4,4);
T2(1:3,end)=-V1_m;

T=T1*R*T2;

varargout{1}=T;
varargout{2}=R(1:3,1:3);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [xdot] = secondOrderDynamics(obj, t, x, controller, params, logger)
    % calculate the dynamical equation of the second order dynamical system
    %
    % Parameters:
    % t: the time instant @type double
    % x: the states @type colvec
    % controller: the controller @type Controller
    % params: the parameter structure @type struct
    % logger: the data logger object @type SimLogger
    %
    % Return values:
    % xdot: the derivative of the system states @type colvec

    % extract the state variables into x and dx
    nx = obj.numState;
    q = x(1:nx);
    dq = x(nx+1:end);

    % store time and states into object private data for future use
    obj.t_ = t;
    obj.states_.x = q;
    obj.states_.dx = dq;

    % compute the mass matrix and drift vector (internal dynamics)
    M = calcMassMatrix(obj, q);
    Fv = calcDriftVector(obj, q, dq);

    %% get the external input
    % initialize the Gv_ext vector
    Gv_ext = zeros(nx,1);
    f_ext_name = fieldnames(obj.Inputs.External);
    if ~isempty(f_ext_name)              % if external inputs are defined
        n_ext = length(f_ext_name);

        for i=1:n_ext
            f_name = f_ext_name{i};
            % get the Gvec function object
            % g_fun = obj.Gvec.External.(f_name);
            % call the callback function to get the external input
            f_ext = obj.ExternalInputFun(obj, f_name, t, q, dq, params, logger);
            % compute the Gvec, and add it up
            Gv_ext = Gv_ext + feval(obj.GvecName_.External.(f_name),q,f_ext);

            % store the external inputs into the object private data
            obj.inputs_.External.(f_name) = f_ext;
        end
    end

    %% holonomic constraints
    h_cstr_name = fieldnames(obj.HolonomicConstraints);
    if ~isempty(h_cstr_name)           % if holonomic constraints are defined
        h_cstr = struct2array(obj.HolonomicConstraints);
        n_cstr = length(h_cstr);
        % determine the total dimension of the holonomic constraints
        cdim = sum([h_cstr.Dimension]);
        % initialize the Jacobian matrix
        Je = zeros(cdim,nx);
        Jedot = zeros(cdim,nx);

        idx = 1;
        for i=1:n_cstr
            cstr = h_cstr(i);

            % calculate the Jacobian
            [Jh,dJh] = calcJacobian(cstr,q,dq);
            cstr_indices = idx:idx+cstr.Dimension-1;
            tol = 1e-2;
            if norm(Jh*dq) > tol
                warning('The holonomic constraint %s violated.', h_cstr_name{i});
            end
            Je(cstr_indices,:) = Jh;
            Jedot(cstr_indices,:) = dJh;
            idx = idx + cstr.Dimension;
        end
    else
        Je = [];
        Jedot = [];
    end

    %% calculate the constrained vector fields and control inputs
    control_name = fieldnames(obj.Inputs.Control);
    Gv_u = zeros(nx,1);
    if ~isempty(control_name)
        Be = feval(obj.GmapName_.Control.(control_name{1}),q);
        Ie    = eye(nx);

        if isempty(Je)
            vfc = [
                dq;
                M \ (Fv)];
            gfc = [
                zeros(size(Be));
                M \ Be];
        else
            XiInv = Je * (M \ transpose(Je));
            % compute vector fields
            % f(x)
            vfc = [
                dq;
                M \ ((Ie-transpose(Je) * (XiInv \ (Je / M))) * (Fv) - transpose(Je) * (XiInv \ Jedot * dq))];

            % g(x)
            gfc = [
                zeros(size(Be));
                M \ (Ie - transpose(Je)* (XiInv \ (Je / M))) * Be];
        end
        % compute control inputs
        if ~isempty(controller)
            u = calcControl(controller, t, x, vfc, gfc, obj, params, logger);
        else
            u = zeros(size(Be,2),1);
        end
        Gv_u = Be*u;
        obj.inputs_.Control.(control_name{1}) = u;
    end
    %% calculate constraint wrench of holonomic constraints
    Gv = Gv_ext + Gv_u;
    % Calculate constrained forces
    Gv_c = zeros(nx,1);
    if ~isempty(h_cstr_name)
        lambda = -XiInv \ (Jedot * dq + Je * (M \ (Fv + Gv)));
        % the constrained wrench inputs
        Gv_c = transpose(Je)*lambda;

        % extract and store
        idx = 1;
        for i=1:n_cstr
            cstr = h_cstr(i);
            hval.(h_cstr_name{i}) = calcConstraint(cstr,q);
            cstr_indices = idx:idx+cstr.Dimension-1;
            input_name = cstr.InputName;
            obj.inputs_.ConstraintWrench.(input_name) = lambda(cstr_indices);
            idx = idx + cstr.Dimension;
        end
    else
        hval = [];
    end

    Gv = Gv + Gv_c;

    % the system dynamics
    xdot = [dq;
        M \ (Fv + Gv)];
    obj.states_.ddx = xdot(nx+1:end);

    if ~isempty(logger)
        calc = logger.calc;

        calc.t       = t;
        calc.states  = obj.states_;
        calc.inputs  = obj.inputs_;
        calc.hval     = hval;

        logger.calc  = calc;
    end
end

author: ayonga

repo: frost-dev

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


text:
classdef SubPixel_Conv < dagnn.ElementWise
%%% Sub-pixel convolution layer %%%
%
% Performs sub-pixel convolution or pixel shuffle as specified in [1].
% From an input of size [H, W, C, N] stored in inputs{1}, this layer
% produces the output of size [scale*H, scale*W, C/(scale*scale), N] in outputs{1}.
% *Back-propagation (backward function) implemented*
%
% [1] Wenzhe Shi et al. "Real-Time Single Image and Video Super-Resolution
% Using an Efficient Sub-Pixel Convolutional Neural Network", CVPR, 2016.

  properties
    scale = 2;
  end

  methods
    function outputs = forward(obj, inputs, params)
        input = inputs{1};
        output=gpuArray(zeros(size(input, 1)*obj.scale, size(input, 2)*obj.scale, size(input, 3)/obj.scale/obj.scale, size(input, 4)));
        for channel = 1:size(input, 3)
            ch = floor((channel-1)/obj.scale/obj.scale)+1;
            c = mod(channel,obj.scale*obj.scale);
            if c == 0, c = obj.scale*obj.scale; end
            q = floor((c-1)/obj.scale)+1;
            r = mod(c, obj.scale);
            if r == 0, r = obj.scale; end
            output(q:obj.scale:end, r:obj.scale:end, ch, :) = input(:, :, channel, :);
        end
        outputs{1} = gpuArray(single(output));
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
        output = derOutputs{1};
        input = gpuArray(zeros(size(output, 1)/obj.scale, size(output, 2)/obj.scale, size(output, 3)*obj.scale*obj.scale, size(output, 4)));
        for channel = 1:size(input, 3)
            ch = floor((channel-1)/obj.scale/obj.scale)+1;
            c = mod(channel,obj.scale*obj.scale);
            if c == 0, c = obj.scale*obj.scale; end
            q = floor((c-1)/obj.scale)+1;
            r = mod(c, obj.scale);
            if r == 0, r = obj.scale; end
            input(:, :, channel, :) = output(q:obj.scale:end, r:obj.scale:end, ch, :);
        end
      derInputs{1} = gpuArray(single(input));
      derParams = {} ;
    end

    function obj = SubPixel_Conv(varargin)
      obj.load(varargin) ;
      obj.scale = obj.scale;
    end
  end
end

author: sooyekim

repo: Deep-SR-ITM


text:
function K = set_params_kernel(K, name, val)
% SET_PARAMS_KERNEL Accessor function for a field (table_kernel)
% K = set_params_kernel(K, name, val)
%
% e.g., K = set_params_kernel(K, 'table', rand(2,3,2)) for a kernel on 3 nodes with 2,3,2 values each

% We should check if the arguments are valid...

switch name
 case 'table', K.table = val;
 otherwise,
  error(['invalid field name ' name]);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function varargout = scatter(q,varargin)
% plot function
%
% Input
%  q - Quaternion
%
% Options
%  rodrigues - plot in rodrigues space
%  axisAngle - plot in axis / angle
%

oP = newOrientationPlot(specimenSymmetry,specimenSymmetry,varargin{:});

[varargout{1:nargout}] = oP.plot(q,varargin{:});

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function plotChannels(imgIn, imgOut)
%
%
%        plotChannels(imgIn, imgOut)
%
%        This function plots colors channels
%
%     Copyright (C) 2015  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

[r_i, c_i, col_i] = size(imgIn);
[r_o, c_o, col_o] = size(imgOut);

if(col_o ~= col_i)
    error('different images!');
end

hold on;
figure(1);

for i=1:col_i
    tmpX = imgIn(:,:,i);
    tmpX = imresize(tmpX, [16, 16], 'bilinear');
    tmpX = tmpX(:);

    tmpY = imgOut(:,:,i);
    tmpY = imresize(tmpY, [16, 16], 'bilinear');
    tmpY = tmpY(:);

    [tmpX, ind] = sort(tmpX(:), 'ascend');
    windowSize = 16;
    tmpY = tmpY(ind);
    tmpY = filter( (1/windowSize)*ones(1,windowSize), 1, tmpY);

    plot(tmpX, tmpY);
end

end

author: banterle


text:
function add_sampling_menu(obj)
    % add grid menu for modifying grid in a ZmapMainWindow
    % obj is a ZmapMainWindow object.
    parent = uimenu(obj.fig,'Label','Sampling');
    MenuSelectedFcn='MenuSelectedFcn';

    %% Grid menu items

    uimenu(parent,'Label','Quick-Grid (auto)',MenuSelectedFcn,@cb_autogrid);
    uimenu(parent,'Label','Define Grid',MenuSelectedFcn,@cb_gridfigure);
    uimenu(parent,'Label','Redraw Grid',MenuSelectedFcn,@cb_refresh);
    uimenu(parent,'Label','Clear Grid (Delete)',MenuSelectedFcn,@cb_clear);

    %% Sampling menu items
    uimenu(parent,'Separator','on','Label','Choose Sample Radius',MenuSelectedFcn,@cb_manualradius);
    XYfun.sample_preview.AddMenuItem(parent, @()obj.map_zap);

    %% Polygon menu items
    uimenu(parent,'Separator','on',...
        'Label','Select events in CIRCLE',MenuSelectedFcn,@cb_makecircle);
    uimenu(parent,'Label','Select events in BOX', MenuSelectedFcn,@cb_makebox);
    uimenu(parent,'Label','Select events in POLYGON', MenuSelectedFcn,@cb_makepolygon);
    uimenu(parent,'Label','Wrap events in POLYGON (complex hull)', MenuSelectedFcn,@cb_makehull);
    uimenu(parent,'Label','about editing polygons...',MenuSelectedFcn,@(~,~)moveable_item('help'));
    uimenu(parent,'Separator','on',...
        'Label','Delete polygon', MenuSelectedFcn, @cb_clear_shape);

    shapeiomenu=uimenu(parent,'Separator','on','Label','Polygon IO...');
    uimenu(shapeiomenu,'Label','get default', MenuSelectedFcn, @(~,~)cb_get_default_shape);
    uimenu(shapeiomenu,'Label','set as default', MenuSelectedFcn, @(~,~)ShapeGeneral.ShapeStash(obj.shape));
    uimenu(shapeiomenu,'Separator','on',...
        'Label','Load a polygon','MenuSelectedFcn',@cb_load_shape)
    uimenu(shapeiomenu,'Label','Save a polygon','MenuSelectedFcn',@cb_save_shape)
    % uimenu(shapeiomenu,'Label','save', MenuSelectedFcn, @(~,~)obj.shape.save(ZmapGlobal.Data.Directories.data));

    function cb_makecircle(src,ev)
        bringToForeground(obj.map_axes);
        %bringToForeground(findobj(obj.fig,'Tag','mainmap_ax'));
        sh=ShapeCircle.selectUsingMouse(obj.map_axes, obj.refEllipsoid);
        set_my_shape(obj,sh);
    end

    function cb_makebox(src,ev)
        bringToForeground(obj.map_axes);
        %bringToForeground(findobj(obj.fig,'Tag','mainmap_ax'));
        sh=ShapePolygon('box');
        set_my_shape(obj,sh);
    end

    function cb_makepolygon(src,ev)
        bringToForeground(obj.map_axes);
        %bringToForeground(findobj(obj.fig,'Tag','mainmap_ax'));
        sh=ShapePolygon('polygon');
        set_my_shape(obj,sh);
    end

    function cb_makehull(src,ev)
        bringToForeground(obj.map_axes);
        %bringToForeground(findobj(obj.fig,'Tag','mainmap_ax'));
        eqs = [obj.catalog.X, obj.catalog.Y];
        ch=convhull(eqs,'simplify',true);
        sh=ShapePolygon('polygon',eqs(ch,:));
        set_my_shape(obj,sh);
    end
    function cb_clear_shape(src,ev)
        ShapeGeneral.clearplot();
        delete(obj.shape);
        obj.shape=ShapeGeneral();
        %obj.replot_all();
    end

    function cb_save_shape(src,ev)
        if isempty(obj.shape)
            errordlg('No shape is currently selected.');
        else
            obj.shape.save();
        end
    end
    function cb_load_shape(src,ev)
        bringToForeground(obj.map_axes);
        %bringToForeground(findobj(obj.fig,'Tag','mainmap_ax'));
        sh = load_shape();
        if isempty(sh)
            errordlg('Unable to load shape, or operation was cancelled');
        else
            obj.set_my_shape(sh);
        end
    end

    function cb_get_default_shape(src,ev)
        sh=ShapeGeneral.ShapeStash();
        if isempty(sh)
            warndlg('No default shape exists');
        else
            set_my_shape(obj,sh);
        end
    end

    function cb_autogrid(~,~)
        % following assumes grid from main map

        if ~isempty(obj.Grid)
            todel=findobj(obj.map_axes,'Tag',['grid_', obj.Grid.Name]);
        else
            todel=[];
        end
        delete(todel);

        [tmpgrid,obj.gridopt]=autogrid(obj.catalog, obj.refEllipsoid, obj.map_axes);
        obj.Grid = tmpgrid.MaskWithShape(obj.shape);
        obj.Grid.plot(obj.map_axes,'ActiveOnly');

    end

    function cb_gridfigure(src,ev)
        watchon
        drawnow
        [gr, gro] = GridOptions.fromDialog(obj.gridopt, obj.refEllipsoid,obj.shape);
        if ~isempty(gr)
            obj.Grid = gr;
            obj.gridopt = gro;
        end
        mygr=findobj(obj.map_axes.Children,'flat','-regexp','Tag','grid_\w.*');
        delete(mygr);
        obj.Grid.plot(obj.map_axes,'ActiveOnly');
        watchoff
    end

    function cb_manualradius(~,~)
        ev = obj.get_event_selection;
        if isempty(ev)
            [evselch, okpressed] = EventSelectionChoice.quickshow(true);
        else

            [evselch, okpressed] = EventSelectionChoice.quickshow(false, ev);
        end
        if okpressed
            obj.set_event_selection(EventSelectionParameters.fromStruct(evselch));
        end

        % send message that global grid changed
    end

    function cb_refresh(~,~)
        if isempty(obj.Grid)
            warning('ZMAP:grid:missingGrid','no grid exists to refresh')
            warndlg('no grid exists to refresh')
            return
        end
        delete(findobj(obj.fig,'Tag',['grid_',obj.Grid.Name]))
        obj.Grid=obj.Grid.MaskWithShape(obj.shape);
        obj.Grid.plot(obj.map_axes,'ActiveOnly');
    end

    function cb_clear(~,~)
        try
            obj.Grid = obj.Grid.delete();
        catch ME
            warning(ME.message)
        end
    end

end

author: CelsoReyes


text:
classdef GaussianD
%%GAUSSIAND Functions to handle the scalar and multivariate Gaussian
%           distribution.
%Implemented methods are: mean,cov, PDF, PDFI, PDFIDerivs (for multivariate
%                         derivatives of the PDF), PDFS, logPDFS,
%                         PDFSGradHessVechS (for the gradient and Hessian
%                         of the elements of a lower-triangular square-root
%                         of the covariance matrix), CDF (for scalar
%                         distributions), invCDF (for scalar
%                         distributions), normProdDist, normConvDist (for
%                         scalar distributions), momentGenFun
%                         (multivariate, including derivatives), cumGenFun
%                         (multivariate, including derivatives), rand,
%                         randS, integralOverRegion, entropy
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)

function val=mean(mu)
%%MEAN Obtain the mean of the Gaussian distribution.
%
%INPUTS: mu The mean of the PDF. If the PDF is multivariate, then this is
%           a numDimX1 column vector.
%
%OUTPUTS: val The numDimX1 mean of the Gaussian distribution.
%
%The Gaussian distribution is parameterized by its mean and covariance
%matrix. Thus, this function just returns the mean it is given.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=mu;
end

function val=cov(Sigma)
%%COV Obtain the covariance matrix of the Gaussian distribution (the
%     variance if scalar).
%
%INPUTS: Sigma The variance (if scalar) or covariance matrix (if
%               multidimensional) of the PDF. The variance cannot be zero
%               and the covariance matrix cannot be singular.
%
%OUTPUTS: val The covariance matrix of the Gaussian distribution.
%
%The Gaussian distribution is parameterized by its mean and covariance
%matrix. Thus, this function just returns the covariance matrix it is
%given.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=Sigma;
end

function vals=PDF(z,mu,Sigma)
%%PDF Evaluate a scalar or multivariate Gaussian (normal) PDF at specified
%     points given the mean and the covariance matrix.
%
%INPUTS: z The points at which the PDF should be evaluated. If the PDF is
%          multivariate, then this is a column vector. If evaluation at
%          multiple points are desired, then this is a numDimXN matrix with
%          each column being the a point (a vector).
%       mu The mean of the PDF. If the PDF is multivariate, then this is a
%          numDimX1 column vector. If omitted or an empty matrix is passed,
%          a zero mean is used.
%    Sigma The variance (if scalar) or numDimXnumDim covariance matrix
%          (if multidimensional) of the PDF. The variance cannot be zero
%          and the covariance matrix cannot be singular. If omitted or an
%          empty matrix is passed, the identity matrix is used as the
%          covariance matrix.
%
%OUTPUTS: vals The scalar values of the normal PDF with mean mu and
%              covariance matrix Sigma evaluated at the points in z. If
%              multiple points are passed (z is a matrix), then val is a
%              row vector.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    numDim=size(z,1);
    if(nargin<2||isempty(mu))
        mu=zeros(numDim,1);
    end

    if(nargin<3||isempty(Sigma))
       Sigma=eye(numDim,numDim);
    end

    diff=bsxfun(@minus,z,mu);
    vals=(2*pi)^(-numDim/2)*(det(Sigma))^(-1/2)*exp(-0.5*invSymQuadForm(diff,Sigma));
end

function vals=logPDF(z,mu,Sigma)
%%LOGPDF Evaluate the natural logarithm of a scalar or multivariate
%        Gaussian (normal) PDF at a certain points given the mean and the
%        covariance matrix.
%
%INPUTS: z The points at which the PDF should be evaluated. If the PDF is
%          multivariate, then this is a column vector. If evaluation at
%          multiple points are desired, then this is a numDimXN matrix with
%          each column being the a point (a vector).
%       mu The mean of the PDF. If the PDF is multivariate, then this is a
%          numDimX1 column vector. If omitted or an empty matrix is passed,
%          a zero mean is used.
%    Sigma The variance (if scalar) or numDimXnumDim covariance matrix
%          (if multidimensional) of the PDF. The variance cannot be zero
%          and the covariance matrix cannot be singular. If omitted or an
%          empty matrix is passed, the identity matrix is used as the
%          covariance matrix.
%
%OUTPUTS: vals The scalar values of the natural logarithm of the normal PDF
%              with mean mu and covariance matrix Sigma evaluated at the
%              points in z. If multiple points are passed (z is a matrix),
%              then val is a row vector.
%
%March 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.

    numDim=size(z,1);
    if(nargin<2||isempty(mu))
        mu=zeros(numDim,1);
    end

    if(nargin<3||isempty(Sigma))
       Sigma=eye(numDim,numDim);
    end

    diff=bsxfun(@minus,z,mu);
    vals=-(1/2)*log(det(2*pi*Sigma))-(1/2)*invSymQuadForm(diff,Sigma);
end

function vals=PDFI(z,mu,SigmaInv,SigmaInvDet)
%%PDFI Evaluate a scalar or multivariate Gaussian (normal) PDF at specified
%      points given the mean and the inverse of the covariance matrix.
%
%INPUTS: z The points at which the PDF should be evaluated. If the PDF is
%          multivariate, then this is a column vector. If evaluations at
%          multiple points are desired, then this is a numDimXN matrix with
%          each column being the a point (a vector).
%       mu The mean of the PDF. If the PDF is multivariate, then this is a
%          numDimX1 column vector. If omitted or an empty matrix is passed,
%          a zero mean is used.
% SigmaInv The inverse variance (if scalar) or numDimXnumDim inverse
%          covariance matrix (if multidimensional) of the PDF. SigmaInv can
%          be singular. If omitted or an empty matrix is passed, the
%          identity matrix is used as the covariance matrix.
% SigmaInvDet Optionally, a length-N set of determinants of the matrices
%          in SigmaInv can be passed so as to speed up the computation. If
%          omitted or an empty matrix is passed, determinants will be taken
%          as needed.
%
%OUTPUTS: val The scalar value of the normal PDF with mean mu and inverse
%             covariance matrix SigmaInv evaluated at the point z. If
%             multiple points are passed (z is a matrix), then val is a row
%             vector.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    numPoints=size(z,2);
    vals=zeros(1,numPoints);
    n=size(z,1);
    if(nargin<2||isempty(mu))
        mu=zeros(numDim,1);
    end

    if(nargin<3||isempty(SigmaInv))
       SigmaInv=eye(numDim,numDim);
    end

    if(nargin<4||isempty(SigmaInvDet))
        constVal=sqrt(det(SigmaInv)/(2*pi)^(n));
    else
        constVal=sqrt(SigmaInvDet/(2*pi)^(n));
    end
    for curPoint=1:numPoints
        diff=z(:,curPoint)-mu;
        %Note that det(A^(-1))=1/det(A) and that det(a*A)=a^n*det(A), where
        %a is a scalar and A is an nXn matrix.

        vals(curPoint)=constVal*exp(-0.5*(diff'*SigmaInv*diff));
    end
end

function [PDFDerivVal,coeffPolyPart]=PDFIDerivs(mu,SigmaInv,numDerivs,x)
%%PDFIDERIVS Compute derivatives of the multivariate Gaussian normal PDF at
%            given the mean and the inverse of the covariance matrix.
%            Derivatives are taken with respect to components of the
%            argument of the PDF, x.
%
%INPUTS: mu The mean of the PDF. If the PDF is multivariate, then this is a
%           numDimX1 column vector.
%  SigmaInv The inverse variance (if scalar) or numDimXnumDim inverse
%           covariance matrix (if multidimensional) of the PDF.
%           SigmaInv can be singular.
% numDerivs A numDimX1 or 1XnumDim vector indicating the number of
%           derivatives to take with respect to each of the dimensions of
%           the state.numDerivs>=0.
%         x The numDimXnumPoints argument of the PDF at which the
%           derivatives of the PDF should be evaluated. If this parameter
%           is omitted or an empty matrix is passed, then a default of
%           x=zeros(numDim,1) is used.
%
%OUTPUTS: PDFDerivVal A numPointsX1 vector of the values of the derivatives
%                   of the PDF function given at the points in x or at x=0
%                   if x is omitted.
%     coeffPolyPart A hypermatrix taking numDim indices that can be
%                   evaluated using the polyValMultiDim at different values
%                   of x to get the polynomial coefficient that multiplies
%                   the exponential term in the given set of derivatives
%                   of the PDF.
%
%All derivatives of the multivariate normal distribution have the form
%sqrt((2*pi)^(-2)*det(SigmaInv))*exp(-1/2*((x-mu)'*SigmaInv*(x-mu))*...
%(polynomial)
%The polynomial can be found using the chain rule each time a derivative is
%taken. The derivative of the exponential term with respect to x(i), the
%ith component of the state is
%sqrt((2*pi)^(-2)*det(SigmaInv))*exp(-1/2*((x-mu)'*SigmaInv*(x-mu)) times
%SigmaInv(i,:)*mu-SigmaInv(i,:)*x.
%This function returns coeffPolyPart, the coefficients of the multivariate
%polynomial that multiplies
%sqrt((2*pi)^(-2)*det(SigmaInv))*exp(-1/2*((x-mu)'*SigmaInv*(x-mu)) when
%computing the derivatives in addition to returning the value of the PDF at
%any desired points.
%
%December 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

numIdx=size(mu,1);

coeffPolyPart=1;

%i is the current dimension we are differentiating.
for i=1:numIdx
    %basePoly shall be the multivariate polynomial that is added every time
    %the exponential term in the moment generating function is
    %differentiated with respect to the ith index. This makes basePoly the
    %derivative with respect to the ith element of t of the argument of the
    %exponent in the multivariate normal PDF.

    %Allocate space for the polynomial and make it have the correct shape.
    basePoly=reshape(zeros(2^numIdx,1),2*ones(1,numIdx));

    %The additive term
    basePoly(1)=SigmaInv(i,:)*mu;

    %The term multiplied by x consists of elements from the inverse
    %covariance matrix. This just consists of terms in the ith row of the
    %inverse covariance matrix (the matrix is symmetric).
    idxVec=ones(numIdx,1);
    %The dimensionalities of all of the variables for the nDim2Index
    %function.
    dims=2*ones(numIdx,1);
    for curDim=1:numIdx
        idxVec(curDim)=2;
        basePoly(nDim2Index(dims,idxVec))=-SigmaInv(i,curDim);
        idxVec(curDim)=1;
    end

    %Now, we enter into a loop to evaluate derivatives of the PDF with
    %respect to the current dimension.
    for derivsLeft=numDerivs(i):-1:1
        %The chain rule means that there are two terms to consider (both of
        %which are multiplied by the same exponential term).
        %The first term is the current coeffPolyPart times the derivative
        %of the exponential term with respect to the ith dimensions. This
        %is the product of coeffPolyPart and basePoly. The second term is
        %the exponential term times the derivative of coeffPolyPart. The
        %two terms then must be added.
        term1=convn(basePoly,coeffPolyPart);%Multiply the polynomials.
        term2=polyDerMultiDim(coeffPolyPart,i);

        %Add the multivariate polynomials.
        coeffPolyPart=polySumMultiDim(term1,term2);
    end
end

%Get rid of redundant parts that arose due to the multiplication.
coeffPolyPart=shrinkMultiDimPoly2Fit(coeffPolyPart);

%coeffPolyPart now contains the multivariate polynomial that is multiplied
%by the exponential term. If no value of t is given, then just evaluate it
%at x=0. In this instance, the exponential term is zero and only the
%constant term from the polynomial appears.
if(nargin<4||isempty(x))
    PDFDerivVal=GaussianD.PDFI([0;0],mu,SigmaInv)*coeffPolyPart(1);
else
    numPoints=size(x,2);
    PDFDerivVal=zeros(numPoints,1);
    for curPoint=1:numPoints
        xCur=x(:,curPoint);
        PDFDerivVal(curPoint)=GaussianD.PDFI(xCur,mu,SigmaInv)*polyValMultiDim(coeffPolyPart,xCur);
    end
end
end

function val=PDFS(z,mu,S)
%%PDFS Evaluate a scalar or multivariate Gaussian (normal) PDF at specifed
%      points given the mean and the lower-triangular square root of the
%      covariance matrix.
%
%INPUTS: z The points at which the PDF should be evaluated. If the PDF is
%          multivariate, then this is a column vector. If evaluations at
%          multiple points are desired, then this is a numDimXN matrix with
%          each column being the a point (a vector).
%       mu The mean of the PDF. If the PDF is multivariate, then this is a
%          numDimX1 column vector.
%        S The square root of the variance (if scalar) or the numDimXnumDim
%          lower-triangular square root of the covariance matrix (if
%          multidimensional) of the PDF such that S*S'=Sigma, where Sigma
%          is the covariance matrix. S cannot be a singular matrix. If
%          omitted or an empty matrix is passed, the identity matrix is
%          used.
%
%OUTPUTS: val The scalar value(s) of the normal PDF with mean mu and square
%             root covariance matrix S evaluated at the points in z. If
%             multiple points are passed (z is a matrix), then val is a row
%             vector.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<2||isempty(mu))
        mu=zeros(numDim,1);
    end

    if(nargin<3||isempty(S))
    	S=eye(numDim,numDim);
    end

%Note that (S*S')^(-1)=(S')^(-1)*S^(-1)
    diff=S\bsxfun(@minus,z,mu);
%Note that det(S*S')=det(S)*det(S') and that det(S)=det(S') so
%det(S*S')=det(S)^2. Also, det(a*S)=a^ndet(S), where a is a scalar and S is
%an nXn matrix. Thus,
%det(2*pi*S*S')=det(sqrt(2*pi)*S)^2=(2*pi)^n*det(S)^2
    n=size(z,1);
    %The abs in the determinant is necessary if the main diagonal of S has
    %negative terms. S can still be such that S*S'=Sigma as the sign of
    %those terms is not unique due to the squaring.
    val = (1/((2*pi)^(n/2)*abs(det(S))))*exp(-0.5*sum(diff.*diff,1));
end

function val=logPDFS(z,mu,S)
%%LOGPDFS Evaluate the natural logarithm of a scalar or multivariate
%         Gaussian (normal) PDF at specified points given the mean and the
%         lower-triangular square root of the covariance matrix.
%
%INPUTS: z The points at which the PDF should be evaluated. If the PDF is
%          multivariate, then this is a column vector. If evaluations at
%          multiple points are desired, then this is a numDimXN matrix with
%          each column being the a point (a vector).
%       mu The mean of the PDF. If the PDF is multivariate, then this is a
%          numDimX1 column vector.
%        S The square root of the variance (if scalar) or the numDimXnumDim
%          lower-triangular square root of the covariance matrix (if
%          multidimensional) of the PDF such that S*S'=Sigma, where Sigma
%          is the covariance matrix. S cannot be a singular matrix. If
%          omitted or an empty matrix is passed, the identity matrix is
%          used.
%
%OUTPUTS: val The scalar value(s) of the natural logarithm of the normal
%             PDF with mean mu and square root covariance matrix S
%             evaluated at the points in z. If multiple points are passed
%             (z is a matrix), then val is a row vector.
%
%March 2019 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<2||isempty(mu))
        mu=zeros(numDim,1);
    end

    if(nargin<3||isempty(S))
    	S=eye(numDim,numDim);
    end

%Note that (S*S')^(-1)=(S')^(-1)*S^(-1)
    diff=S\bsxfun(@minus,z,mu);
%Note that det(S*S')=det(S)*det(S') and that det(S)=det(S') so
%det(S*S')=det(S)^2. Also, det(a*S)=a^ndet(S), where a is a scalar and S is
%an nXn matrix. Thus,
%det(2*pi*S*S')=det(sqrt(2*pi)*S)^2=(2*pi)^n*det(S)^2
    n=size(z,1);
    %The abs in the determinant is necessary if the main diagonal of S has
    %negative terms. S can still be such that S*S'=Sigma as the sign of
    %those terms is not unique due to the squaring.
    % val=-(n/2)*log((2*pi))-log(abs(det(S)))-0.5*sum(diff.*diff,1);
    val=log((1/((2*pi)^(n/2)*abs(det(S))))*exp(-0.5*sum(diff.*diff,1)));
end

function [grad,CDetGrad,Hess,CDetHess]=PDFSGradHessVechS(x,mu,C)
%%PDFGRADHESSVECHS Find the gradient and (if requested Hessian) of the
%            normal (Gaussian) probability density function (PDF), taken
%            with respect to the vech(C), where C is the lower-triangular
%            square root of the covariance matrix of the distribution.
%
%INPUTS: x The dXnumPoints points at which the normal PDF is considered.
%       mu The dX1 mean of the normal PDF.
%        C The dXd lower-triangular square root covariance matrix of the
%          normal PDF. This cannot be singular.
%
%OUTPUTS: grad The gradient of the normal PDF with respect to vech(C)
%              (vector of first partial derivatives  with respect to the
%              elements of vech(C)). This is an (n*(n+1)/2)XnumPoints set
%              of vectors.
%     CDetGrad The (n*(n+1)/2)X1 gradient of 1/sqrt(C*C') with respect to
%              vech(C). This term is needed to compute grad and Hess and is
%              often needed in algorithms using grad and Hess.
%         Hess The Hessian of the normal PDF with respect to vech(C)
%              (vector of second partial derivatives  with respect to the
%              elements of vech(C)). This is an
%              (n*(n+1)/2)X(n*(n+1)/2)XnumPoints set of symmetric matrices.
%              Element i,j in a matrix is the second derivative with
%              respect to elements i and j of vech(C).
%     CDetHess The (n*(n+1)/2)X1 Hessian of 1/sqrt(C*C') with respect to
%              vech(C). This term is needed to compute grad and Hess and is
%              often needed in algorithms using grad and Hess.
%
%The gradient and Hessian provided by this function play a role in
%multivariate kernel bandwidth estimation algorithms, such as the
%cross-validation bandwidth estimation algorithms of [1], which optimize
%over the components of a lower-triangular square root matrix to avoid
%obtaining invalid covariance matrix estimates.
%
%This function is implemented based on the rules of matrix calculus.
%
%A 5-dimensional example.
% c=[17;23;4;10;11;5;6;12;18;13;19;25;21;2;9];
% C=vech2Mat(c,0);
% x=[5;-10;33;24;-18];
% mu=[0;0;0;0;0];
% [grad,Hess]=GaussianD.PDFSGradHessVechS(x,mu,C)
%
%REFERENCES:
%[1] T. Duong and M. L. Hazelton, "Cross-validation bandwidth matrices for
%    multivariate kernel density estimation," Scandinavian Journal of
%    Statistics, vol. 32, no. 3, pp. 485-506, Sep. 2005.
%
%January 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Since the derivatives are not with respect to x, it does not matter if the
%distribution is not zero-mean; just shift the value.
x=bsxfun(@minus,x,mu);

%Problem dimensionality
d=size(x,1);
numPoints=size(x,2);

vechC=vech(C);
numVechEls=length(vechC);

CDetRoot=1/sqrt(det(C*C'));

%Find the gradient of 1/sqrt(C*C') with respect to vech(C). Only the
%derivatives corresponding to the diagonals are non-zero.
CDetGrad=vech(diag(-CDetRoot./diag(C)));

CInv=inv(C);

%The exponential term
expTerm=exp(-0.5*invSymQuadForm(x,C,1));

%Find the gradient of C^(-1) with respect to vech(C).
CInvGrad=zeros(d,d,numVechEls);
CDeriv=zeros(d,d);
curEl=1;
for j=1:d
    for i=j:d
        CDeriv(i,j)=1;

        CInvGrad(:,:,curEl)=CInv*CDeriv*CInv;
        curEl=curEl+1;
        CDeriv(i,j)=0;
    end
end

%Find the gradient of (C*C')^(-1) with respect to vech(C).
CCpInvGrad=zeros(d,d,numVechEls);
for i=1:numVechEls
    temp=CInv'*CInvGrad(:,:,i);
    CCpInvGrad(:,:,i)=temp+temp';
end

%Find the gradient of the exponential term with respect to vech(C).
expTermGrad=zeros(numVechEls,numPoints);
for i=1:numVechEls
    expTermGrad(i,:)=(1/2)*expTerm.*sum(bsxfun(@times,x,CCpInvGrad(:,:,i)*x),1);
end

%Now, find the derivatives of the PDF with respect to every element of
%of vech(C).
grad=zeros(numVechEls,numPoints);
for i=1:numVechEls
    grad(i,:)=CDetGrad(i)*expTerm+CDetRoot*expTermGrad(i,:);
end
grad=grad/(2*pi)^(d/2);

if(nargout>2)%If the Hessian is desired.
    Hess=zeros(numVechEls,numVechEls,numPoints);
    CDetHess=zeros(numVechEls,numVechEls);
    for m=1:numVechEls
        CDerivM=zeros(d,d);
        [i,j]=vechInd2Sub(d,m);%--derivative indices
        CDerivM(i,j)=1;

        %This term is used in the Hessian of exp((1/2)*x'*inv(C*C')*x)
        term1H=sum(bsxfun(@times,x, CCpInvGrad(:,:,m)*x),1);
        for n=m:numVechEls
            CDerivN=zeros(d,d);
            [k,l]=vechInd2Sub(d,n);%--derivative indices
            CDerivN(k,l)=1;

            %The Hessian of inv(C) with respect to indices i,j and k,l.
            CInvHess=CInv*(CDerivN*CInv*CDerivM+CDerivM*CInv*CDerivN)*CInv;

            term1=CInv'*CInvHess;
            term2=CInvGrad(:,:,m)'*CInvGrad(:,:,n);
            %The Hessian of inv(C*C') with respect to indices i,j and k,l.
            CCpInvHess=term1+term1'+term2+term2';
            %The Hessian of exp((1/2)*x'*inv(C*C')*x)
            term2=sum(bsxfun(@times,x, CCpInvGrad(:,:,n)*x),1);
            term3=sum(bsxfun(@times,x,CCpInvHess*x),1);
            expTermHess=(expTerm/2).*((1/2)*term1H.*term2-term3);

            Hess(m,n,:)=expTermGrad(m,:)*CDetGrad(n)+expTermGrad(n,:)*CDetGrad(m)+CDetRoot*expTermHess;
            if(i==j&&l==k)%The second derivative term of the Hessian
                CDetHess(m,n)=CDetRoot/(C(i,i)*C(k,k));
                if(i==k)
                    CDetHess(m,n)=CDetHess(m,n)*2;
                end
                %Due to the symmetry of the Hessian.
                CDetHess(n,m)=CDetHess(m,n);

                Hess(m,n,:)=Hess(m,n,:)+reshape(expTerm*CDetHess(m,n),1,1,numPoints);
            end

            Hess(m,n,:)=Hess(m,n,:)*(2*pi)^(-d/2);
            %Due to the symmetry of the Hessian, the upper triangular part
            %is also known.
            Hess(n,m,:)=Hess(m,n,:);
        end
    end
end

end

function val=CDF(z,mu,varVal)
%%CDF Evaluate cumulative distribution function (CDF) of a a scalar
%     Gaussian (normal) distribution at a specified points given the mean
%     and the variance, or for a normal(0,1) distribution if the mean and
%     variance are omitted.
%
%INPUTS: z A matrix of the point(s) at which the CDF should be evaluated.
%       mu The mean of the distribution. If omitted or an empty matrix is
%          passed, a mean of 0 is used.
%   varVal The variance of the distribution. If omitted or an empty matrix
%          is passed, a variance of 1 is used.
%
%OUTPUTS: val The scalar value(s) of the normal CDF with mean mu and
%             variance varVal evaluated at the point(s) z.
%
%This just uses the relation between the normal CDF and the error function
%along with the erf function in Matlab.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

	if(nargin<2||isempty(mu))
        mu=0;
	end

	if(nargin<3||isempty(varVal))
        varVal=1;
	end

    x=(z-mu)/sqrt(varVal);
    val=(1+erf(x/sqrt(2)))/2;
end

function val=invCDF(prob,mu,varVal)
%%CDF Evaluate the inverse CDF of a scalar Gaussian (normal) distribution
%     at a given point given the mean and the variance, or for a
%     normal(0,1) distribution if the mean and variance are omitted. When
%     considering a normal(0,1) distribution, this is also known as the
%     probit function.
%
%INPUTS: prob The probability or probabilities (0<=prob<=1) at which the
%             argument of the CDF is desired.
%          mu The mean of the distribution. If omitted or an empty matrix
%             is passed, a mean of 0 is used.
%      varVal The variance of the distribution. If omitted or an empty
%             matrix is passed, a variance of 1 is used.
%
%OUTPUTS: val The argument(s) of the CDF that would give the probability or
%             probabilities in prob.
%
%This just uses the relation between the normal CDF and the error function
%along with the erfinv (inverse error function) command in Matlab.
%
%June 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<2||isempty(mu))
       mu=0;
    end

    if(nargin<3||isempty(varVal))
       varVal=1;
    end

    val=sqrt(2*varVal)*erfinv(2*prob-1)+mu;
end

function [mu,SigmaInv,multiConst]=normProdDist(mu1,SigmaInv1,mu2,SigmaInv2)
%%NORMPRODDIST The product of two multivariate normal distributions is an
%              unnormalized Gaussian distribution. This finds the
%              parameters of the product distribution.
%
%INPUTS: mu1,SigmaInv1 The mean and inverse of the covariance matrix
%                      (inverse of the variance for a scalar distribution)
%                      of the first normal distribution.
%        mu2,SigmaInv2 The mean and inverse of the covariance matrix
%                      (inverse of the variance for a scalar distribution)
%                      of the second normal distribution.
%
%OUTPUTS: mu, SigmaInv The mean and inverse covariance matrix of the
%                      product distribution.
%           multiConst The product distribution is not normalized. This is
%                      the multiplicative constant that is multiplied by a
%                      normalized distribution.
%
%The derivation of the product of two normal PDFs is a standard exercise in
%many statistics classes. The product distribution is also given
%explicitly in [1].
%
%Note that while SigmaInv1 and SigmaInv2 can each be singular, the sum must
%be non-singular. Note that for products involving distributions with
%distant means and small variances, multiConst might be numerically zero.
%
%REFERENCES:
%[1] K. B. Petersen and M. S. Pedersen, "The matrix cookbook," Technical
%    University of Denmark, Tech. Rep., 15 Nov. 2012. [Online]. Available:
%    http://www2.imm.dtu.dk/pubdb/views/publication details.php?id=3274
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

SigmaInv=SigmaInv1+SigmaInv2;
mu=SigmaInv\(SigmaInv1*mu1+SigmaInv2*mu2);

multiConst=GaussianD.PDFI(mu1,mu2,SigmaInv);
end

function [mu,Sigma]=normConvDist(mu1,Sigma1,mu2,Sigma2)
%%NORMCONVDIST The convolution of two scalar normal distributions is also a
%              normal distribution. This provides the parameters for the
%              convoluted distribution.
%
%INPUTS: mu1,Sigma1 The mean and variance of the first scalar normal
%                   distribution.
%        mu2,Sigma2 The mean and variance of the second scalar normal
%                   distribution.
%
%OUTPUTS: mu, SigmaInv The mean and variance of the product distribution.
%
%Note that the product distribution is normalized. The derivation of the
%product distribution is straightforward nothing that the Fourier transform
%of a normal distribution is also a fully-normalized (integrates to one)
%normal distribution.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

mu=mu1+mu2;
Sigma=Sigma1+Sigma2;

end

function [momentVal,coeffPolyPart]=momentGenFun(mu,Sigma,numDerivs,t)
%%MOMENTGENFUN Evaluate the moment generating function (or one of its
%              derivatives) of the multivariate normal distribution. Taking
%              the ith, jth, kth... derivative of the moment generating
%              function with respect to the first, second, third...
%              components of the argument and evaluating it at t=0 provides
%              the noncentral moment of the multivariate normal
%              distribution involving the ith, jth, kth power of the first
%              second, third... components of the random vector.
%
%INPUTS: mu The mean of the PDF. If the PDF is multivariate, then this is a
%           numDimX1 column vector.
%     Sigma The variance (if scalar) or numDimXnumDim covariance matrix
%           (if multidimensional) of the PDF.
% numDerivs A numDimX1 or 1XnumDim vector indicating the number of
%           derivatives to take with respect to each of the components of
%           the argument of the moment generating function. numDerivs>=0.
%         t The numDimXnumPoints argument of the moment generating
%           function at which the derivatives of the moment generating
%           function should be evaluated. If this parameter is omitted or
%           an empty matrix is passed, then a default of
%           t=zeros(numDim,1) is used.
%
%OUTPUTS: momentVal A numPointsX1 vector of the values of the derivatives
%                   of the moment generating function given at the points
%                   in t or at t=0 if t is omitted.
%     coeffPolyPart A hypermatrix taking numDim indices that can be
%                   evaluated using the polyValMultiDim at different values
%                   of t to get the polynomial coefficient that multiplies
%                   the exponential term in the given set of derivatives
%                   of the moment generating function.
%
%The moment generating function of a random vector is defined to be
%E(exp(t'*x)) where E is the expected value operator, x is the random
%variable and t is a real parameter having the same dimensionality as the
%random variable. It can be shown that the moment generating function of a
%multivariate normal distribution is
%E(exp(t'*x))=exp(t'*mu+(1/2)*t'*Sigma*t)
%Derivatives of this can be evaluated systematically. First,
%differentiating the exponential with respect to the ith component of t
%leads to the original exponential term times
%mu(i)+sum(Sigma(:,i).*t(:))
%Thus, all derivatives include the original exponential term times a
%multivariate polynomial. This function keeps track of the polynomial
%through all of the derivatives, using the chain rule, the fact that
%multivariate polynomial multiplication can be performed using the convn
%function, and using polyDerMultiDim and polySumMultiDim for
%multidimensional polynomial differentiation and addition.
%
%As an example, consider finding the noncentral third moment E(x1*x2*x3) of
%a 3D multivariate normal distribution. This can be done using
%If one solves for it by hand, one gets
%mu(1)*Sigma(2,3)+mu(2)*Sigma(1,3)+mu(3)*Sigma(1,2)+prod(mu)
%However, the moment can also be found by evaluating the derivatives of the
%moment generating function with respect to the first, second and third
%variables at t=[0;0;0]. For example, consider
% mu=[1;2;3];
% Sigma=[8, 3, 2;
%        3,11, 9;
%        2, 9,18];
% numDerivs=[1;1;1];
% momentVal=GaussianD.momentGenFun(mu,Sigma,numDerivs)
%One will find the value of the noncentral moment to be 28, which is
%correct.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

%The number of dimensions.
numIdx=length(mu);

%Derivative of the moment generating function produce a polynomial times
%the original exponential term. We shall keep track of that polynomial
%while differentiating. The loops below take successive derivatives with
%respect to the different dimensions of the argument vector of the moment
%generating function. Each component is differentiated for the required
%number of derivatives. coeffPolyPart holds the accumulated polynomial that
%multiplies the exponential part of the differentiated moment generating
%function. Initially, coeffPolyPart=1 to signify that no derivatives have
%been taken. It is initially implemented with more elements than necessary
%to simplify the addition of the terms after the use of the polyDerMultiDim
%function below.
coeffPolyPart=1;

%i is the current dimension we are differentiating.
for i=1:numIdx
    %basePoly shall be the multivariate polynomial that is added every time
    %the exponential term in the moment generating function is
    %differentiated with respect to the ith index. This makes basePoly the
    %derivative with respect to the ith element of t of the argument of the
    %exponent in the multivariate normal moment generating function

    %Allocate space for the polynomial and make it have the correct shape.
    if(numIdx>1)
        basePoly=reshape(zeros(2^numIdx,1),2*ones(1,numIdx));
    else
        basePoly=zeros(2^numIdx,1);
    end

    %The additive term is the component of the mean that was multiplied by
    basePoly(1)=mu(i);

    %The term multiplied by x consists of elements from the covariance
    %matrix. This just consists of terms in the ith row of the
    %covariance matrix (the matrix is symmetric).
    idxVec=ones(numIdx,1);
    %The dimensionalities of all of the variables for the nDim2Index
    %function.
    dims=2*ones(numIdx,1);
    for curDim=1:numIdx
        idxVec(curDim)=2;
        basePoly(nDim2Index(dims,idxVec))=Sigma(i,curDim);
        idxVec(curDim)=1;
    end

    %Now, we enter into a loop to evaluate derivatives of the moment
    %generating function with respect to the current dimension.
    for derivsLeft=numDerivs(i):-1:1
        %The chain rule means that there are two terms to consider (both of
        %which are multiplied by the same exponential term).
        %The first term is the current coeffPolyPart times the derivative
        %of the exponential term with respect to the ith dimensions. This
        %is the product of coeffPolyPart and basePoly. The second term is
        %the exponential term times the derivative of coeffPolyPart. The
        %two terms then must be added.
        term1=convn(basePoly,coeffPolyPart);%Multiply the polynomials.
        term2=polyDerMultiDim(coeffPolyPart,i);

        %Add the multivariate polynomials.
        coeffPolyPart=polySumMultiDim(term1,term2);
    end
end

%coeffPolyPart now contains the multivariate polynomial that is multiplied
%by the exponential term. If no value of t is given, then just evaluate it
%at t=0. In this instance, the exponential term is zero and only the
%constant term from the polynomial appears.
if(nargin<4||isempty(t))
    momentVal=coeffPolyPart(1);
else
    numPoints=size(t,2);
    momentVal=zeros(numPoints,1);
    for curPoint=1:numPoints
        tCur=t(:,curPoint);
        momentVal(curPoint)=exp(tCur'*mu+tCur'*Sigma*tCur)*polyValMultiDim(coeffPolyPart,tCur);
    end
end

end

function cumVal=cumGenFun(mu,Sigma,numDerivs,t)
%%CUMGENFUN Evaluate the cumulant generating function (or one of its
%           derivatives) of the multivariate normal distribution. Taking
%           the ith, jth, kth... derivative of the cumulant generating
%           function with respect to the first, second, third...
%           components of the argument and evaluating it at t=0 provides
%           the cumulant of the multivariate normal distribution involving
%           the ith, jth, kth power of the first second, third...
%           components of the random vector. The cumulant generating
%           function is the natural logarithm of the moment generating
%           function.
%
%INPUTS: mu The mean of the PDF. If the PDF is multivariate, then this
%           is a numDimX1 column vector.
%     Sigma The variance (if scalar) or numDimXnumDim covariance matrix
%           (if multidimensional) of the PDF.
% numDerivs A numDimX1 or 1XnumDim vector indicating the number of
%           derivatives to take with respect to each of the components of
%           the argument of the cumulant generating function.
%           numDerivs>=0.
%         t The numDimXnumPoints argument of the cumulant generating
%           function at which the derivatives of the cumulant generating
%           function should be evaluated. If this parameter is omitted or
%           an empty matrix is passed, then a default of
%           t=zeros(numDim,1) is used.
%
%OUTPUTS: cumVal A numPointsX1 vector of the values of the derivatives
%                of the cumulant generating function given at the points
%                in t or at t=0 if t is omitted.
%
%Cumulants are useful in interpolating probability distributions through
%the use of, for example, an Edgeworth series. The cumulant generating
%function is defined as the natural logarithm of the moment generating
%function. It can be shown that the moment generating function of the
%multivariate normal distribution is
%E(exp(t'*x))=exp(t'*mu+(1/2)*t'*Sigma*t)
%Thus, the cumulant generating function is just
%t'*mu+(1/2)*t'*Sigma*t
%Consequently, for derivatives higher than two, the cumulant generating
%function (and hence the cumulants) of the multivariate normal distribution
%are all zero.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

if(nargin<4||isempty(t))
    numDim=length(mu);
    t=zeros(numDim,1);
end

numPoints=size(t,2);
cumVal=zeros(numPoints,1);

sumVal=sum(numDerivs);
switch(sumVal)
    case 0%No derivatives.
        for curPoint=1:numPoints
            tCur=t(:,curPoint);
            cumVal(curPoint)=tCur'*mu+(1/2)*tCur'*Sigma*tCur;
        end
    case 1
        %Find the nonzero term.
        derivIdx=find(numDerivs);
        for curPoint=1:numPoints
            tCur=t(:,curPoint);
            cumVal(curPoint)=mu(derivIdx)+sum(Sigma(:,derivIdx).*tCur(:));
        end
    case 2
        derivIdx=find(numDerivs);
        cumVal(:)=Sigma(derivIdx(1),derivIdx(2));%The same for all t.
    otherwise%No third or higher order cumulants.
       cumVal(:)=0;
end
end

function x=rand(N,mu,P)
%%RAND Generate multivariate Gaussian random variables with a given mean
%      vector and covariance matrix.
%
%INPUTS: N The number of random variables to generate.
%       mu The xDim X1 mean of the multivariate Gaussian to generate.
%        P The xDim X xDim positive definite covariance matrix of the
%          multivariate Gaussian to generate. If this parameter is omitted
%          or an empty matrix is passed, then the identity matrix will be
%          used.
%
%OUTPUTS: x An xDimXN matrix of random instances of the multivariate
%           Gaussian distribution.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    xDim=size(mu,1);
    if(nargin<3||isempty(P))
        P=zeros(xDim,xDim);
    end

    xDim=size(mu,1);
    x=bsxfun(@plus,mu,chol(P,'lower')*randn(xDim,N));
end

function x=randS(N,mu,S)
%%RANDS Generate multivariate Gaussian random variables with a given mean
%       vector and lower-triangular square root covariance matrix.
%
%INPUTS: N The number of random variables to generate.
%       mu The xDim X1 mean of the multivariate Gaussian to generate.
%        S The xDim X xDim lower triangular square root covariance matrix
%          of the multivariate Gaussian to generate.
%
%OUTPUTS: x An xDimXN matrix of random instances of the multivariate
%           Gaussian distribution.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    xDim=size(mu,1);
    x=bsxfun(@plus,mu,S*randn(xDim,N));
end

function [P,error,curIter]=integralOverRegion(mu, Sigma,minVals,maxVals,epsVal,alpha,maxIter)
%%INTEGRALOPVERREGION Compute the probability of a Gaussian probability
%                  density function (PDF) within a (hyper-)rectangular
%                  region. The probability is computed using a transformed
%                  Monte Carlo method designed for this specific integral
%                  that converges significantly faster than generic,
%                  textbook Monte Carlo integration techniques.
%
%INPUTS:  mu The NX1 mean vector of the multivariate normal distribution.
%      Sigma The NXN positive definite covariance matrix of the
%            distribution.
%    minVals An NX1 or 1XN vector of the lower integration bounds for all
%            of the dimensions.
%    maxVals An NX1 or 1XN vector of the upper integration bounds. Note
%            that minVals(i)<maxVals(i) for all elements.
%     epsVal The desired error tolerance for the probability computation.
%            If omitted or an empty matrix is passed, the defaut of 1e-6
%            (1e-4%) is used.
%      alpha The confidence factor for the standard error tolerance. Thus,
%            for the estimate of epsVal used to determine termination to be
%            correct 99% of the time, one should use
%            alpha=GaussianD.invCDF(0.99); If this parameter is omitted or
%            an empty matrix is passed, the default value of 3.5 is used.
%    maxIter The maximum number of iterations to use. If this parameter is
%            omitted or an empty matrix is passed, the default value of
%            1000 is used.
%
%OUTPUTS: P The approximate probability within the region (0-1).
%     error The estimated error in P, with a confidence interval determined
%           by alpha.
%   curIter The number of the last iteration performed before termination.
%
%The algorithm is that of [1]. Various transformations are applied to the
%function to achieve better convergence than when performing typical Monte
%Carlo integration.
%
%EXAMPLE:
%This is the 3D example in the paper. The probability should be about
%0.8279.
% minBounds=[-Inf;-Inf;-Inf];
% maxBounds=[1;4;2];
% R=[1,3/5,1/3;
%    3/5,1,11/15;
%    1/3,11/15,1];
% mu=[0;0;0];
% PApprox=GaussianD.integralOverRegion(mu,R,minBounds,maxBounds)
%
%REFERENCES:
%[1] A.Genz, "Numerical computation of multivariate normal probabilities,"
%    Journal of Computational and Graphical Statistics, vol. 1, no. 2, pp.
%    141-149, June 1992.
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

    if(nargin<7||isempty(maxIter))
       maxIter=1000;
    end

    if(nargin<6||isempty(alpha))
       alpha=3.5;
    end

    if(nargin<5||isempty(epsVal))
       epsVal=1e-6;
    end

    %Recenter things for a zero-mean distribution.
    minVals=minVals-mu;
    maxVals=maxVals-mu;

    %First, we reorder the variables so that the largest integration
    %regions come last. This is suggested in [1] to make the algorithm
    %faster.
    [~,idxSort]=sort(maxVals-minVals,'ascend');
    a=minVals(idxSort);
    b=maxVals(idxSort);
    C2=Sigma(idxSort,idxSort);

    C=chol(C2,'lower');

    numDim=size(C,1);

    intSum=0;
    varSum=0;

    %Allocate space
    d=zeros(numDim,1);
    e=zeros(numDim,1);
    f=zeros(numDim,1);

    d(1)=GaussianD.CDF(a(1)/C(1,1));
    e(1)=GaussianD.CDF(b(1)/C(1,1));
    f(1)=e(1)-d(1);

    y=zeros(numDim-1,1);%Allocate space.
    for curIter=0:(maxIter-1)

        w=rand(numDim-1,1);

        for i=2:numDim
            y(i-1)=GaussianD.invCDF(d(i-1)+w(i-1)*(e(i-1)-d(i-1)));
            deltaVal=C(i,1:(i-1))*y(1:(i-1));

            if(a(i)==-Inf)
                d(i)=0;
            else
                d(i)=GaussianD.CDF((a(i)-deltaVal)/C(i,i));
            end

            if(b(i)==Inf)
                e(i)=1;
            else
                e(i)=GaussianD.CDF((b(i)-deltaVal)/C(i,i));
            end
            f(i)=(e(i)-d(i))*f(i-1);
        end

        N=curIter+1;

        delta=(f(numDim)-intSum)/N;
        intSum=intSum+delta;
        varSum=(N-2)*varSum/N+delta^2;
        error=alpha*sqrt(varSum);

        if(error<epsVal)
            break;
        end
    end

    %The probability estimate.
    P=min(1,max(0,intSum));
end

function entropyVal=entropy(Sigma)
%%ENTROPY Obtain the differential entropy of the multivariate Gaussian
%         distribution given in nats. The differential entropy of a
%         continuous distribution is entropy=-int_x p(x)*log(p(x)) dx where
%         the integral is over all values of x. Units of nats mean that the
%         natural logarithm is used in the definition. Unlike the Shannon
%         entropy for discrete variables, the differential entropy of
%         continuous variables can be both positive and negative.
%
%INPUTS: Sigma The NXN positive definite covariance matrix of the
%              distribution.
%
%OUTPUTS: entropyVal The value of the differential entropy in nats.
%
%Differential entropy is defined in Chapter 8 of [1].
%
%REFERENCES:
%[1] T. M. Cover and J. A. Thomas, Elements of Information Theory, 2nd ed.
%    Hoboken, NJ: Wiley-Interscience, 2006.
%
%April 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.

    entropyVal=(1/2)*log(det(2*pi*exp(1)*Sigma));
end

end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [MW] = TW2MW(TW)
% Convert power from terawatts to megawatts.
% Chad A. Greene 2012
MW = TW*1000000;

author: Sable

repo: mcbench-benchmarks


text:
function []=plot_reprojectVSreal_points(CBimagesInfo,cameraCBparameters)
%% function for plot each image with reprojected points vs. real points and straight lines, and reprojection error statistics in STEP0
%
% INPUTS:
% * CBimagesInfo
% * cameraCBparameters:  a structure containing all the calibration parameters created in STEP0
%
%%
imageIndecesToDisplay=1:numel(cameraCBparameters.imagesUsed);

% plot each image with reprojected points vs. real points and straight lines
I=CBimagesInfo.I;
Nimages=numel(imageIndecesToDisplay);
imagePoints=cameraCBparameters.imagePoints(:,:,imageIndecesToDisplay);
reprojectedPoints=cameraCBparameters.cameraParameters.ReprojectedPoints(:,:,imageIndecesToDisplay);
reprojectionErrors=cameraCBparameters.cameraParameters.ReprojectionErrors(:,:,imageIndecesToDisplay);
Nrows=cameraCBparameters.boardSize(1);
Ncols=cameraCBparameters.boardSize(2);
icam=cameraCBparameters.icam;

reprojectedPointsMat=reshape(reprojectedPoints,Nrows-1,Ncols-1,2,Nimages);
reprojectedPointsHorizonal=reprojectedPointsMat([1 size(reprojectedPointsMat,1)],:,:,:);
reprojectedPointsVertical=reprojectedPointsMat(:,[1 size(reprojectedPointsMat,2)],:,:);

% PLOT
% define tabs
f = figure('name',['Reprojected points on all images taken by camera ' num2str(icam) '. Scroll between tabs to view the different images'],'units','normalized','outerposition',[.1 .1 .8 .8]);;
tabgp = uitabgroup(f);

icount=0;
for iplot=cameraCBparameters.imagesUsed(imageIndecesToDisplay)
    icount=icount+1;
    %plot
    tab(icount) = uitab(tabgp,'Title',['IM' num2str(iplot)]);

    axes('Parent',tab(icount)); % somewhere to plot

    subplot(1,6,1:5)

    imshow(I(:,:,:,iplot)); hold all;
    plot(imagePoints(:,1,icount), imagePoints(:,2,icount),'go','linewidth',1.5);
    plot(reprojectedPoints(:,1,icount),reprojectedPoints(:,2,icount),'r+','linewidth',1.5);
    title(['Camera ' num2str(icam) ' Image ' num2str(iplot)]);
    drawnow
    % plot straight lines
    plot(squeeze(reprojectedPointsHorizonal(:,1,1,icount)),squeeze(reprojectedPointsHorizonal(:,1,2,icount)),'-c','linewidth',1.5);
    plot(squeeze(reprojectedPointsVertical(1,:,1,icount)),squeeze(reprojectedPointsVertical(1,:,2,icount)),'-m','linewidth',1.5);
    for icol=2:Ncols-1
        plot(squeeze(reprojectedPointsHorizonal(:,icol,1,icount)),squeeze(reprojectedPointsHorizonal(:,icol,2,icount)),'-c','linewidth',1.5);
    end
    for irow=1:Nrows-1
        plot(squeeze(reprojectedPointsVertical(irow,:,1,icount)),squeeze(reprojectedPointsVertical(irow,:,2,icount)),'-m','linewidth',1.5);
    end
    legend('Detected Points','Reprojected Points','straight horizontal lines','straight vertical lines');
    hold off;

    reprojectionErrorsNow=reprojectionErrors(:,:,icount);
    reprojectionErrorsMgnNow=sqrt(sum(reprojectionErrorsNow.^2,2));
    reprojectionErrorsNow=[reprojectionErrorsNow reprojectionErrorsMgnNow];

    subplot(1,6,6)
    boxplot(reprojectionErrorsNow,'Labels',{'X','Y','Mgn'});
    ylim([-max(abs(reprojectionErrorsNow(:))) max(abs(reprojectionErrorsNow(:)))]);
    title({'Reprojection error'; 'statistics [pix]'});

    drawnow
end

end

%%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function [ eigenvalues, kpoints, nelect ] = import_eigenval( filename )
%IMPORT_EIGENVAL Import a VASP EIGENVAL file.
%   [eigenvalues,kpoints,nelect] = IMPORT_EIGENVAL(filename) imports a VASP
%   EIGENVAL file. Optional parameter filename specifies the name of the
%   file. eigenvalues is a NKPOINTS x NBANDS x ISPIN array of eigenvalues,
%   kpoints is a NKPOINTS x 4 array of k-point coordinates and weights, and
%   nelect is the number of electrons.

% todo:
% check compatibility with non-spin-polarized files

  if nargin == 0
      filename='EIGENVAL';
  end

  fid = fopen(filename);
  if fid==-1
    error(['File ' filename ' not found']);
  end

  buffer = fscanf(fid, '%d', 4); % various data
  ispin = buffer(4); % 1 = non-polarized, 2 = polarized
  fgetl(fid); % empty string
  fgetl(fid); % various data
  fgetl(fid); % various data
  fgetl(fid); % various data
  fgetl(fid); % comment
  nelect = fscanf(fid, '%d', 1); % number of electrons
  nkpoints = fscanf(fid, '%d', 1); % number of k-points
  nbands = fscanf(fid, '%d', 1); % number of bands

  fgetl(fid); % empty string

  kpoints = zeros(nkpoints,4);
  eigenvalues = zeros(nkpoints, nbands, ispin);

  for kpoint = 1:nkpoints
    fgetl(fid); % blank line
    kpoints(kpoint,:) = fscanf(fid, '%f', 4)';
    fgetl(fid); % empty string
    for band = 1:nbands
        buffer = fscanf(fid,'%f',ispin+1);
        eigenvalues(kpoint,band,1:ispin) = buffer(2:(ispin+1));
        fgetl(fid); % empty string
    end
  end

  eigenvalues = sort(eigenvalues,2); % sort the bands by energy

end


text:
function [fds, ics, averFD, averIC] = fdvolfft(im)
% FDVOLFFT Compute fractal dimensions FDS and intercepts ICS of a 3D image
%          IM using FFT and draw rose plots of FDS and ICS
% IM: input 3D array of a 3D image (grayvalue)
% FDS: a 2D array of size 24x12 which stores the FDs in 24x12 directions
% ICS: a 2D array of size 24x12 which stores the average intercepts in
%      24x12 directions
% AVERFD: average fractal dimension for all directions
% AVERIC: average intercept for all directions
%
%
% This script is free to use except for commercial purpose
% Written by  Mr. Jianbo ZHANG, J.Zhang@ewi.utwente.nl
% 19 Jan, 2005
%
%
% for example, the following script can compute the fractal dimension of
%    an MRI image volume
%
% load mri
% D = squeeze(D) ;
% image_num = 8;
% image(D(:,:,image_num))
% axis image
% colormap(map)
% title('Image of MRI slice No. 8 ')
% [fds ics averFD averIC]= fdvolfft(D)
%
% Note: Since an MRI image is a self-affine, not self-similiar volume, its
%   fractal dimensions do not necessarily lie between 3 and 4.

tic
NUM_AZI = 24; % number of azimuth directions
              % that the frequency space is evenly divided
NUM_ZEN = 12; % number of zenith directions
              % that the frequency space is evenly divided
NUM_RAD = 30; % number of points that the radial line is evenly divided

if nargin < 1,
    error('Missed input argument which must be a 3D array!')
end

[M N P] = size(im);
xctr = 1 + bitshift(N, -1); % x coordinate of center point
yctr = 1 + bitshift(M, -1); % y coordinate of center point
zctr = 1 + bitshift(P, -1); % z coordinate of center point
imMean = mean(im(:));
fim = fftshift(fftn(double(im) - imMean));
psd = log(fim .* conj(fim) + 10^(-6));
sumBrite = zeros(NUM_AZI, NUM_ZEN, NUM_RAD); %accumulation PSD
                                          %for along each direction and
                                          %radius
nCount = zeros(NUM_AZI, NUM_ZEN, NUM_RAD); % PSD count
radius = zeros(2 * NUM_RAD,1); %accumulation PSD for all directions
radCount = zeros(2 * NUM_RAD,1); % number of PSD for all directions

%Compute phase distribution
phase = zeros(180);
for k = 1:P
    for j = 1:M
        for i = 1:N
            realv = real(fim(j,i,k));
            imagv = imag(fim(j,i,k));
            if realv == 0
                value = pi/2;
            else
                value = atan((imagv / realv));
                ang = floor(180 * (pi / 2 + value) / pi);
            end
            if ang < 0
                ang = 0;
            end
            if ang > 179
                ang = 179;
            end
            phase(ang + 1) = phase(ang + 1) + 1;
        end
    end
end
maxphase = max(phase);
figure, plot(phase / maxphase);
title('Phase histogram (0...2 \pi)');%
axis off

%accumulation of PSD for each direction and radius
rmax = log(min(min(M,N),P)/2);% maximum radius
for k = 1:P
    if k ~= zctr
        zval = zctr - k;
        z2 = zval * zval;
        for j = 1:M
            if j ~= yctr
                yval = yctr - j;
                y2 = yval * yval;
                for i = 1:N
                    if i ~= xctr
                        xval = i - xctr;
                        r = sqrt(z2 + y2 + xval * xval);
                        rho = log(r);
                        if rho > 0 & rho <= rmax
                            mval = psd(j,i,k);
                            temp1 = yval / xval;
                            temp2 = zval / r;
                            theta = atan(temp1);
                            phi = acos(temp2);
                            if xval < 0
                                theta = theta + pi;
                            end
                            if theta < 0
                                theta = theta + 2 * pi;
                            end

                            aziSN = floor(NUM_AZI * theta /(2 * pi));
                            if aziSN > NUM_AZI - 1 | aziSN < 0
                                aziSN = NUM_AZI - 1;
                            end
                            zenSN = floor(NUM_ZEN * phi / pi);
                            if zenSN > NUM_ZEN - 1 | zenSN < 0
                                zenSN = NUM_ZEN - 1;
                            end
                            radSN = floor(2 * NUM_RAD * rho / rmax);
                            h = floor(radSN/2);
                            if radSN > 2 * NUM_RAD - 1
                                h = NUM_RAD - 1;
                                radSN = 2 * NUM_RAD - 1;
                            end

                            if h >= 5
                                if zenSN == 0 | zenSN == NUM_ZEN-1
                                    sumBrite(:, zenSN + 1, h + 1) = sumBrite(:, zenSN + 1, h + 1) + mval;
                                    nCount(:, zenSN + 1, h + 1) = nCount(:,zenSN + 1, h + 1) + 1;
                                else
                                    sumBrite(aziSN + 1,zenSN + 1, h + 1) = sumBrite(aziSN + 1, zenSN + 1, h + 1) + mval;
                                    nCount(aziSN + 1, zenSN + 1, h + 1) = nCount(aziSN + 1,zenSN + 1, h + 1) + 1;
                                end
                            end
                            if radSN >= 5
                                radius(radSN + 1) = radius(radSN + 1) + mval;
                                radCount(radSN + 1) = radCount(radSN + 1) + 1;
                            end
                        end
                    end
                end
            end
        end
    end
end

%linear regression
for aziSN = 1:NUM_AZI
    for zenSN = 1: NUM_ZEN
        sumx = 0;
        sumy = 0;
        sumx2 = 0;
        sumxy = 0;
        sumn = 0;
        for range = 6:NUM_RAD
            if nCount(aziSN, zenSN, range) > 0
                yval = sumBrite(aziSN, zenSN, range)/nCount(aziSN, zenSN, range);
                xval = (range -1) * rmax / NUM_RAD;
                sumx = sumx + xval;
                sumy = sumy + yval;
                sumx2 = sumx2 + xval * xval;
                sumxy = sumxy + xval * yval;
                sumn = sumn + 1;
            end
        end
        slope(aziSN, zenSN) = (sumn * sumxy - sumx * sumy) / (sumn * sumx2 - sumx * sumx);
        ics(aziSN, zenSN) = (sumy - slope(aziSN, zenSN) * sumx) / sumn;
        fds(aziSN, zenSN) = (11 + slope(aziSN, zenSN))/2;
    end
end

%compute average slope over all directions and radius
sumn = 0;
for radSN = 6:(2 * NUM_RAD)
    if radCount(radSN) > 0
        sumn = sumn + 1;
        yval(sumn) = radius(radSN) / radCount(radSN);
        tempr(sumn) = (radSN -1) * rmax / (2 * NUM_RAD);
    end
end
p = polyfit(tempr,yval,1);
averFD =(11+ p(1))/2; % Overall average fractal dimension
averIC = p(2); % Overall average intercept
fitln = polyval(p,tempr);
figure; plot(tempr,yval,tempr,fitln,'r-');
title('Log Log plot of PSD. vs Freq.');
ylabel('Log PSD');
xlabel('Log Frequency');
legend('data','LSE fit line');

% draw hedgehog plot of fractal dimension and intercept
% I personally prefer to call these 3D visulization plots of fractal
% dimension and intercept as 'hedgehog' plot
fd = [fds; fds(1,:)];
fd = [fd, fd(:,1)];
aziSN = (1:NUM_AZI + 1)';
zenSN = 1:(NUM_ZEN + 1);
[ph th] = meshgrid(pi/2 - ((zenSN - 1) * pi / NUM_ZEN) , 2 * pi / NUM_AZI / 2 + (aziSN - 1) * 2 * pi / NUM_AZI);
[x,y,z] = sph2cart(th,ph, fd);
figure, surf(x,y,z,fd); colorbar
title('Hedgehog plot of fractal dimension');

ic = [ics; ics(1,:)];
ic = abs([ic, ic(:,1)]); % ic stores the absolue value of ics
aziSN = (1:NUM_AZI + 1)';
zenSN = 1:(NUM_ZEN + 1);
[ph th] = meshgrid(pi/2 - ((zenSN - 1) * pi / NUM_ZEN) , 2 * pi / NUM_AZI / 2 + (aziSN - 1) * 2 * pi / NUM_AZI);
[x,y,z] = sph2cart(th,ph, ic);
figure, surf(x,y,z, ic); colorbar
title('Hedgehog plot of intercept');

disp(['Elapsed time: ' num2str(toc)]);


text:
% function hs = tplot_power(Tx, t, fs, opt)
% function hs = tplot_power(Wx, t, as, opt)
%
% Plots the 2D magnitude of either the synchrosqueezing or
% wavelet transform of signal x via tplot.
% Also plots the median power estimate for each frequency/scale on
% the side, scaled to the magnitude of the original signal.
%
% Input (see help tplot):
%   Tx, t, fs (e.g. output of synsq_cwt_fw, or properly chosen slices)
%     Use opt.style = 'freq'.
%  or
%   Wx, t, as (e.g., output of cwt_fw, or properly chosen slices)
%     Use opt.style = 'scale'.
%
%  Additional options to tplot's opt structure:
%   opt.filter: remove/filter side lobes before estimating
%               power? (values: 1 or 0, default: 1)
%   opt.type: wavelet type (necessary to scale power plots
%             accurately, see help wfiltfn, help cwt_fw)
%
% Output:
%   hs - 2-dim vector of handles to the tplot, and the power plot
%    hs(1) - handle to the tplot axis
%    hs(2) - handle to the 90 degree angled power plot axis
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------
function hs = tplot_power(Tx, t, fs, opt)
  if nargin<4, opt = struct(); end

  if ~isfield(opt, 'type'),
    warning(['tplot_power: opt.type not known, db magnitudes only ' ...
             'accurate up to scale']);
    Css = 1;
  else
    Css = synsq_adm(opt.type);
  end

  [na0, n] = size(Tx);
  hs = zeros(1,2);
  hs(1) = subplot(1,2,1);
  [tmp,opt] = tplot(Tx,t,fs,opt);
  xlabel('t');
  grid on;

  % Find frequency axis limits
  if fs(2)>fs(1) % Standard
      flim = [max(opt.flim(1), fs(1)), min(opt.flim(2), fs(end))];
      flimi(1) = find(fs >= flim(1), 1, 'first');
      flimi(2) = find(fs <= flim(2), 1, 'last');
  else % Periodicity
      flim = [max(opt.flim(1), fs(end)), min(opt.flim(2), fs(1))];
      flimi(1) = find(fs >= flim(2), 1, 'last');
      flimi(2) = find(fs <= flim(1), 1, 'first');
  end

  % Restrict ourselves to this region
  aTx = abs(Tx); clear Tx;
  fs = fs(flimi(1):flimi(2));
  aTx = aTx(flimi(1):flimi(2), :);

  % Calculate the median of the absolute values, filtered properly.
  [muTx,Lbdi,Rbdi] = synsq_filtered_time_quantile(aTx, t, fs, opt, .5);

  mlog2 = @(x) x;
  if opt.clog
      mlog2 = @(x) log2(x);
  end

  if (~opt.bd)
      hold on;
      plot(t(Lbdi), mlog2(fs), '--k');
      plot(t(Rbdi), mlog2(fs), '--k');
  end

  % Mean of aTx plot
  hs(2) = subplot(1,2,2);

  plot(mlog2(fs), 10*log10(1/Css*muTx));
  axis tight;
  xlim([min(mlog2(fs)), max(mlog2(fs))]);
  %  %ylim([min(10*log10(muTx)), max(10*log10(muTx))]);
  set(gca, 'XTick', cellfun(@(x)mlog2(str2num(x)), opt.ticklabels));
  set(gca, 'XTickLabel', '');
  xlabel('');
  grid on;

  if (fs(2)>fs(1))
      view([90 -90]);
  else
      view([90 90]);
  end
  ylabel('\mu(| . |) [db]');

  set(hs(1), 'Units', 'normalized', 'OuterPosition', [0 0 .75 1.05])
  set(hs(2), 'Units', 'normalized', 'OuterPosition', [.75 0 .25 1.05])
end

author: HeLiangHIT


text:
function calpak_test695 ( )

%*****************************************************************************80
%
%% CALPAK_TEST695 tests YMDF_DIF_ENGLISH.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 April 2013
%
%  Author:
%
%    John Burkardt
%
  dhi = 1;
  dlo = 1;
  fhi = 0.0;
  flo = 0.0;
  mhi = 1;
  mlo = 1;
  seed = 123456789;
  yhi = 1960;
  ylo = 1970;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST695\n' );
  fprintf ( 1, '  YMDF_DIF_ENGLISH gets the day difference\n' );
  fprintf ( 1, '  between YMDF dates.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  YMDF1        YMDF2        (YMDF2 - YMDF1)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10

    [ y1, m1, d1, f1, seed ] = ymdf_uniform_english ( ylo, mlo, dlo, flo, ...
      yhi, mhi, dhi, fhi, seed );

    s1 = ymdf_to_s_english ( y1, m1, d1, f1 );

    [ y2, m2, d2, f2, seed ] = ymdf_uniform_english ( ylo, mlo, dlo, flo, ...
      yhi, mhi, dhi, fhi, seed );

    s2 = ymdf_to_s_english ( y2, m2, d2, f2 );

    [ days, ierror ] = ymdf_dif_english ( y1, m1, d1, f1, y2, m2, d2, f2 );

    fprintf ( 1, '  %10s  %10s  %11.2f\n', s1, s2, days );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function svn(varargin)
%svn Control TortoiseSVN from within matlab.
% Wrapper for TortoiseSVN.
% Useage: SVN command dir/file(s)
% Supported commands: about, checkout*, import*, update*, commit*, add,
% revert*, cleanup*, resolve*, repocreate*, switch*, export*, mergeall*,
% settings, remove, rename, diff, conflicteditor, relocate, help,
% repostatus, repobrowser, ignore, blame, createpatch, revisiongraph, lock,
% unlock, properties.
%
% Commands with an * only take one directory as an argument.
% Other commands that take a path can accept more than one argument. If no
% directory is given, the current working directory is used.
%
% This function is not intended to completely replace SVN. It is a quick
% way to do commits, checkouts and adds to a directory while in Matlab.
%
% Examples:
%   svn checkout [Dialog]
%   edit newfile.m [Dialog]
%   svn add newfile.m [Dialog]
%   svn commit [Dialog]
%
% For a full description of what each command does, visit the tortoisesvn
% website: http://tortoisesvn.net/docs/nightly/TortoiseSVN_en/tsvn-automation.html
%
% See also: add, checkout, commit
if nargin<1
    error(sprintf('You must give at least one command\nSee: http://tortoisesvn.net/docs/nightly/TortoiseSVN_en/tsvn-automation.html'));
else
    command=varargin{1};
end
if nargin<2
    varargin{2}=pwd;
end
switch command
    case 'about'
    case 'checkout' % Check out one path at a time
        varargin=singleDirOnly(command,varargin);
    case 'import'% Import one path at a time
        varargin=singleDirOnly(command,varargin);
    case 'update'
        varargin=singleDirOnly(command,varargin);
    case 'commit' % Commit path at a time
        varargin=singleDirOnly(command,varargin);
    case 'add'
    case 'revert'
        varargin=singleDirOnly(command,varargin);
    case 'cleanup'
        varargin=singleDirOnly(command,varargin);
    case 'resolve'
        varargin=singleDirOnly(command,varargin);
    case 'repocreate'
        varargin=singleDirOnly(command,varargin);
    case 'switch'
        varargin=singleDirOnly(command,varargin);
    case 'export'
        varargin=singleDirOnly(command,varargin);
    case 'merge'
        error('Unsupported command.');
    case 'mergeall'
        varargin=singleDirOnly(command,varargin);
    case 'copy'
        error('Unsupported command.');
    case 'settings'
    case 'remove'
    case 'rename'
    case 'diff'
    case 'showcompare'
        error('Unsupported command.');
    case 'conflicteditor'
    case 'relocate'
    case 'help'
    case 'repostatus'
    case 'repobrowser'
    case 'ignore'
    case 'blame'
    case 'cat'
        error('Unsupported command.');
    case 'createpatch'
    case 'revisiongraph'
    case 'lock'
    case 'unlock'
    case 'properties'
    otherwise
        error(sprintf('Unknown TortoiseSVN command: %s\nSee: %s',command,'http://tortoisesvn.net/docs/nightly/TortoiseSVN_en/tsvn-automation.html'));
end
% Use double quotes if there is a space in the path.
% No UNC as this is called using the dos command
tsvn='"C:\Program Files\TortoiseSVN\bin\TortoiseProc.exe"';
for i=2:numel(varargin)
    cmd=sprintf('%s /command:%s /closeonend:1 /path:"%s"',tsvn,command,abspath(varargin{i}));
    dos(cmd);
end
end

% For functions which should only take a single directory
function in=singleDirOnly(cmd,in)
if (~isdir(abspath(in{2})))
    error('You can only commit directories');
end
if length(in)>2
    warning('SVN:DIRLIMIT',['You can only ' cmd ' one directory at a time. Only using first directory']);
    in(3:end)='';
end
end

function [absolutepath]=abspath(partialpath)
% Taken from xlswrite.
% parse partial path into path parts
[pathname filename ext] = fileparts(partialpath);
% no path qualification is present in partial path; assume parent is pwd, except
% when path string starts with '~' or is identical to '~'.
if isempty(pathname) && isempty(strmatch('~',partialpath))
    Directory = pwd;
elseif isempty(regexp(partialpath,'(.:|\\\\)','once')) && ...
        isempty(strmatch('/',partialpath)) && ...
        isempty(strmatch('~',partialpath));
    % path did not start with any of drive name, UNC path or '~'.
    Directory = [pwd,filesep,pathname];
else
    % path content present in partial path; assume relative to current directory,
    % or absolute.
    Directory = pathname;
end
% construct absulute filename
absolutepath = fullfile(Directory,[filename,ext]);
end


text:
function [seq, ground_truth] = load_video_info(video_path, video_name)

ground_truth = dlmread([video_path '/groundtruth_rect.txt']);
% ground_truth = dlmread([video_path, '/', video_name, '_gt', '.txt']);

seq.len = size(ground_truth, 1);
seq.init_rect = ground_truth(1,:);

img_path = [video_path '/img/'];

if exist([img_path num2str(1, '%04i.png')], 'file'),
    img_files = num2str((1:seq.len)', [img_path '%04i.png']);
elseif exist([img_path num2str(1, '%04i.jpg')], 'file'),
    img_files = num2str((1:seq.len)', [img_path '%04i.jpg']);
elseif exist([img_path num2str(1, '%04i.bmp')], 'file'),
    img_files = num2str((1:seq.len)', [img_path '%04i.bmp']);
else
    error('No image files to load.')
end

seq.s_frames = cellstr(img_files);

end

repo: UDT