text:
classdef prtRvVq < prtRv
    % prtRvVq  Vector quantization random variable
    %
    %   RV = prtRvVq creates a prtRvVq object with empty means and
    %   probabilities. The means and probabilties must be set either
    %   directly, or by calling the MLE method.
    %
    %   Vector quanitization uses k-means to discretize the data space
    %   using nCategories means. The means are the discrete points in space and
    %   have probabilies representing their prominence in the data. The
    %   pdf is calculated by mapping to the nearest entry of the means and
    %   giving the data point the corresponding entry in probabilities.
    %
    %   RV = prtRvVq(PROPERTY1, VALUE1,...) creates a prtRvVq object RV
    %   with properties as specified by PROPERTY/VALUE pairs.
    %
    %   A prtRvVq object inherits all properties from the prtRv class. In
    %   addition, it has the following properties:
    %
    %   nCategories   - The number of categories
    %   means         - The means of each catergory that are used to
    %                   approximate the density
    %   probabilities - The probabilities of each category
    %
    %  A prtRvVq object inherits all methods from the prtRv class. The MLE
    %  method can be used to estimate the distribution parameters from
    %  data.
    %
    %  Example:
    %
    %  dataSet = prtDataGenUnimodal;        % Load a dataset consisting of
    %                                       % 2 features
    %  dataSet = retainFeatures(dataSet,1); % Retain only the first feature
    %
    %  RV = prtRvVq;                        % Create a prtRvVq object
    %  RV = RV.mle(dataSet);                % Compute the VQ parameters
    %                                       % form the data
    %  RV.plotPdf                           % Plot the pdf
    %
    %   See also: prtRv, prtRvMvn, prtRvGmm, prtRvMultinomial,
    %   prtRvUniform, prtRvUniformImproper

    properties (SetAccess = private)
        name = 'Vector Quantization Random Variable'
        nameAbbreviation = 'RVVQ';
    end

    properties (SetAccess = protected)
        isSupervised = false;
        isCrossValidateValid = true;
    end

    properties (Dependent = true)
        probabilities   % The probabilities
        means           % The means
        nCategories     % The number of categories
    end

    properties (Dependent = true, Hidden=true)
        InternalKMeansPrototypes
        InternalDiscrete
        nDimensions
    end

    properties (SetAccess = 'private', GetAccess = 'private', Hidden=true)
        %InternalKMeansPrototypesDepHelp = prtClassKmeansPrototypes('kMeansHandleEmptyClusters','random');
        meanDepHelper = [];
        InternalDiscreteDepHelp = prtRvDiscrete();
        nCategoriesDepHelp = 2;
    end

    methods
        % The Constructor
        function R = prtRvVq(varargin)
            R = constructorInputParse(R,varargin{:});
        end
        function val = get.means(R)
            val = R.meanDepHelper;
        end
        function val = get.probabilities(R)
            val = R.InternalDiscreteDepHelp.probabilities;
        end
        function val = get.nCategories(R)
            val = R.nCategoriesDepHelp;
        end

        function val = get.InternalDiscrete(R)
            val = R.InternalDiscreteDepHelp;
        end

        function R = set.InternalDiscrete(R,val)
            R.InternalDiscreteDepHelp = val;
        end

        function R = set.InternalDiscreteDepHelp(R,val)
            assert(isa(val,'prtRvDiscrete'),'InternalDiscrete must be a prtRvDiscrete.')
            R.InternalDiscreteDepHelp = val;
        end

        function R = set.nCategories(R,val)
            assert(numel(val)==1 && val==floor(val) && val > 0,'nCategories must be a scalar positive integer.');

            R.nCategoriesDepHelp = val;
        end

        function R = set.probabilities(R,val)
            assert(isnumeric(val) && isvector(val),'probabilities must be numer vector whose values sum to one.');

            if ~isempty(R.InternalDiscreteDepHelp.symbols)
                assert(size(R.InternalDiscreteDepHelp.symbols,1) == numel(val),'size mismatch between probabilities and means.')
            end
            R.InternalDiscreteDepHelp.probabilities = val(:);
            R.nCategories = numel(val);
        end

        function R = set.means(R,val)
            assert(ndims(val)==2 && isnumeric(val),'means must be a 2D numeric matrix.')

            if ~isempty(R.InternalDiscreteDepHelp.probabilities)
                if  numel(R.InternalDiscreteDepHelp.probabilities) ~= size(val,1)
                    if isvector(val)
                        error('prt:prtRvVq','means must be an NxD matrix where D is the dimensionality and N is the number of components. You may need to transpose the supplied value.')
                    else
                        error('prt:prtRvVq','size mismatch between probabilities and means.')
                    end
                end
            end

            R.nCategories = size(val,1);

            R.InternalDiscrete.symbols = val;
            R.meanDepHelper = val;
        end

        function val = get.nDimensions(R)
            if ~isempty(R.means)
                val = size(R.means,2);
            else
                val = [];
            end
        end

        function R = mle(R,X)
            X = R.dataInputParse(X); % Basic error checking etc
            assert(isnumeric(X) && ndims(X)==2,'X must be a 2D numeric array.');

            R.means = prtUtilKmeans(X,R.nCategories,'handleEmptyClusters','random');

            trainingOutput = R.data2ClosestMeanInd(X);

            R.InternalDiscrete = prtRvDiscrete('symbols',R.means,'probabilities',histc(trainingOutput,1:R.nCategories)/size(X,1));
        end

        function vals = pdf(R,X)
            [isValid, reasonStr] = R.isValid;
            assert(isValid,'PDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            X = R.dataInputParse(X); % Basic error checking etc
            assert(size(X,2) == R.nDimensions,'Data, RV dimensionality missmatch. Input data, X, has dimensionality %d and this RV has dimensionality %d.', size(X,2), R.nDimensions)
            assert(isnumeric(X) && ndims(X)==2,'X must be a 2D numeric array.');

            trainingOutput = R.data2ClosestMeanInd(X);
            vals = R.probabilities(trainingOutput);
            vals = vals(:);
        end

        function vals = logPdf(R,X)
            [isValid, reasonStr] = R.isValid;
            assert(isValid,'LOGPDF cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            vals = log(pdf(R,X));
        end

        function vals = draw(R,N)
            if nargin < 2 || isempty(N)
                N = 1;
            end

            [isValid, reasonStr] = R.isValid;
            assert(isValid,'DRAW cannot yet be evaluated. This RV is not yet valid %s.',reasonStr);

            assert(numel(N)==1 && N==floor(N) && N > 0,'N must be a positive integer scalar.')

            vals = R.InternalDiscrete.draw(N);
        end

        function varargout = plotPdf(R,varargin)
            h = plotPdf(R.InternalDiscrete);

            varargout = {};
            if nargout
                varargout = {h};
            end
        end

        function varargout = plotCdf(R,varargin)
            h = plotCdf(R.InternalDiscrete);

            varargout = {};
            if nargout
                varargout = {h};
            end
        end

        function quantizedData = vq(R,X)
            X = R.dataInputParse(X); % Basic error checking etc
            quantizedData = data2ClosestMeanInd(R,X);
        end
    end

    methods (Hidden = true)
        function [val, reasonStr] = isValid(R)
            if numel(R) > 1
                val = false(size(R));
                for iR = 1:numel(R)
                    [val(iR), reasonStr] = isValid(R(iR));
                end
                return
            end

            [val, reasonStr] = isValid(R.InternalDiscrete);
            reasonStr = strrep(reasonStr,'symbols','means');
        end

        function val = plotLimits(R)
            val = plotLimits(R.InternalDiscrete);
        end

        function val = isPlottable(R)
            val = isPlottable(R.InternalDiscrete);
        end
    end

    methods (Hidden = true, Access=protected)
        function closestMeanInds = data2ClosestMeanInd(R,X)
            distance = prtDistanceEuclidean(X,R.means);
            [dontNeed, closestMeanInds] = min(distance,[],2); %#ok<ASGLU>
        end
    end

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function FigSize(xdim,ydim)
% =======================================================================
% Sets the figure window to a desired size. The default is full screen
% =======================================================================
% FigSize(xdim,ydim)
% -----------------------------------------------------------------------
% OPTIONAL INPUT
%   - xdim: horizontal size (in cm)
%   - ydim: vertical size (in cm)
% -----------------------------------------------------------------------
% EXAMPLE
%   FigSize(12,6)
%   plot(1:10)
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2015. Updated November 2020
% -----------------------------------------------------------------------

if nargin==0
    figure('units','normalized','outerposition',[0 0.1 1 0.9])
else
    % Get the position of the figure
    % [left bottom width height] [standard in px: 403 246   560   420]
    position = get(gcf, 'Position');

    % Place the figure on the south west corner of the screen (in cm)
    position(1) = 1;
    position(2) = 2;

    % xdim dimension in cm
    if exist('xdim','var')
        position(3) = xdim;
    end

    % ydim dimension in cm
    if exist ('ydim','var')
        position(4) = ydim;
    end

    % Apply rescaling
    set(gcf, 'units', 'centimeters', 'Position', position)
end

author: ambropo

repo: VAR-Toolbox

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


text:
function [net, info] = cnn_cifar(varargin)
% CNN_CIFAR   Demonstrates MatConvNet on CIFAR-10
%    The demo includes two standard model: LeNet and Network in
%    Network (NIN). Use the 'modelType' option to choose one.

run(fullfile(fileparts(mfilename('fullpath')), ...
  '..', '..', 'matlab', 'vl_setupnn.m')) ;

opts.modelType = 'lenet' ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.expDir = fullfile('data', sprintf('cifar-%s', opts.modelType)) ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.dataDir = fullfile('data','cifar') ;
opts.imdbPath = fullfile(opts.expDir, 'imdb.mat');
opts.whitenData = true ;
opts.contrastNormalization = true ;
opts.networkType = 'simplenn' ;
opts.train = struct() ;
opts = vl_argparse(opts, varargin) ;

% -------------------------------------------------------------------------
%                                                    Prepare model and data
% -------------------------------------------------------------------------

switch opts.modelType
  case 'lenet'
    net = cnn_cifar_init('networkType', opts.networkType) ;
  case 'nin'
    net = cnn_cifar_init_nin('networkType', opts.networkType) ;
  otherwise
    error('Unknown model type ''%s''.', opts.modelType) ;
end

if exist(opts.imdbPath, 'file')
  imdb = load(opts.imdbPath) ;
else
  imdb = getCifarImdb(opts) ;
  mkdir(opts.expDir) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

net.meta.classes.name = imdb.meta.classes(:)' ;

% -------------------------------------------------------------------------
%                                                                     Train
% -------------------------------------------------------------------------

switch opts.networkType
  case 'simplenn', trainfn = @cnn_train ;
  case 'dagnn', trainfn = @cnn_train_dag ;
end

[net, info] = trainfn(net, imdb, getBatch(opts), ...
  'expDir', opts.expDir, ...
  net.meta.trainOpts, ...
  opts.train, ...
  'val', find(imdb.images.set == 3)) ;

% -------------------------------------------------------------------------
function fn = getBatch(opts)
% -------------------------------------------------------------------------
switch lower(opts.networkType)
  case 'simplenn'
    fn = @(x,y) getSimpleNNBatch(x,y) ;
  case 'dagnn'
    bopts = struct('numGpus', numel(opts.train.gpus)) ;
    fn = @(x,y) getDagNNBatch(bopts,x,y) ;
end

% -------------------------------------------------------------------------
function [images, labels] = getSimpleNNBatch(imdb, batch)
% -------------------------------------------------------------------------
images = imdb.images.data(:,:,:,batch) ;
labels = imdb.images.labels(1,batch) ;
if rand > 0.5, images=fliplr(images) ; end

% -------------------------------------------------------------------------
function inputs = getDagNNBatch(opts, imdb, batch)
% -------------------------------------------------------------------------
images = imdb.images.data(:,:,:,batch) ;
labels = imdb.images.labels(1,batch) ;
if rand > 0.5, images=fliplr(images) ; end
if opts.numGpus > 0
  images = gpuArray(images) ;
end
inputs = {'input', images, 'label', labels} ;

% -------------------------------------------------------------------------
function imdb = getCifarImdb(opts)
% -------------------------------------------------------------------------
% Preapre the imdb structure, returns image data with mean image subtracted
unpackPath = fullfile(opts.dataDir, 'cifar-10-batches-mat');
files = [arrayfun(@(n) sprintf('data_batch_%d.mat', n), 1:5, 'UniformOutput', false) ...
  {'test_batch.mat'}];
files = cellfun(@(fn) fullfile(unpackPath, fn), files, 'UniformOutput', false);
file_set = uint8([ones(1, 5), 3]);

if any(cellfun(@(fn) ~exist(fn, 'file'), files))
  url = 'http://www.cs.toronto.edu/~kriz/cifar-10-matlab.tar.gz' ;
  fprintf('downloading %s\n', url) ;
  untar(url, opts.dataDir) ;
end

data = cell(1, numel(files));
labels = cell(1, numel(files));
sets = cell(1, numel(files));
for fi = 1:numel(files)
  fd = load(files{fi}) ;
  data{fi} = permute(reshape(fd.data',32,32,3,[]),[2 1 3 4]) ;
  labels{fi} = fd.labels' + 1; % Index from 1
  sets{fi} = repmat(file_set(fi), size(labels{fi}));
end

set = cat(2, sets{:});
data = single(cat(4, data{:}));

% remove mean in any case
dataMean = mean(data(:,:,:,set == 1), 4);
data = bsxfun(@minus, data, dataMean);

% normalize by image mean and std as suggested in `An Analysis of
% Single-Layer Networks in Unsupervised Feature Learning` Adam
% Coates, Honglak Lee, Andrew Y. Ng

if opts.contrastNormalization
  z = reshape(data,[],60000) ;
  z = bsxfun(@minus, z, mean(z,1)) ;
  n = std(z,0,1) ;
  z = bsxfun(@times, z, mean(n) ./ max(n, 40)) ;
  data = reshape(z, 32, 32, 3, []) ;
end

if opts.whitenData
  z = reshape(data,[],60000) ;
  W = z(:,set == 1)*z(:,set == 1)'/60000 ;
  [V,D] = eig(W) ;
  % the scale is selected to approximately preserve the norm of W
  d2 = diag(D) ;
  en = sqrt(mean(d2)) ;
  z = V*diag(en./max(sqrt(d2), 10))*V'*z ;
  data = reshape(z, 32, 32, 3, []) ;
end

clNames = load(fullfile(unpackPath, 'batches.meta.mat'));

imdb.images.data = data ;
imdb.images.labels = single(cat(2, labels{:})) ;
imdb.images.set = set;
imdb.meta.sets = {'train', 'val', 'test'} ;
imdb.meta.classes = clNames.label_names;

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function NSxToHL(fname)

% NSxToHL
%
% Opens and reads an NSx file without the header information and saves
% the binary data into a .dat file with the same name. This can be used for
% specific applications that require this type of data, e.g. Klusters.
% Works with File Spec 2.1, 2.2 and 2.3.
%
% It does not support pauses at this time.
%
% Use OUTPUT = NSxToHL(fname)
%
% All input arguments are optional.
%
%   fname:        Name of the file to be opened. If the fname is omitted
%                 the user will be prompted to select a file.
%                 DEFAULT: Will open Open File UI.
%
%   OUTPUT:       Contains the binary data.
%
%   Example 1:
%   NSxToHL('c:\data\sample.ns5');
%
%   In the example above, the file c:\data\sample.ns5 will be opened and
%   the data will be read and saved without the header information in file
%   sample.ns5.dat.
%
%   Example 2:
%   NSxToHL;
%
%   In the example above, the file user will be prompted for the file. The
%   selected file (e.g. FILENAME.NSx) will be opened and the data will be
%   read and saved in the file FILENAME.NSx.dat.
%
%   Kian Torab
%   kian@blackrockmicro.com
%   Blackrock Microsystems
%   Version 1.0.0.0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Version History
%
% 1.0.0.0:
%   - Initial release.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Opening the file
% Popup the Open File UI. Also, process the file name, path, and extension
% for later use, and validate the entry.
if ~exist('fname', 'var')
    if ~ismac
        [fname, path] = getFile('*.ns*', 'Choose an NSx file...');
    else
        [fname, path] = getFile('*.*', 'Choose an NSx file...');
    end
    if fname == 0
        disp('No file was selected.');
        if nargout
            clear variables;
        end
        return;
    end
else
    if isempty(fileparts(fname))
        fname = which(fname);
    end
    [path,fname, fext] = fileparts(fname);
    fname = [fname fext];
    path  = [path '/'];
end
if fname==0
    return;
end

%% Reading the headerless file
data = openNSxHL([path fname]);

%% Writing to file
% Determining the filename for the converted file
newFilename = [path fname '.dat'];

% Opening the output file for saving
FIDw = fopen(newFilename, 'w+', 'ieee-le');

% Writing data into file
disp('Writing the converted data into the new .dat file...');
fwrite(FIDw, data, 'int16');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


%==============================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% function [Sc,dS,d2S] = curvatureST(uc,omega,m,varargin)
%
% Matrix-free spatio-temporal curvature regularization energy for vc
% where vc is cell-centered
%
% Sv) = 0.5 * \int_{\omega}\int_0^1
%               alpha(1)*v(x,t)'*A*v(x,t)+ alpha(2)*v(x,t)'*B*v(x,t) dx dt,
%
% where A is the curvature operator and B the first-order time derivative
% operator.
%
% Input:
%
%   vc          instationary velocity field (cell-centered)
%   omega       spatial domain
%   m           number of discretization points in space
%   varargin    optional parameters (see below)
%
% Optional Input:
%   tspan       time span (default: [0 1])
%   nt          number of time points for velocity
%
%
% Output:
%
%   Sc          current value  (0.5 * hd * uc'* A *uc)
%   dS          derivative     (hd * uc'*A )
%   d2S         Hessian        A
%  if ~matrixFree,  d2S is sparse matrix; else, d2S is struct; endif
%
% see also curvature.m
% ==================================================================================

function [Sc,dS,d2S] = mfCurvatureST(vc,omega,m,varargin)
if nargin == 0
    help(mfilename);
    return;
end

if strcmp(vc,'para')
    Sc = 'cell-centered';       % grid
    dS = 1;                     % matrixFree
    d2S = @spectralPrecondPCG;  % solver
    return;
end

persistent A omegaOld mOld alphaOld ntOld tspanOld

if ~exist('mOld','var'),     mOld = [];     end;
if ~exist('omegaOld','var'), omegaOld = []; end;
if ~exist('alphaOld','var'), alphaOld = []; end;

matrixFree  = 0;
alpha       = [1 1e-3];
tspan       = [0 1];
nt          = [];
for k=1:2:length(varargin) % overwrites default parameter
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

dim = numel(omega)/2;

if isempty(nt) % roughly estimate nt
    nt = round(numel(vc)/(prod(m)*dim))-1;
end

d2S.regularizer = regularizer;
d2S.alpha  = alpha;
d2S.B      = @(omega,m)     getCurvatureMatrixST( omega,tspan,m,nt,alpha);
d2S.d2S    = @(u,omega,m)   curvatureOperatorST(u,omega,tspan,m,nt,alpha);
d2S.diag   = @(omega,m)     getCurvatureDiag(omega,tspan,m,nt,alpha); %getDiag(omega,tspan,m,nt,alpha);
d2S.solver = @spectralPrecondPCG;
d2S.res    = vc;
dS         = d2S.d2S(vc,omega,m)';
Sc         = .5*dS*vc;

function C = getCurvatureMatrixST(omega,tspan,m,nt,alpha)

h   = (omega(2:2:end)-omega(1:2:end))./m;
hd  = prod(h);
% compute time-stepsize
dt = abs(tspan(2)-tspan(1))/nt;

% build gradient matrix for one transformation
C =  getCurvatureMatrix(omega,m);
a = sqrt(alpha(1).*hd.*dt*[1/2;ones(nt-1,1);1/2]);
% apply spatial regularization to all transformations and sum up in
% time
C =  kron(sdiag(a(:)),C);

% get time regularization matrix
if alpha(2)>0
    b = sqrt(alpha(2)*hd.*dt);
    Bt = b * getTimeDiffusionMatrix(nt,dt,prod(m),length(omega)/2);
    C  = [C;Bt];
end

function D = sdiag(v)
	D = diag(sparse(v(:)));

% matrix free implementation of spatio-temporal curvature operator
function Ay = curvatureOperatorST(uc,omega,tspan,m,nt,alpha)
dim = numel(omega)/2;
h   = (omega(2:2:end)-omega(1:2:end))./m;
hd  = prod(h(1:dim));
dt  = abs(tspan(1)-tspan(2))/nt;
w   = dt*[1/2;ones(nt-1,1);1/2];
n   = prod(m);

uc   = reshape(uc,dim*n,[]);
ntp1 = size(uc,2);
Ay = 0*uc;

switch dim
    case 2
        D2x =  D2(1,omega,m);
        D2y =  D2(2,omega,m);

        % space: curvature
        for k=1:ntp1
            uct = reshape(uc(:,k),[m dim]);
            t1  = D2x*uct(:,:,1) + uct(:,:,1)*D2y;
            t1  = D2x*t1 + t1*D2y;
            t2  = D2x*uct(:,:,2) + uct(:,:,2)*D2y;
            t2  = D2x*t2 + t2*D2y;

            Ay(:,k) = (w(k)*hd*alpha(1))*[t1(:);t2(:)];
        end
        Ay = Ay(:);
    case 3
        % space: curvature
        d2  = @(i) D2(i,omega,m); % this is a shortcut to the discrete second derivative

        % the following line is a shortcut for
        %  - permuting the 3d-array using the permutation J
        %  - reshape it to a 2D-array of size q-by-prod(m)/q, where q=m(J(1))
        %  - multiply by A (which is q-by-q)
        %  - undo the reshape, i.e. make the result to m(J(1))-by-m(J(2))-by-m(J(3))
        %  - undo the permutation
        operate = @(A,z,J) ipermute(reshape(A*reshape(permute(z,J),m(J(1)),[]),m(J)),J);
        for k=1:ntp1 % loop over all time steps
            uct = reshape(uc(:,k),[m,dim]);
            for ell=1:3 % run over all components y^ell of Y=(Y^1,Y^2,Y^3)
                % compute
                % (I_3\otimes I_2\otimes d2(1) + I_3\otimes d2(2)\otimes I_1 ...
                % + d2(3)\otimes I_2\otimes I_1) y^ell
                %                 utt = 0;
                for rep=1:2 % go twice to have the bi-Laplacian
                    for j=1:3
                        z = operate(d2(j),uct(:,:,:,ell),[j,setdiff(1:3,j)]);
                        Ay((ell-1)*n+(1:n),k) = Ay((ell-1)*n+(1:n),k) + reshape(z,[],1);
                    end;
                    if rep==1 % overwrite uct and reset Ay. Compare to t1,t2 in dim=2 implementation
                        uct(:,:,:,ell) = reshape(Ay((ell-1)*n+(1:n),k),m);
                        Ay((ell-1)*n+(1:n),k) = 0;
                    end
                end;
            end
            Ay(:,k) = w(k)*hd*alpha(1)*Ay(:,k);
        end
        Ay = Ay(:);
    otherwise
        error('%s - dimension %d not supported.',mfilename,dim);
end
% time: diffusion
if alpha(2)>0
    d3  = @(Y) (Y(:,2:end)-Y(:,1:end-1))/dt;
    d3T = @(Y) ([-Y(:,1),Y(:,1:end-1)-Y(:,2:end),Y(:,end)])/dt;
    At = d3T(d3(uc));
    Ay = Ay(:) + (alpha(2)* hd*dt)*At(:);
end

function D = getCurvatureDiag(omega,tspan,m,nt,alpha)
% just for testing - not very efficient
dim = numel(omega)/2;
h   = (omega(2:2:end)-omega(1:2:end))./m;
hd  = prod(h(1:dim));
dt  = abs(tspan(1)-tspan(2))/nt;

C = getCurvatureMatrix(omega,m);
D = full(diag(C'*C'));

D = hd*dt*alpha(1)*repmat(D(:),1,1);
if alpha(2)>0
    errory('nyi')
end

function D = D2(i,omega,m)
h = (omega(2:2:end)-omega(1:2:end))./m;
D = spdiags(ones(m(i),1)*[1,-2,1],-1:1,m(i),m(i))/h(i)^2;
D([1,end]) = -D([2,end-1]);

function A = getTimeDiffusionMatrix(nt,dt,n,dim)
Dt=spdiags(ones(nt+1,1)*[-1/dt 1/dt],0:1,nt,nt+1);
A = kron(Dt, speye(n*dim)); % 2nd deriveative over time for two components


text:
function [zs] = ssampler(z,op1,sc)
% PURPOSE: Systematic sampling of a high-frequency time series
% ------------------------------------------------------------
% SYNTAX: zs = ssampler(z,op1,sc)
% ------------------------------------------------------------
% OUTPUT: zs: nx1 sampled time series
% ------------------------------------------------------------
% INPUT:  z: nx1 ---> vector of high frequency data
%         op1: type of temporal aggregation
%         op1=1 ---> sum (flow)
%         op1=2 ---> average (index)
%         op1=3 ---> last element (stock) ---> interpolation
%         op1=4 ---> first element (stock) ---> interpolation
%         sc: number of high frequency data points
%            for each low frequency data points
% ------------------------------------------------------------
% LIBRARY: copylow, temporal_agg
% ------------------------------------------------------------

% written by:
%  Enrique M. Quilis
%  Macroeconomic Research Department
%  Ministry of Economy and Competitiveness
%  <enrique.quilis@mineco.es>

% Version 1.0 [May 2009]

aux = temporal_agg(z,op1,sc);
zs = copylow(aux,3,sc);

author: Sable

repo: mcbench-benchmarks


text:
function gui_cg()
%
% This function is needed by FASTICAG

% This file just removes the global variables
% that are used in FASTICAG from the memory

% @(#)$Id: gui_cg.m,v 1.2 2003/04/05 14:23:57 jarmo Exp $

clear global c_FastICA_appr_strD;
clear global c_FastICA_appr_strV;
clear global c_FastICA_dMod_strD;
clear global c_FastICA_dMod_strV;
clear global c_FastICA_finetune_strD;
clear global c_FastICA_finetune_strV;
clear global c_FastICA_g1_strD;
clear global c_FastICA_g1_strV;
clear global c_FastICA_g2_strD;
clear global c_FastICA_g2_strV;
clear global c_FastICA_iSta_strD;
clear global c_FastICA_iSta_strV;
clear global c_FastICA_stabili_strD;
clear global c_FastICA_stabili_strV;
clear global c_FastICA_verb_strD;
clear global c_FastICA_verb_strV;
clear global g_FastICA_a1;
clear global g_FastICA_a2;
clear global g_FastICA_approach;
clear global g_FastICA_displayIn;
clear global g_FastICA_displayMo;
clear global g_FastICA_epsilon;
clear global g_FastICA_finetune;
clear global g_FastICA_g;
clear global g_FastICA_ica_A;
clear global g_FastICA_ica_W;
clear global g_FastICA_ica_sig;
clear global g_FastICA_initGuess;
clear global g_FastICA_initState;
clear global g_FastICA_interrupt;
clear global g_FastICA_loadType;
clear global g_FastICA_maxFinetune;
clear global g_FastICA_maxNumIte;
clear global g_FastICA_mixedmean;
clear global g_FastICA_mixedsig;
clear global g_FastICA_myy;
clear global g_FastICA_numOfIC;
clear global g_FastICA_pca_D;
clear global g_FastICA_pca_E;
clear global g_FastICA_sampleSize;
clear global g_FastICA_stabilization;
clear global g_FastICA_verbose;
clear global g_FastICA_white_dwm;
clear global g_FastICA_white_sig;
clear global g_FastICA_white_wm;
clear global hb_FastICA_initGuess;
clear global he_FastICA_a1;
clear global he_FastICA_a2;
clear global he_FastICA_displayInterval;
clear global he_FastICA_epsilon;
clear global he_FastICA_file;
clear global he_FastICA_suffix;
clear global he_FastICA_maxFinetune;
clear global he_FastICA_maxIterations;
clear global he_FastICA_myy;
clear global he_FastICA_numOfIC;
clear global he_FastICA_sampleSize;
clear global hf_FastICA_MAIN;
clear global hf_FastICA_AdvOpt;
clear global hf_FastICA_Load;
clear global hf_FastICA_Save;
clear global hpm_FastICA_approach;
clear global hpm_FastICA_displayMode;
clear global hpm_FastICA_finetune;
clear global hpm_FastICA_g;
clear global hpm_FastICA_initState;
clear global hpm_FastICA_stabilization;
clear global hpm_FastICA_verbose;
clear global ht_FastICA_dim;
clear global ht_FastICA_icaStatus;
clear global ht_FastICA_initGuess;
clear global ht_FastICA_mixedStatus;
clear global ht_FastICA_newDim;
clear global ht_FastICA_numOfSamp;
clear global ht_FastICA_whiteStatus;

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function [Ncoef, L]=frameclength(F,Ls)
%FRAMECLENGTH  Number of coefficients from length of signal
%   Usage: Ncoef=frameclength(F,Ls);
%          [Ncoef,L]=frameclength(...);
%
%   `Ncoef=frameclength(F,Ls)` returns the total number of coefficients
%   obtained by applying the analysis operator of frame *F* to a signal
%   of length *Ls* i.e. `size(frana(F,f),1)` for `Ls=length(f)`.
%
%   `[Ncoef,L]=frameclength(F,Ls)` additionally returns *L*, which is the
%   same as returned by |framelength|.
%
%   If the frame length *L* is longer than the signal length *Ls*, the
%   signal will be zero-padded to *L* by |frana|.
%
%   See also: frame, framelengthcoef

callfun = upper(mfilename);
complainif_notposint(Ls,'Ls',callfun);
complainif_notvalidframeobj(F,callfun);

L = F.length(Ls);

% Some frames need special function
if isfield(F,'clength')
    Ncoef = F.clength(L);
else
    % Generic, works for any non-realonly frame and for
    % all representaions not having any extra coefficients

    Ncoef = L*F.red;

    if F.realinput
        Ncoef=Ncoef/2;
    end

    assert(abs(Ncoef-round(Ncoef))<1e-3,...
           sprintf('%s: There is a bug. L=%d should be an integer.',...
           upper(mfilename),Ncoef));

    Ncoef=round(Ncoef);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
classdef ChiD
%%CHID Functions to handle the chi distribution.
%These distributions arise when taking the square root of the sum of the
%squares of independent normal random variables.
%Implemented methods are: mean, var, moments, PDF, CDF, (supports
%                         noncentral distributions) rand, entropy
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)

function val=mean(nu)
%%MEAN Obtain the mean of central chi probability distribution.
%
%INPUTS: nu The number of degrees of freedom of the central chi
%           distribution. Note that nu>0.
%
%OUTPUTS: val The mean of the central chi distribution.
%
%The mean of the central chi distribution is given in [1].
%
%REFERENCES:
%[1] Weisstein, Eric W. "Chi Distribution." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/ChiDistribution.html
%
%May 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=sqrt(2)*exp(gammaln((nu+1)/2)-gammaln(nu/2));
end

function val=var(nu)
%%VAR Obtain the variance of the central chi probability distribution.
%
%INPUTS: nu The number of degrees of freedom of the chi distribution. Note
%           that nu>0.
%
%OUTPUTS: val  The variance of the central chi distribution.
%
%The variance of the central chi distribution is given in [1].
%
%REFERENCES:
%[1] Weisstein, Eric W. "Chi Distribution." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/ChiDistribution.html
%
%May 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=2*exp(gammaln(1+nu/2)-gammaln(nu/2))-ChiD.mean(nu)^2;
end

function val=moments(m,nu)
%%MOMENTS Obtain noncentral moments of the chi distribution.
%
%INPUTS: m The order of the moment m>=0. The zeroth order moment is just 1.
%       nu The number of degrees of freedom of the chi distribution. Note
%          that nu>0.
%
%OUTPUTS: val  The mth-degree noncentral moments of the chi distribution.
%
%Moments of the chi distribution are given in [1].
%
%REFERENCES:
%[1] Weisstein, Eric W. "Chi Distribution." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/ChiDistribution.html
%
%May 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=2^(m/2)*exp(gammaln((1/2)*(nu+m))-gammaln(nu/2));
end

function val=PDF(x,nu)
%%PDF Evaluate the central chi probability distribution function (PDF).
%
%INPUTS: x The point(s) at which the central chi PDF is to be evaluated.
%          Note that x>=0 for a nonzero probability.
%       nu The number of degrees of freedom of the chi distribution. Note
%          that nu>0. Numerical precision problems can arise for very large
%          values of nu.
%
%OUTPUTS: val The value(s) of the central chi PDF evaluated at x.
%
%The PDF is given in [1].
%
%REFERENCES:
%[1] Weisstein, Eric W. "Chi Distribution." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/ChiDistribution.html
%
%May 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=zeros(size(x));
    sel=(x>=0);

    y=x(sel);
    val(sel)=2^(1-nu/2)*y.^(nu-1).*exp(-y.^2/2)/gamma(nu/2);
end

function prob=CDF(x,nu)
%%CDF Evaluate the cumulative distribution function (CDF) of the central
%     chi distribution at desired points.
%
%INPUTS: x The point(s) at which the chi CDF is to be evaluated.
%       nu The number of degrees of freedom of the chi-square
%          distribution. Note that nu>0.
%
%OUTPUTS: prob The value(s) of the CDF of the chi distribution with nu
%              degrees of freedom evaluated at x.
%
%The CDF of the central chi distribution is given in [1].
%
%REFERENCES:
%[1] Weisstein, Eric W. "Chi Distribution." From MathWorld--A Wolfram Web
%    Resource. http://mathworld.wolfram.com/ChiDistribution.html
%
%May 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.

    prob=gammainc(x.^2/2,nu/2);
end

function vals=rand(N,nu,lambda)
%%RAND Generate chi distributed random variables with the given parameters.
%
%INPUTS: N If N is a scalar, then rand returns an NXN matrix of random
%          variables. If N=[M,N1] is a two-element row vector, then rand
%          returns an MXN1 matrix of random variables.
%       nu The number of degrees of freedom of the chi distribution. Note
%          that nu>0. If lambda is not zero, nu >=1.
%   lambda The non-centrality parameter of the distribution. In the central
%          chi distribution, this is zero. If this parameter is omitted or
%          an empty matrix is passed, then 0 is used.
%
%OUTPUTS: vals A matrix whose dimensions are determined by N of the
%              generated chi-squared random variables.
%
%A chi distributed random variable is just the square root of a chi-squared
%random variable. Thus, this function takes the square root of
%ChiSquareD.rand(N,nu,lambda).
%
%EXAMPLE:
%Here, we demonstrate that a histogram of many random samples will tend to
%agree with the PDF of the central distribution.
% nu=4.4;
% lambda=0;
% numSamp=[10000,1];
% samp=ChiD.rand(numSamp,nu);
% figure(1)
% clf
% hold on
% h=histogram(samp,'Normalization','pdf');
% %We will plot the PDF.
% numPoints=100;
% x=linspace(0,20,numPoints);
% PDFVals=ChiD.PDF(x,nu);
% plot(x,PDFVals,'-r','linewidth',2)
%
%May 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<3||isempty(lambda))
        lambda=0;
    end

    if(isscalar(N))
        dims=[N N];
    else
        dims=N;
    end

    vals=sqrt(ChiSquareD.rand(dims,nu,lambda));
end

function entropyVal=entropy(nu)
%%ENTROPY Obtain the differential entropy of the central chi 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: nu The number of degrees of freedom of the chi distribution. Note
%           that nu>0.
%
%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=gammaln(nu/2)+(1/2)*(nu-log(2)-(nu-1)*psi(nu/2));
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 [imgRgb, overlayImg, overlayMaskImg, anatImg] = mrAnatOverlayMontage(overlayImg, xform, anatImg, anatXform, cmap, overlayClipRng, acpcSlices, fname, plane, alpha, labelFlag, upsamp, numAcross, clusterThresh, autoCropBorder)
% imgRgb = mrAnatOverlayMontage(overlayImg, xform, anatImg, anatXform, cmap, overlayClipRng, acpcSlices, [fname], [plane=3], [alpha=1], [labelFlag=true], [upsamp=0], [numAcross=[]], [clusterThresh=0], [autoCropBorder=0])
%
% Quick hack for displaying nice overlay montages.
% plane: presents sagittal (1) coronal (2) or axial (3) slices. default is
% axials.   STILL A PROBLEM WITH THE SLICE SELECTION!!!
%
% HISTORY:
% 2006.08.07 RFD: wrote it.
% 2006.08.10 MBS: added plane flag

if(~exist('cmap','var')) cmap = []; end
if(~exist('fname','var')) fname = ''; end
if(~exist('plane','var')||isempty(plane)) plane = 3; end
if(~exist('alpha','var')||isempty(alpha)) alpha = 1; end
if(~exist('labelFlag','var')||isempty(labelFlag)) labelFlag = true; end
if(size(overlayImg,4)~=3&&size(overlayImg,4)~=1)
    error('overlay must be XxYxZ or XxYxZx3!');
end
if(~exist('upsamp','var')||isempty(upsamp)) upsamp = 0; end
if(~exist('numAcross','var')) numAcross = []; end
overlayInterp = [0 0 0 0 0 0];
if(~exist('clusterThresh','var')||isempty(clusterThresh)) clusterThresh = 0; end
clusterConnectivity = 18;

% Set to Inf for no autocropping
if ~exist('autoCropBorder', 'var') || isempty(autoCropBorder)
autoCropBorder = 0;
end

[t,r,s,k] = affineDecompose(xform);
if(size(anatXform,1)==1 || size(anatXform,2)==1)
    % assume it's not a full xform but just mmPerVox. To convert it to a
    % full xform, we have to assume that it is the same as the overlay
    % xform except for maybe a scale difference.
    anatMmPerVox = anatXform(:)';
    anatXform = xform*diag([anatXform./s 1]);
else
    [at,ar,anatMmPerVox,ak] = affineDecompose(anatXform);
end

if(isfinite(autoCropBorder))
    tmp = sum(anatImg,3);
    x = find(sum(tmp,1)); x = [x(1) x(end)];
    y = find(sum(tmp,2)); y = [y(1) y(end)];
    tmp = squeeze(sum(anatImg,1));
    z = find(sum(tmp,1)); z = [z(1) z(end)];
    clear tmp;
    pad = [-autoCropBorder autoCropBorder];
    x = x+pad; y = y+pad; z = z+pad;
    anatImg = anatImg(y(1):y(2),x(1):x(2),z(1):z(2));
    anatXform = inv(inv(anatXform)*[1 0 0 -y(1); 0 1 0 -x(1); 0 0 1 -z(1); 0 0 0 1]);
end

anatSz = size(anatImg);
if(alpha>0)
    overlaySz = size(overlayImg);
    if(size(overlayImg,4)>1)
        overlayMaskImg = repmat(double(mean(overlayImg,4)>overlayClipRng(1)),[1 1 1 size(overlayImg,4)]);
    else
        overlayMaskImg = double(overlayImg>overlayClipRng(1));
    end
    if(clusterThresh>0)
        overlayMaskImg = double(bwareaopen(overlayMaskImg,clusterThresh,clusterConnectivity));
    end
    if(~all(xform(:)==anatXform(:)) || ~all(anatSz(1:3)==overlaySz(1:3)))
        %bb = dtiGet(0,'defaultBoundingBox');
        bb = anatXform*[1,1,1,1;[anatSz,1]]';
        bb = bb(1:3,:)';
        disp('Resampling overlay image to match background image resolution...');
        overlayImg = mrAnatResliceSpm(overlayImg, inv(xform), bb, anatMmPerVox, overlayInterp, false);
        overlayMaskImg = mrAnatResliceSpm(overlayMaskImg, inv(xform), bb, anatMmPerVox, overlayInterp, false);
    end
    overlayImg(overlayImg>overlayClipRng(2)) = overlayClipRng(2);
    overlayImg = (overlayImg-overlayClipRng(1))./(overlayClipRng(2)-overlayClipRng(1));
else
    overlayImg = [];
end
slXform = inv(anatXform);
if plane==3
    % axials, reorient so that the eyes point up
    perm = [2 1 3 4];
    flip = 1;
    for(ii=1:length(acpcSlices)) slLabel{ii} = sprintf('Z = %d',acpcSlices(ii)); end
    slImg = slXform*[zeros(length(acpcSlices),2) acpcSlices' ones(length(acpcSlices),1)]';
    slImg = round(slImg(3,:));
elseif plane==1
    % sagittals, looking left
    perm = [3 2 1 4];
    flip = [2,1];
    for(ii=1:length(acpcSlices)) slLabel{ii} = sprintf('X = %d',acpcSlices(ii)); end
    slImg = slXform*[acpcSlices' zeros(length(acpcSlices),2) ones(length(acpcSlices),1)]';
    slImg = round(slImg(1,:));
elseif plane==2
    % coronals
    perm = [3 1 2 4];
    flip = 1;
    for(ii=1:length(acpcSlices)) slLabel{ii} = sprintf('Y = %d',acpcSlices(ii)); end
    slImg = slXform*[zeros(length(acpcSlices),1) acpcSlices' zeros(length(acpcSlices),1) ones(length(acpcSlices),1)]';
    slImg = round(slImg(2,:));
end
anatImg = flipdim(permute(anatImg,perm),flip(1));
if(numel(flip)>1), anatImg = flipdim(anatImg,flip(2)); end
if(~isempty(overlayImg))
    overlayImg = flipdim(permute(overlayImg,perm),flip(1));
    overlayMaskImg = flipdim(permute(overlayMaskImg,perm),flip(1));
    if(numel(flip)>1)
        overlayImg = flipdim(overlayImg,flip(2));
        overlayMaskImg = flipdim(overlayMaskImg,flip(2));
    end
end

if(~labelFlag) slLabel = ''; end

[imgRgb,junk,numColsRows] = makeMontage3(anatImg, slImg, anatMmPerVox(1), upsamp, [], numAcross, 0);
if(~isempty(overlayImg))
    overlayMask = makeMontage3(overlayMaskImg, slImg, anatMmPerVox(1), upsamp, [], numAcross, 0);
    overlayMask = overlayMask>=0.5;
    overlay = makeMontage3(overlayImg, slImg, anatMmPerVox(1), upsamp, [], numAcross, 0);
    if(size(overlayImg,4)==3)
        if(alpha>1)
            % alpha>1 applies a gamma to the transparency map
            g = 1/alpha;
            alpha = 1;
            lum = repmat(mean(overlay,3).^g,[1,1,3]);
            imgRgb(overlayMask) = (1-lum(overlayMask).*alpha).*imgRgb(overlayMask) + lum(overlayMask).*alpha.*overlay(overlayMask);
        else
            imgRgb(overlayMask) = (1-alpha).*imgRgb(overlayMask) + alpha.*overlay(overlayMask);
        end
    else
        if(isempty(cmap)), cmap = autumn(256); end
        overlay = mean(overlay,3);
        %n = isnan(overlay);
        %overlay(n) = 0;
        %if(size(overlay,3)>1)
        %    n = sum(n,3); n(n==0) = NaN;
        %    overlay = sum(overlay,3)./sum(n,3);
        %    overlay(isnan(n)) = 0;
        %end
        inds = round(overlay*255+1);
        inds(isnan(inds)) = 1;
        overlay = reshape(cmap(inds,:),[size(overlay) 3]);
        overlay(~overlayMask|isnan(overlay)) = 0;
        imgRgb(overlayMask) = (1-alpha).*imgRgb(overlayMask) + alpha.*overlay(overlayMask);
    end
end
figure;
if(exist('imshow.m','file'))
  iptsetpref('ImshowBorder','tight')
  imshow(imgRgb);
else
  image(imgRgb); axis equal; axis off;
  set(gca,'Position',[0,0,1,1]);
end
set(gcf,'userData',imgRgb);
disp('im=get(gcf,''userData''); % get the raw image data');
mrUtilLabelMontage(slLabel, numColsRows, gcf, gca);

if(~isempty(fname)) mrUtilPrintFigure(fname); end
if(~isempty(overlayImg))
    legendLabels = explode(',',sprintf('%0.0f,',linspace(overlayClipRng(1),overlayClipRng(2),5)));
    legendLabels = legendLabels(1:end-1);
    legendLabels{end} = ['>=' num2str(overlayClipRng(2),'%0.3g')];
    if(~isempty(cmap))
        %if(~isempty(fname)) mrUtilMakeColorbar(cmap, legendLabels, '', [fname '_legend']);
        %else  mrUtilMakeColorbar(cmap, legendLabels, ''); end
    end
end

% be kind to those who forget a semicolon:
if(nargout<1) clear imgRgb; end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% load 'RSS_data.mat';%original rss signals
% load 'B.mat';% smoothing rss signals
[action_size,sample_size]=size(RSS_data);
for i=1:action_size
    for j=1:1
               rss_sample=RSS_data{i,j};
%                rss_smooth_sample_medfilt=medfilt1(RSS_data{i,j},5);
%                rss_smooth_sample =B{i,j};
    end
    subplot(4,4,i);
%     plot(rss_sample,'r');
%     hold on;
    plot(rss_sample,'k','LineWidth',2,'MarkerSize',12);
      xlabel('Packets');
     ylabel('RSS','fontsize',8);
%      legend('Original Signals','Smooth Signals','location','best');

end
figure;

for i=1:16
    for j=1:1
               csi_sample=CSI_data{i,j};
%                rss_smooth_sample_medfilt=medfilt1(RSS_data{i,j},5);
%                rss_smooth_sample =B{i,j};
              csi_sample_re=permute(csi_sample,[2,1,3]);
    end
    subplot(4,4,i);
%     plot(rss_sample,'r');
%     hold on;
    plot(csi_sample_re(:,10,1),'k','LineWidth',2,'MarkerSize',12);
      xlabel('Packets');
     ylabel('CSI','fontsize',8);
%      legend('Original Signals','Smooth Signals','location','best');

end

author: linteresa

repo: WiAR

save_path: github-repos/MATLAB/linteresa-WiAR


text:
% RBio: MATLAB toolbox for reading/writing sparse matrices in the Rutherford/
%   Boeing format, and for reading/writing problems in the UF Sparse Matrix
%   Collection from/to a set of files in a directory.
%
%   RBread    - read a sparse matrix from a Rutherford/Boeing file
%   RBreade   - read a symmetric finite-element matrix from a R/B file
%   RBtype    - determine the Rutherford/Boeing type of a sparse matrix
%   RBwrite   - write a sparse matrix to a Rutherford/Boeing file
%   RBraw     - read the raw contents of a Rutherford/Boeing file
%   RBfix     - read a possibly corrupted matrix from a R/B file
%   RBinstall - install the RBio toolbox for use in MATLAB
%   RBmake    - compile the RBio toolbox for use in MATLAB
%
% Example:
%
%   load west0479
%   C = west0479 ;
%   RBwrite ('mywest', C, 'WEST0479 chemical eng. problem', 'west0479')
%   A = RBread ('mywest') ;
%   norm (A-C,1)
%
% See also UFget, mread, mwrite.

% Copyright 2007, Timothy A. Davis


text:
%  Figure 10.11      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
%  fig10_11.m is a script to generate Fig. 10.11
%  the frequency response of the notch network
clf;
nnotch=[1/.81 0  1] ;
dnotch=[1/625 2/25  1];
% define frequency range
w=logspace(-1,1);
w(36)=1;
subplot(211)
w=logspace(-1,1);
w(24)=0.89;
% compute Bode
[magn phn]=bode(nnotch,dnotch,w);
% phn=php-360*ones(phn);
magn1=[magn, ones(size(magn))];
subplot(211); loglog(w,magn1); grid;
xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 10.11 Bode plot of a notch filter')
subplot(212); semilogx(w,phn); grid;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');


text:
classdef String
  % String utils

  % comma-separated list of items
  methods (Static)

    function c = commaSepList2CellStr(s)
      c = regexp(s,',','split');
      c = strtrim(c);
    end

    function s = cellstr2CommaSepList(c)
      if isempty(c)
        s = '';
      else
        s = [sprintf('%s,',c{1:end-1}) c{end}];
      end
    end

    function s = cellstr2DelimList(c,d)
      if isempty(c)
        s = '';
      else
        pat = sprintf('%%s%s',d);
        s = [sprintf(pat,c{1:end-1}) c{end}];
      end
    end

    % see civilizedStringFromCellArrayOfStrings

    function s = niceUpperCase(s)
      s = lower(s);
      if ~isempty(s)
        s(1) = upper(s(1));
      end
    end

  end

end

author: kristinbranson

repo: JAABA


text:
% Jiao Xianjun (putaoshu@msn.com; putaoshu@gmail.com)
% Generate pss time domain sequences under different frequency offsets.
% A script of project: https://github.com/JiaoXianjun/rtl-sdr-LTE

function pss_fo_set = pss_fo_set_gen(pss, fo_search_set)
num_pss = size(pss, 2);
len_pss = size(pss, 1);

sampling_rate = 1.92e6; % LTE spec
num_fo = length(fo_search_set);

pss_fo_set = zeros(len_pss, num_fo*num_pss);
for i=1:num_pss
    sp = (i-1)*num_fo + 1;
    ep = sp + num_fo - 1;
    pss_fo_set(:,sp:ep) = kron(ones(1, num_fo), pss(:,i)).*exp(1i.*2.*pi.*(1./sampling_rate).*(0:(len_pss-1)).'*fo_search_set);
    pss_fo_set(:,sp:ep) = conj(pss_fo_set(:,sp:ep))./len_pss;
end

% % normalize along with column
% pss_fo_set = sqrt(len_pss).*pss_fo_set./kron( ones(len_pss,1), sqrt( sum(abs(pss_fo_set).^2,1) ) );

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function obj = updateSurfl(obj, surfaceIndex)

%-AXIS INDEX-%
axIndex = obj.getAxisIndex(obj.State.Plot(surfaceIndex).AssociatedAxis);

%-CHECK FOR MULTIPLE AXES-%
[xsource, ysource] = findSourceAxis(obj,axIndex);

%-SURFACE DATA STRUCTURE- %
meshData = get(obj.State.Plot(surfaceIndex).Handle);
figureData = get(obj.State.Figure.Handle);

%-AXIS STRUCTURE-%
axisData = get(ancestor(meshData.Parent,'axes'));

%-SCENE DATA-%
eval( sprintf('scene = obj.layout.scene%d;', xsource) );

%-GET CONTOUR INDEX-%
obj.PlotOptions.nPlots = obj.PlotOptions.nPlots + 1;
contourIndex = obj.PlotOptions.nPlots;
obj.PlotOptions.contourIndex(surfaceIndex) = contourIndex;

%-------------------------------------------------------------------------%

%-associate scene-%
obj.data{surfaceIndex}.scene = sprintf('scene%d', xsource);
obj.data{contourIndex}.scene = sprintf('scene%d', xsource);

%-------------------------------------------------------------------------%

%-surface type for face color-%
obj.data{surfaceIndex}.type = 'surface';

%-scatter3d type for contour mesh lines-%
obj.data{contourIndex}.type = 'scatter3d';
obj.data{contourIndex}.mode = 'lines';

%-------------------------------------------------------------------------%

%-get plot data-%
xData = meshData.XData;
yData = meshData.YData;
zData = meshData.ZData;

%-reformat data to mesh-%
xDataSurface = xData;
yDataSurface = yData;
zDataSurface = zData;

xDataContourDir1 = [xDataSurface; NaN(1, size(xDataSurface, 2))];
yDataContourDir1 = [yDataSurface; NaN(1, size(yDataSurface, 2))];
zDataContourDir1 = [zDataSurface; NaN(1, size(zDataSurface, 2))];

xDataContourDir2 = xDataContourDir1(1:end-1,:)';
yDataContourDir2 = yDataContourDir1(1:end-1,:)';
zDataContourDir2 = zDataContourDir1(1:end-1,:)';

xDataContourDir2 = [xDataContourDir2; NaN(1, size(xDataContourDir2, 2))];
yDataContourDir2 = [yDataContourDir2; NaN(1, size(yDataContourDir2, 2))];
zDataContourDir2 = [zDataContourDir2; NaN(1, size(zDataContourDir2, 2))];

xDataContour = [xDataContourDir1(:); xDataContourDir2(:)];
yDataContour = [yDataContourDir1(:); yDataContourDir2(:)];
zDataContour = [zDataContourDir1(:); zDataContourDir2(:)];

%-------------------------------------------------------------------------%

%-set data on surface-%
obj.data{surfaceIndex}.x = xDataSurface;
obj.data{surfaceIndex}.y = yDataSurface;
obj.data{surfaceIndex}.z = zDataSurface;

%- setting grid mesh by default -%
% x-direction
xData = xData(1, :);
obj.data{surfaceIndex}.contours.x.start = xData(1);
obj.data{surfaceIndex}.contours.x.end = xData(end);
obj.data{surfaceIndex}.contours.x.size = mean(diff(xData));
obj.data{surfaceIndex}.contours.x.show = true;

% y-direction
yData = yData(:, 1);
obj.data{surfaceIndex}.contours.y.start = yData(1);
obj.data{surfaceIndex}.contours.y.end = yData(end);
obj.data{surfaceIndex}.contours.y.size = mean(diff(yData));;
obj.data{surfaceIndex}.contours.y.show = true;

%-------------------------------------------------------------------------%

%-set data on scatter3d-%
obj.data{contourIndex}.x = xDataContour(:);
obj.data{contourIndex}.y = yDataContour(:);
obj.data{contourIndex}.z = zDataContour(:);

%-------------------------------------------------------------------------%

%-COLORING-%

%-------------------------------------------------------------------------%

%-get colormap-%
cMap = figureData.Colormap;
fac = 1/(length(cMap)-1);
colorScale = {};

for c = 1: length(cMap)
    colorScale{c} = { (c-1)*fac , sprintf('rgb(%f,%f,%f)', 255*cMap(c, :))};
end

%-------------------------------------------------------------------------%

%-get edge color-%
if isnumeric(meshData.EdgeColor)
    cDataContour = sprintf('rgb(%f,%f,%f)', 255*meshData.EdgeColor);

elseif strcmpi(meshData.EdgeColor, 'interp')
    cDataContour = zDataContour(:);
    obj.data{contourIndex}.line.colorscale = colorScale;

    obj.data{surfaceIndex}.contours.x.show = false;
    obj.data{surfaceIndex}.contours.y.show = false;

elseif strcmpi(meshData.EdgeColor, 'flat')
    cData = meshData.CData;

    if size(cData, 3) ~= 1
        cMap = unique( reshape(cData, ...
            [size(cData,1)*size(cData,2), size(cData,3)]), 'rows' );
        cData = rgb2ind(cData, cMap);

        edgeColorScale = {};
        fac = 1/(length(cMap)-1);

        for c = 1: length(cMap)
            edgeColorScale{c} = { (c-1)*fac , sprintf('rgb(%f,%f,%f)', 255*cMap(c, :))};
        end

        obj.data{surfaceIndex}.line.cmin = 0;
        obj.data{surfaceIndex}.line.cmax = 255;
        obj.data{contourIndex}.line.colorscale = edgeColorScale;
    else
        obj.data{contourIndex}.line.cmin = axisData.CLim(1);
        obj.data{contourIndex}.line.cmax = axisData.CLim(2);
        obj.data{contourIndex}.line.colorscale = colorScale;
    end

    cDataContourDir1 = [cData; NaN(1, size(cData, 2))];
    cDataContourDir2 = cDataContourDir1(1:end-1,:)';
    cDataContourDir2 = [cDataContourDir2; NaN(1, size(cDataContourDir2, 2))];
    cDataContour = [cDataContourDir1(:); cDataContourDir2(:)];

    obj.data{surfaceIndex}.contours.x.show = false;
    obj.data{surfaceIndex}.contours.y.show = false;

elseif strcmpi(meshData.EdgeColor, 'none')
    cDataContour = 'rgba(0,0,0,0)';
    obj.data{surfaceIndex}.contours.x.show = false;
    obj.data{surfaceIndex}.contours.y.show = false;

end

%-set edge color-%
obj.data{contourIndex}.line.color = cDataContour;
obj.data{surfaceIndex}.contours.x.color = cDataContour;
obj.data{surfaceIndex}.contours.y.color = cDataContour;

%-------------------------------------------------------------------------%

%-get face color-%
faceColor = meshData.FaceColor;

if isnumeric(faceColor)

    if all(faceColor == [1, 1, 1])
        faceColor = [0.96, 0.96, 0.96];
    end

    for n = 1:size(zDataSurface, 2)
        for m = 1:size(zDataSurface, 1)
            cDataSurface(m, n, :) = faceColor;
        end
    end

    [cDataSurface, cMapSurface] = rgb2ind(cDataSurface, 256);
    cDataSurface = double(cDataSurface) + axisData.CLim(1);

    for c = 1: size(cMapSurface, 1)
        colorScale{c} = { (c-1)*fac , sprintf('rgba(%f,%f,%f, 1)', cMapSurface(c, :))};
    end

    obj.data{surfaceIndex}.cmin = axisData.CLim(1);
    obj.data{surfaceIndex}.cmax = axisData.CLim(2);

elseif strcmpi(faceColor, 'interp')
    cDataSurface = zDataSurface;

    if surfaceIndex > xsource
        cData = [];

        for idx = xsource:surfaceIndex
            cData = [cData; obj.data{idx}.z];
        end

        cMin = min(cData(:));
        cMax = max(cData(:));

        for idx = xsource:surfaceIndex
            obj.data{idx}.cmin = cMin;
            obj.data{idx}.cmax = cMax;
        end
    end

elseif strcmpi(faceColor, 'flat')
    cData = meshData.CData;

    if size(cData, 3) ~= 1
        cMap = unique( reshape(cData, ...
            [size(cData,1)*size(cData,2), size(cData,3)]), 'rows' );
        cDataSurface = rgb2ind(cData, cMap);

        colorScale = {};
        fac = 1/(length(cMap)-1);

        for c = 1: length(cMap)
            colorScale{c} = { (c-1)*fac , sprintf('rgb(%f,%f,%f)', 255*cMap(c, :))};
        end
    else
        cDataSurface = cData;
        obj.data{surfaceIndex}.cmin = axisData.CLim(1);
        obj.data{surfaceIndex}.cmax = axisData.CLim(2);
    end

end

%-set face color-%
obj.data{surfaceIndex}.colorscale = colorScale;
obj.data{surfaceIndex}.surfacecolor = cDataSurface;

%-lighting settings-%

if isnumeric(meshData.FaceColor) && all(meshData.FaceColor == [1, 1, 1])
    obj.data{surfaceIndex}.lighting.diffuse = 0.5;
    obj.data{surfaceIndex}.lighting.ambient = 0.725;
else
    % obj.data{surfaceIndex}.lighting.diffuse = 1.0;
    % obj.data{surfaceIndex}.lighting.ambient = 0.9;
end

if meshData.FaceAlpha ~= 1
    obj.data{surfaceIndex}.lighting.diffuse = 0.5;
    obj.data{surfaceIndex}.lighting.ambient = 0.725 + (1-meshData.FaceAlpha);
end

if obj.PlotlyDefaults.IsLight
    obj.data{surfaceIndex}.lighting.diffuse = 1.0;
    obj.data{surfaceIndex}.lighting.ambient = 0.3;
end

%-opacity-%
obj.data{surfaceIndex}.opacity = meshData.FaceAlpha;

%-------------------------------------------------------------------------%

%-line style-%

obj.data{contourIndex}.line.width = 3*meshData.LineWidth;

if strcmpi(meshData.LineStyle, '-')
    obj.data{contourIndex}.line.dash = 'solid';
else
    obj.data{contourIndex}.line.dash = 'dot';
    obj.data{surfaceIndex}.contours.x.show = false;
    obj.data{surfaceIndex}.contours.y.show = false;
end

%-------------------------------------------------------------------------%

%-SCENE CONFIGUTATION-%

%-------------------------------------------------------------------------%

%-aspect ratio-%
asr = obj.PlotOptions.AspectRatio;

if ~isempty(asr)
    if ischar(asr)
        scene.aspectmode = asr;
    elseif isvector(ar) && length(asr) == 3
        xar = asr(1);
        yar = asr(2);
        zar = asr(3);
    end
else

    %-define as default-%
    xar = max(xData(:));
    yar = max(yData(:));
    xyar = max([xar, yar]);
    zar = 0.75*xyar;
end

scene.aspectratio.x = 1.1*xyar;
scene.aspectratio.y = 1.0*xyar;
scene.aspectratio.z = zar;

%---------------------------------------------------------------------%

%-camera eye-%
ey = obj.PlotOptions.CameraEye;

if ~isempty(ey)
    if isvector(ey) && length(ey) == 3
        scene.camera.eye.x = ey(1);
        scene.camera.eye.y = ey(2);
        scene.camera.eye.z = ey(3);
    end
else

    %-define as default-%
    xey = - xyar; if xey>0 xfac = -0.0; else xfac = 0.0; end
    yey = - xyar; if yey>0 yfac = -0.3; else yfac = 0.3; end
    if zar>0 zfac = -0.1; else zfac = 0.1; end

    scene.camera.eye.x = xey + xfac*xey;
    scene.camera.eye.y = yey + yfac*yey;
    scene.camera.eye.z = zar + zfac*zar;
end

%-------------------------------------------------------------------------%

%-scene axis configuration-%

scene.xaxis.range = axisData.XLim;
scene.yaxis.range = axisData.YLim;
scene.zaxis.range = axisData.ZLim;

scene.xaxis.tickvals = axisData.XTick;
scene.xaxis.ticktext = axisData.XTickLabel;

scene.yaxis.tickvals = axisData.YTick;
scene.yaxis.ticktext = axisData.YTickLabel;

scene.zaxis.tickvals = axisData.ZTick;
scene.zaxis.ticktext = axisData.ZTickLabel;

scene.xaxis.zeroline = false;
scene.yaxis.zeroline = false;
scene.zaxis.zeroline = false;

scene.xaxis.showline = true;
scene.yaxis.showline = true;
scene.zaxis.showline = true;

scene.xaxis.tickcolor = 'rgba(0,0,0,1)';
scene.yaxis.tickcolor = 'rgba(0,0,0,1)';
scene.zaxis.tickcolor = 'rgba(0,0,0,1)';

scene.xaxis.ticklabelposition = 'outside';
scene.yaxis.ticklabelposition = 'outside';
scene.zaxis.ticklabelposition = 'outside';

scene.xaxis.title = axisData.XLabel.String;
scene.yaxis.title = axisData.YLabel.String;
scene.zaxis.title = axisData.ZLabel.String;

scene.xaxis.tickfont.size = axisData.FontSize;
scene.yaxis.tickfont.size = axisData.FontSize;
scene.zaxis.tickfont.size = axisData.FontSize;

scene.xaxis.tickfont.family = matlab2plotlyfont(axisData.FontName);
scene.yaxis.tickfont.family = matlab2plotlyfont(axisData.FontName);
scene.zaxis.tickfont.family = matlab2plotlyfont(axisData.FontName);

%-------------------------------------------------------------------------%

%-SET SCENE TO LAYOUT-%
obj.layout = setfield(obj.layout, sprintf('scene%d', xsource), scene);

%-------------------------------------------------------------------------%

%-surface name-%
obj.data{surfaceIndex}.name = meshData.DisplayName;
obj.data{contourIndex}.name = meshData.DisplayName;

%-------------------------------------------------------------------------%

%-surface showscale-%
obj.data{surfaceIndex}.showscale = false;
obj.data{contourIndex}.showscale = false;

%-------------------------------------------------------------------------%

%-surface visible-%
obj.data{surfaceIndex}.visible = strcmp(meshData.Visible,'on');
obj.data{contourIndex}.visible = strcmp(meshData.Visible,'on');

%-------------------------------------------------------------------------%

leg = get(meshData.Annotation);
legInfo = get(leg.LegendInformation);

switch legInfo.IconDisplayStyle
    case 'on'
        showleg = true;
    case 'off'
        showleg = false;
end

obj.data{surfaceIndex}.showlegend = showleg;

%-------------------------------------------------------------------------%

end

author: plotly


text:
function interp1_table1_import

libname = 'interp1_table1.a';

type11r = ['double[K1] r_ck, double[K1] i_ck, int32 K1, ' ...
	'double[J1*L1+1] r_h1, int32 J1, int32 L1, ' ...
	'double[M] p_tm, int32 M, double[M] &r_fm, double[M] &i_fm'];

myimport(libname, 'interp1_table1_real_per', 'void', type11r);

function myimport(libname, name, type_return, type_call)
import(libname, name, name, type_return, type_call);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function centerfigureonscreen(hFig)

screenSize = get(groot, 'ScreenSize');
hFig.Position = [screenSize(3)/2-hFig.Position(3)/2 ...
                 screenSize(4)/2-hFig.Position(4)/2 ...
                 hFig.Position(3) ...
                 hFig.Position(4)];

end

author: pglira


text:
function [endmembers,I,Y,R_gt,A_gt,names,wl,Y_noise] = prepare_supervised_unmixing(dataset, options)
%PREPARE_SUPERVISED_UNMIXING Summary of this function goes here
%   Detailed explanation goes here
if nargin < 2
    options = [];
end

if nargin < 1
    dataset = '03';
end

force_positive_I = 0;

if isstruct(options)
    arg_set = fieldnames(options);
    for i = 1:length(arg_set)
        eval([arg_set{i},'=options.',arg_set{i},';']);
    end
end

Y_noise = 0;

switch dataset
    case '00'
        load('toy_image_end_var_3.mat');
    case '001'
        load('toy_image_end_var_3.mat');
    case '002'
        load('toy_image_end_var_3_B.mat');
    case '01'
        if 0
            options = [];
            options.chosen = [1 3];
            [I,E_gt,A_gt,names,wl,Y_noise,extra] = generate_toy_image_unified(...
                60,60,'image','dirichlet',60,options);
            pure_spectra = extra.pure_spectra;
            save('toy_image_end_var_pavia_A_spectra_snr60.mat','I','E_gt',...
                'A_gt','names','wl','pure_spectra','Y_noise');
        else
%             load('toy_image_end_var_pavia_A_spectra_13_snr20.mat');
            load('C:\Data\GMM\Y_noise_10\toy_image_end_var_pavia_A_spectra_13_no_1.mat');
        end
    case '02'
        load('../GMM/toy_image_end_var_snr60_1231.mat');
    case '03'
        if 0
            options = [];
            options.chosen = [6 2];
            [I,E_gt,A_gt,names,wl,Y_noise,extra] = generate_toy_image_unified(...
                60,60,'gmm','dirichlet',60,options);
            ws = extra.ws;
            mus = extra.mus;
            sigmas = extra.sigmas;
            save('toy_image_end_var_dirichlet_snr60.mat','I','E_gt',...
                'A_gt','names','wl','ws','mus','sigmas','Y_noise');
        else
            load('toy_image_end_var_dirichlet_snr60.mat');
        end
    case '21'
        load('../../Data/PaviaUniversity_A.mat');
        A_gt = double(A_gt);
    case '62'
        load('../../data/muufl_gulfport_B.mat');
        A_gt = double(A_gt);
    case '7'
        addpath('../../../HyperspectralImageAnalysis/Project2');
        [I,wl,pure_spectra,names] = load_santabar_spectra_lib('91');
        A_gt = zeros(size(I,1), size(I,2), length(pure_spectra));
        R_gt = zeros(length(pure_spectra), size(I,3));
    case '71'
        addpath('../../../HyperspectralImageAnalysis/Project2');
        [I,wl,pure_spectra,names] = load_santabar_spectra_lib('259');
        A_gt = zeros(size(I,1), size(I,2), length(pure_spectra));
        R_gt = zeros(length(pure_spectra), size(I,3));
    otherwise
end

if force_positive_I
    I(I<0) = 0;

    % for BCM, the first band of the gulfport dataset contains too many 0s
    % s.t. one endmember has all the pure pixels give 0 at the first band
    if strcmp(dataset,'62')
        I(:,:,1) = I(:,:,2);
        I(:,:,3) = I(:,:,4);
    end
end

[Y,A_gt,rows,cols] = reshape_hsi(I,A_gt);
[N,M] = size(A_gt);
[N,B] = size(Y);

switch dataset
    case {'0','01','7','71'} % pure_spectra is a cell array of endmember spectra
        endmembers = pure_spectra;
    case {'03'} % slice E_gt to get the endmember spectra
        endmembers = cell(1,M);
        for j = 1:M
            endmembers{j} = reshape(E_gt(j,:,:),B,N)';
        end
    otherwise % A_gt uses 1 to mark pure pixels
        endmembers = cell(1,M);
        for j = 1:M
            X = Y(A_gt(:,j)==1, :);
            endmembers{j} = X;
        end
end

R_gt = zeros(M,B);
for j = 1:M
    R_gt(j,:) = mean(endmembers{j});
end

author: zhouyuanzxcv

repo: Hyperspectral


text:
function W = tess_smooth_sources(Vertices, Faces, VertConn, FWHM, Method)
% TESS_SMOOTH_SOURCES: Gaussian smoothing matrix over a mesh.
%
% USAGE:  W = tess_smooth_sources(Vertices, Faces, VertConn=[], FWHM=0.010, Method='average')
%
% INPUT:
%    - Vertices : Vertices positions ([X(:) Y(:) Z(:)])
%    - Faces    : Triangles matrix
%    - VertConn : Vertices connectivity, logical sparse matrix [Nvert,Nvert]
%    - FWHM     : Full width at half maximum, in meters (default=0.010)
%    - Method   : {'euclidian', 'path', 'average', 'surface'}
% OUPUT:
%    - W: smoothing matrix (sparse)
%
% DESCRIPTION:
%    - The distance between two points is an average of:
%        - the direct euclidian between the two points and
%        - the number of edges between the two points * the average length of an edge
%    - Gaussian smoothing function on the euclidian distance:
%      f(r) = 1 / sqrt(2*pi*sigma^2) * exp(-(r.^2/(2*sigma^2)))
%    - Full Width at Half Maximum (FWHM) is related to sigma by:
%      FWHM = 2 * sqrt(2*log2(2)) * sigma

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2010-2013

% ===== PARSE INPUTS =====
if (nargin < 5) || isempty(Method)
    Method = 'average';
end
if (nargin < 4) || isempty(FWHM)
    FWHM = 0.010;
end
if (nargin < 3) || isempty(VertConn)
    VertConn = tess_vertconn(Vertices, Faces);
end
if ~islogical(VertConn)
    error('Invalid vertices connectivity matrix.');
end
nv = size(Vertices,1);

% ===== ANALYZE INPUT =====
% Calculate Gaussian kernel properties
Sigma = FWHM / (2 * sqrt(2*log2(2)));
% FWTM = 2 * sqrt(2*log2(10)) * Sigma;
% Get the average edge length
[vi,vj] = find(VertConn);
meanDist = mean(sqrt((Vertices(vi,1) - Vertices(vj,1)).^2 + (Vertices(vi,2) - Vertices(vj,2)).^2 + (Vertices(vi,3) - Vertices(vj,3)).^2));
% Guess the number of iterations
nIter = min(10, ceil(FWHM / meanDist));

% ===== COMPUTE DISTANCE =====
switch lower(Method)
    % === METHOD 1: USE EUCLIDIAN DISTANCE ===
    case 'euclidian'
        % Get the neighborhood around each vertex
        VertConn = mpower(VertConn, nIter);
        [vi,vj] = find(VertConn);
        % Use Euclidean distance
        d = sqrt((Vertices(vi,1) - Vertices(vj,1)).^2 + (Vertices(vi,2) - Vertices(vj,2)).^2 + (Vertices(vi,3) - Vertices(vj,3)).^2);
        Dist = sparse(vi, vj, d, nv, nv);

    % === METHOD 2: USE NUMBER OF CONNECTIONS =====
    % === METHOD 3: AVERAGE METHOD 1+2 ===
    case {'path', 'average'}
        % Initialize loop variables
        VertConnGrow = speye(nv);
        VertIter = sparse(nv,nv);
        vall = [];

        % Grow and keep track of the layers
        for iter = 1:nIter
            disp(sprintf('SMOOTH> Iteration %d/%d', iter, nIter));
            % Grow selection of vertices
            VertConnPrev = VertConnGrow;
            VertConnGrow = double(VertConnGrow * VertConn > 0);
            % Find all the new connections
            vind = find(VertConnGrow - VertConnPrev > 0);
            [vi,vj] = ind2sub([nv,nv], vind);
            VertIter = VertIter + iter * sparse(vi, vj, ones(size(vi)), nv, nv);
        end

        % Use distance in number of connections
        Dist = VertIter * meanDist;
        Dist(1:nv+1:nv*nv) = 0;

        % == AVERAGE WITH METHOD 1 ==
        if strcmpi(Method, 'average')
            % Calculate Euclidean distance
            [vi,vj] = find(VertConnGrow);
            d = sqrt((Vertices(vi,1) - Vertices(vj,1)).^2 + (Vertices(vi,2) - Vertices(vj,2)).^2 + (Vertices(vi,3) - Vertices(vj,3)).^2);
            % Average with results of method #2
            Dist = (0.5 .* Dist + 0.5 .* sparse(vi, vj, d, nv, nv));
        end

    % ===== METHOD 4: CALCULATE SURFACE DISTANCE =====
    % WARNING: NOT FINISHED!!!!
    case 'surface'
        % Initialize loop variables
        VertConnGrow = speye(nv);
        Dist = sparse([], [], [], nv, nv, 3*nnz(VertConn));
        vall = [];
        nIter = 2;
        % Grow until we reach an accepteable distance
        for iter = 1:nIter
            disp(sprintf('Iteration %d', iter));
            % Get neighbors
            VertConnGrow = VertConnGrow * VertConn;
            % Get all the existing edges in the surface
            vind = find(VertConnGrow);
            % Remove all the previously processed connections
            vind = setdiff(vind, vall);
            [vi,vj] = ind2sub([nv,nv], vind);
            % Remove diagonal
            iDel = (vi == vj);
            vi(iDel) = [];
            vj(iDel) = [];
            % Calculate the distance to the neighbor nodes
            if (iter == 1)
                % Calculate all the distances for all the pairs of edges
                d = sqrt((Vertices(vi,1) - Vertices(vj,1)).^2 + (Vertices(vi,2) - Vertices(vj,2)).^2 + (Vertices(vi,3) - Vertices(vj,3)).^2);
            else
                % Initialize d matrix
                d = zeros(size(vi));
                % Process each new connection separately
                for i = 1:length(vi)
                    % Find nodes that are connected to both nodes
                    iMid = find((Dist(vi(i),:) & VertConn(vj(i),:)) | (VertConn(vi(i),:) & Dist(vj(i),:)));
                    % Find nodes for which we know one connection at least
                    iMid0 = (Dist(vi(i),iMid) & Dist(vj(i),iMid));
                    iMid1 = (Dist(vi(i),iMid) & ~iMid0);
                    iMid2 = (Dist(vj(i),iMid) & ~iMid0);
                    % Compute distances
                    dMid = 0*iMid;
                    dMid(iMid0) = Dist(vi(i),iMid0) + Dist(vj(i),iMid0);
                    dMid(iMid1) = Dist(vi(i),iMid1) + sqrt((Vertices(vj(i),1) - Vertices(iMid1,1)).^2 + (Vertices(vj(i),2) - Vertices(iMid1,2)).^2 + (Vertices(vj(i),3) - Vertices(iMid1,3)).^2)';
                    dMid(iMid2) = Dist(vj(i),iMid2) + sqrt((Vertices(vi(i),1) - Vertices(iMid2,1)).^2 + (Vertices(vi(i),2) - Vertices(iMid2,2)).^2 + (Vertices(vi(i),3) - Vertices(iMid2,3)).^2)';
                    dMid(dMid == 0) = Inf;
                    % Find the shortest path
                    d(i) = min(dMid);
                    if isinf(d(i))
                        error('???');
                    end
                end
            end
            % Add to processed vertices
            vall = union(vind, vall);
            % Create a sparse distance matrix
            Dist = Dist + sparse(vi, vj, d, nv, nv);
        end
end

% ===== APPLY GAUSSIAN FUNCTION =====
% Gaussian function
fun = inline('1 / sqrt(2*pi*sigma2) * exp(-(x.^2/(2*sigma2)))', 'x', 'sigma2');
% Calculate interpolation as a function of distance
[vi,vj] = find(Dist>0);
vind = sub2ind([nv,nv], vi, vj);
w = fun(Dist(vind), Sigma^2);
% Build final symmetric matrix
%W = sparse([vi;vj], [vj;vi], [w;w], nv, nv);
W = sparse(vi, vj, w, nv, nv);
% Add the diagonal
W = W + fun(0,Sigma^2) * speye(nv);
% Normalize columns
W = bst_bsxfun(@rdivide, W, sum(W,1));
% Remove insignificant values
[vi,vj] = find(W>0.005);
vind = sub2ind([nv,nv], vi, vj);
W = sparse(vi, vj, W(vind), nv, nv);

% ===== FIX BAD TRIANGLES =====
% Only for methods including neighbor distance
if ismember(lower(Method), {'path', 'average'})
    % Configurations to detect:
    %    - One face divided in 3 with a point in the middle of the face
    %    - Square divided into 4 triangles with one point in the middle
    % Calculate face-vertex connectivity
    VertFacesConn = tess_faceconn(Faces);
    % Find vertices connected to three or four faces
    sumVert  = sum(VertConn,2);
    sumFaces = sum(VertFacesConn,2);
    % Three/Four vertices: average the values of their neighbors
    iVert = find(((sumVert == 3) & (sumFaces == 3)) | ((sumVert == 4) & (sumFaces == 4)));
    AvgConn = bst_bsxfun(@rdivide, W * VertConn(:,iVert), sumVert(iVert)');
    W(iVert,:) = AvgConn';
    W(:,iVert) = AvgConn;
end

author: brainstorm-tools


text:
function a = eq( x, y )

%Disciplined convex programming information for EQ (==):
%   Both the left- and right-hand sides of an equality constraint must
%   be affine (or constant). If either side of the constraint is complex,
%   then the real and imaginary portions are constrained separately.
%
%Disciplined geometric programming information for EQ (>):
%   Both the left- and right-hand sides of an equality constraint must
%   be log-affine, which includes positive constants and monomials.

evalin( 'caller', 'cvx_verify' );
b = cvx_pushcnstr( x, y, '==' );
if nargout, a = b; end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function pass = test_dst(pref)

if ( nargin == 0 )
    pref = chebfunpref();
end

% TODO: include more extensive tests.

n = 15;
seedRNG(42)
R = rand(n);

tol = 50*eps;

for k = 1:8
    pass(k) = norm(chebfun.idst(chebfun.dst(R,1),1) - R) < tol;
end

N = 15;
x = rand(N, 2);
tol = 100*eps;

% DST-I
DST1 = sin(pi/(N+1)*((0:N-1)+1)'*((0:N-1)+1));
y = DST1 * x;
pass(9) = norm( y - chebfun.dst(x, 1) ) < tol;

% DST-II
DST2 = sin(pi/N*((0:N-1)+1)'*((0:N-1)+1/2));
y = DST2 * x;
pass(10) = norm( y - chebfun.dst(x, 2) ) < tol;

% DST-III
DST3 = sin(pi/N*((0:N-1)+1/2)' * ((0:N-1)+1));
DST3(:, end) = .5*DST3(:, end);
y = DST3 * x;
pass(11) = norm( y - chebfun.dst(x, 3) ) < tol;

% DST-IV
DST4 = sin(pi/N*((0:N-1)+1/2)'*((0:N-1)+1/2));
y = DST4 * x;
pass(12) = norm( y - chebfun.dst(x, 4) ) < 10*tol;

% DST-V
DST5 = sin(pi/(N+1/2)*((0:N-1)+1)'*((0:N-1)+1));
y = DST5 * x;
pass(13) = norm( y - chebfun.dst(x, 5) ) < 10*tol;

% DST-VI
DST6 = sin(pi/(N+1/2)*((0:N-1)+1)'*((0:N-1)+1/2));
y = DST6 * x;
pass(14) = norm( y - chebfun.dst(x, 6) ) < 10*tol;

% DST-VII
DST7 = sin(pi/(N+1/2)*((0:N-1)+1/2)'*((0:N-1)+1));
y = DST7 * x;
pass(15) = norm( y - chebfun.dst(x, 7) ) < 10*tol;

% DST-VIII
DST8 = sin(pi/(N+1/2)*((0:N-1)+1/2)'*((0:N-1)+1/2));
y = DST8 * x;
pass(16) = norm( y - chebfun.dst(x, 8) ) < 10*tol;

% IDST-I
y = DST1 \ x;
pass(17) = norm( y - chebfun.idst(x, 1) ) < tol;

% IDST-II
y = DST2 \ x;
pass(18) = norm( y - chebfun.idst(x, 2) ) < tol;

% IDST-III
y = DST3 \ x;
pass(19) = norm( y - chebfun.idst(x, 3) ) < tol;

% IDST-IV
y = DST4 \ x;
pass(20) = norm( y - chebfun.idst(x, 4) ) < tol;

% IDST-V
y = DST5 \ x;
pass(21) = norm( y - chebfun.idst(x, 5) ) < tol;

% IDST-VI
y = DST6 \ x;
pass(22) = norm( y - chebfun.idst(x, 6) ) < tol;

% IDST-VII
y = DST7 \ x;
pass(23) = norm( y - chebfun.idst(x, 7) ) < tol;

% % IDST-VIII
pass(24) = true;
% y = DST8 \ x;
% pass(24) = norm( y - chebfun.idst(x, 8) ) < tol;

% Check complex inputs
c = rand(10,1) + 1i*rand(10,1);
for j = 1:8
    v1 = chebfun.dst(c, j);
    v2 = chebfun.dst(real(c),j) + 1i*chebfun.dst(imag(c),j);
    pass(24+j) = norm(v1 - v2) < tol;
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function automobile_plot ( )

%*****************************************************************************80
%
%% AUTOMOBILE_PLOT reads the automobile dataset and makes a scatterplot.
%
%  Discussion:
%
%    The hardest part of this exercise was reading the comma-separated data file,
%    which contains integer, real and string data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    A Frank, A Asuncion,
%    UCI Machine Learning Repository,
%    http://archive.ics.uci.edu/ml,
%    School of Information and Computer Science,
%    University of California, Irvine, California.
%

%
%  Read the data from the file.
%
  auto_cell = automobile_read ( 'automobile.txt' );
%
%  Price and weight are specific columns of the cell array.
%  Extract them, and convert them to numeric vectors.
%
  price = auto_cell(:,26);
  price = cell2mat ( price );
  weight = auto_cell(:,14);
  weight = cell2mat ( weight );
%
%  Ignore missing data.
%
  index = ( ( price ~= -1 ) & ( weight ~= -1 ) );
%
%  Create a scatter plot.
%
  figure ( 1 )
  clf

  plot ( price(index), weight(index), 'bo', 'MarkerSize', 5 )
  grid on
  xlabel ( 'Price in 1985 Dollars' );
  ylabel ( 'Curb weight in pounds' );
  title ( 'Scatter plot of price versus weight.' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function [BB2 Conf Valid tld] = tldTracking(tld,BB1,I,J)
% Estimates motion of bounding box BB1 from frame I to frame J

% initialize output variables
BB2    = []; % estimated bounding
Conf   = []; % confidence of prediction
Valid  = 0;  % is the predicted bounding box valid? if yes, learning will take place ...

if isempty(BB1) || ~bb_isdef(BB1), return; end % exit function if BB1 is not defined

% estimate BB2
xFI    = bb_points(BB1,10,10,5); % generate 10x10 grid of points within BB1 with margin 5 px
xFJ    = lk(2,tld.img{I}.input,tld.img{J}.input,xFI,xFI); % track all points by Lucas-Kanade tracker from frame I to frame J, estimate Forward-Backward error, and NCC for each point
medFB  = median2(xFJ(3,:)); % get median of Forward-Backward error
medNCC = median2(xFJ(4,:)); % get median for NCC
idxF   = xFJ(3,:) <= medFB & xFJ(4,:)>= medNCC; % get indexes of reliable points
BB2    = bb_predict(BB1,xFI(:,idxF),xFJ(1:2,idxF)); % estimate BB2 using the reliable points only

tld.xFJ = xFJ(:,idxF); % save selected points (only for display purposes)

% detect failures
if ~bb_isdef(BB2) || bb_isout(BB2,tld.imgsize), BB2 = []; return; end % bounding box out of image
if tld.control.maxbbox > 0 && medFB > 10, BB2 = []; return; end  % too unstable predictions

% estimate confidence and validity
patchJ   = tldGetPattern(tld.img{J},BB2,tld.model.patchsize); % sample patch in current image
[~,Conf] = tldNN(patchJ,tld); % estimate its Conservative Similarity (considering 50% of positive patches only)

% Validity
Valid    = tld.valid(I); % copy validity from previous frame
if Conf > tld.model.thr_nn_valid, Valid = 1; end % tracker is inside the 'core'

repo: OpenTLD


text:
function AdjustAxes(f,varargin)

%AdjustAxes - Adjust axes limits for all subplots
%
%  USAGE
%
%    AdjustAxes(f,<options>)
%
%    f              optional figure handle (default = gcf)
%    <options>      optional list of property-value pairs (see table below)
%
%    =========================================================================
%     Properties    Values
%    -------------------------------------------------------------------------
%     'x'           [min max], 'uniform' (same for all axes) or 'auto'
%     'y'           [min max], 'uniform' (same for all axes) or 'auto'
%     'xy'          [min max], 'uniform' (same for all axes) or 'auto'
%    =========================================================================
%
%  DEFAULTS
%
%    When no option is provided, all axes are automatically adjusted (same as
%    xlim('xy','uniform')). When 'x' is provided alone, y axes are left unchanged
%    and vice versa.
%

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

% Default values
xlims = [];
ylims = [];

% Optional parameter
if nargin < 1,
	f = gcf;
elseif ischar(f),
	varargin = {f,varargin{:}};
	f = gcf;
elseif ~ishandle(f),
	error('Incorrect value for ''f'' (type ''help AdjustAxes'' for details).');
end

% Parse parameter list
for i = 1:2:length(varargin),
  if ~ischar(varargin{i}),
    error(['Parameter ' num2str(i+1) ' is not a property (type ''help AdjustAxes'' for details).']);
  end
  switch(lower(varargin{i})),
    case 'x',
      xlims = lower(varargin{i+1});
      if ~isdvector(xlims,'<') && ~isstring_FMAT(xlims,'auto','uniform'),
        error('Incorrect value for property ''xlims'' (type ''help AdjustAxes'' for details).');
      end
    case 'y',
      ylims = lower(varargin{i+1});
      if ~isdvector(ylims,'<') && ~isstring_FMAT(ylims,'auto','uniform'),
        error('Incorrect value for property ''ylims'' (type ''help AdjustAxes'' for details).');
      end
    case 'xy',
      xlims = lower(varargin{i+1});
      ylims = lower(varargin{i+1});
      if ~isdvector(xlims,'<') && ~isstring_FMAT(xlims,'auto','uniform'),
        error('Incorrect value for property ''AdjustAxess'' (type ''help AdjustAxes'' for details).');
      end
    otherwise,
      error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help AdjustAxes'' for details).']);
  end
end

if isempty(xlims) && isempty(ylims),
	xlims = 'uniform';
	ylims = 'uniform';
end

% Compute min and max limits for all subplots
sub = get(f,'children');
if strcmp(xlims,'uniform'),
	lims = [];
	for i = 1:length(sub),
		if strcmp(get(sub,'type'),'axes'),
			lims(end+1,:) = xlim(sub(i));
		end
	end
	if isempty(lims), return; end
	xlims = [min(lims(:,1)) max(lims(:,2))];
end
if strcmp(ylims,'uniform'),
	lims = [];
	for i = 1:length(sub),
		if strcmp(get(sub,'type'),'axes'),
			lims(end+1,:) = ylim(sub(i));
		end
	end
	if isempty(lims), return; end
	ylims = [min(lims(:,1)) max(lims(:,2))];
end

% Apply
for i = 1:length(sub),
	if strcmp(get(sub,'type'),'axes'),
		if strcmp(xlims,'auto'),
			xlim(sub(i),'auto');
		elseif ~isempty(xlims),
			xlim(sub(i),xlims);
		end
		if strcmp(ylims,'auto'),
			ylim(sub(i),'auto');
		elseif ~isempty(ylims),
			ylim(sub(i),ylims);
		end
	end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef kMeans < handle
%
% k-means++: The Advantages of Careful Seeding", by David Arthur and
% Sergei Vassilvitskii, SODA 2007.

  properties
    n      % number of clusters
    center % cluster center
  end

  methods

    function km = kMeans(n)
      if nargin == 1, km.n = n; end
    end

    function initSeeds(km,obj)

      % Choose one center uniformly at random from among the data points.
      km.center = obj(randi(length(obj)));
      D = inf(length(obj),1);

      for k = 2:km.n

        %For each data point x, compute D(x), the distance between x and
        %the nearest center that has already been chosen.
        D = min(D,angle(obj,km.center(k-1)));

        %Choose one new data point at random as a new center, using a
        %weighted probability distribution where a point x is chosen with
        %probability proportional to D(x)2.
        id = discretesample(D.^2,1);
        km.center(k) = obj(id);

      end

    end

    function [label, center] = doClustering(km,obj)

      km.initSeeds(obj);

      label = find(km.center,obj);
      label_old = nan(size(label));

      % computes new centers until the labels do not change anymore
      while any(label(:) ~= label_old(:))

        % compute new centers
        for i = 1:km.n
          km.center(i) = mean(obj(label == i));
        end

        % store old labels
        label_old = label;

        % assign new labels
        label = find(km.center,obj);

      end

      center = km.center;

    end

  end

end

author: mtex-toolbox

repo: mtex

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


text:
function igood = check_clusters(hid, ss, rmax)

nmax = max(hid);
igood = zeros(nmax,1);
for j = 1:nmax
    ix = hid==j;
    if sum(ix)>0
        igood(j) = is_good(ss(ix), rmax);
    end
end
end

function igood = is_good(ss, rmax)
dt = 1/1000;
fcontamination = 0.1;
[K, Qi, Q00, Q01, rir] = ccg(ss, ss, 500, dt); % % compute the auto-correlogram with 500 bins at 1ms bins
Q = min(Qi/(max(Q00, Q01))); % this is a measure of refractoriness
R = min(rir); % this is a second measure of refractoriness (kicks in for very low firing rates)
if Q<fcontamination && R<rmax
   igood = 1;
else
    igood = 0;
end

end

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function c = r8ge_mxm ( n, a, b )

%*****************************************************************************80
%
%% R8GE_MXM multiplies two R8GE matrices.
%
%  Discussion:
%
%    The R8GE storage format is used for a general M by N matrix.  A storage
%    space is made for each logical entry.  The two dimensional logical
%    array is mapped to a vector, in which storage is by columns.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrices.
%    N must be positive.
%
%    Input, real A(N,N), B(N,N), the R8GE factor matrices.
%
%    Output, real C(N,N), the R8GE product matrix.
%
  c(1:n,1:n) = a(1:n,1:n) * b(1:n,1:n);

  return
end


text:
% 2-D Transform
% The 2-D Fourier transform is useful for processing 2-D signals and other 2-D data such as images.
% Create and plot 2-D data with repeated blocks.

P = peaks(20);
X = repmat(P,[5 10]);
imagesc(X)

% TODO : Compute the 2-D Fourier transform of the data.
% Shift the zero-frequency component to the center of the output, and
% plot the resulting 100-by-200 matrix, which is the same size as X.

% In the case of Radar signal processing. Convert the signal in MxN matrix, where M is the size of Range FFT samples and N is the size of Doppler FFT samples:
signal  = reshape(X, [100, 200]);

% Run the 2D FFT across both the dimensions.
signal_fft = fft2(signal, 100, 200);

% Shift zero-frequency terms to the center of the array
signal_fft = fftshift(signal_fft);

% Take the absolute value
signal_fft = abs(signal_fft);

% Here since it is a 2D output, it can be plotted as an image. Hence, we use the imagesc function
imagesc(signal_fft);

author: williamhyin


text:
function t = fieldnames(obj)
% Fieldnames of a file-array object
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: fieldnames.m 1143 2008-02-07 19:33:33Z spm $

t = {...
    'fname'
    'dim'
    'dtype'
    'offset'
    'scl_slope'
    'scl_inter'
};

author: fieldtrip

repo: fieldtrip


text:
% OtTableLast Move to the last record of the table.
%
% Usage:
%    OtTableLast(htable)

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-04-11 Created (Johan Meijdam, VORtech)

function iret = OtTableLast(htable)
   iret = Ot('table.last', htable);
end

author: Sable

repo: mcbench-benchmarks


text:
function [bel, niter] = parallel_protocol(engine, evidence, pot_type, local_kernel, msg)

bnet = bnet_from_engine(engine);
ns = bnet.node_sizes;
onodes = find(~isemptycell(evidence));

ndoms = length(engine.gdl.doms);
prod_of_msg = cell(1, ndoms);
bel = cell(1, ndoms);
old_bel = cell(1, ndoms);

converged = 0;
iter = 1;
while ~converged & (iter <= engine.max_iter)

  % each node multiplies all its incoming msgs and computes its local belief
  old_bel = bel;
  for i=1:ndoms
    prod_of_msg{i} = mk_initial_pot(pot_type, engine.gdl.doms{i}, ns, bnet.cnodes, onodes);
    nbrs = engine.gdl.nbrs{i};
    for j=nbrs(:)'
      prod_of_msg{i} = multiply_by_pot(prod_of_msg{i}, msg{j,i});
    end
    bel{i} = normalize_pot(multiply_by_pot(local_kernel{i}, prod_of_msg{i}));
  end

  if ~isempty(engine.fid)
    for i=1:ndoms
      tmp = pot_to_marginal(bel{i});
      %fprintf(engine.fid, '%9.7f ', tmp.T(1));
      fprintf(engine.fid, '%9.7f ', tmp.U(1));
    end
    %fprintf(engine.fid, '  U ');
    %for i=1:ndoms
    %  tmp = pot_to_marginal(bel{i});
    %  fprintf(engine.fid, '%9.7f ', tmp.U(1));
    %end
    fprintf(engine.fid, '\n');
  end

  % converged?
  if iter==1
    converged = 0;
  else
    converged = 1;
    for i=1:ndoms
      if ~approxeq_pot(bel{i}, old_bel{i}, engine.tol)
	converged = 0;
	break;
      end
    end
  end

  if ~converged
    old_msg = msg;
    % each node sends a msg to each of its neighbors
    for i=1:ndoms
      nbrs = engine.gdl.nbrs{i};
      for j=nbrs(:)'
	% multiply all incoming msgs except from j
	temp = prod_of_msg{i};
	temp = divide_by_pot(temp, old_msg{j,i});
	% send msg from i to j
	temp = multiply_by_pot(temp, local_kernel{i});
	temp2 = marginalize_pot(temp, engine.gdl.sepset{i,j}, engine.maximize);
	msg{i,j} = normalize_pot(temp2);
      end
    end
  end

  iter = iter + 1;
end

niter = iter-1;

if 0
for i=1:ndoms
  prod_of_msg{i} = mk_initial_pot(pot_type, engine.gdl.doms{i}, ns, bnet.cnodes, onodes);
  nbrs = engine.gdl.nbrs{i};
  for j=nbrs(:)'
    prod_of_msg{i} = multiply_by_pot(prod_of_msg{i}, msg{j,i});
  end
  bel{i} = normalize_pot(multiply_by_pot(local_kernel{i}, prod_of_msg{i}));
end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [ dim_num, node_num, element_num, element_order, node_data_num, ...
  node_coord, element_node, node_data ] = tec_read ( tec_file_name )

%*****************************************************************************80
%
%% TEC_READ reads finite element data from a TEC file.
%
%  Discussion:
%
%    This program reads a TEC file containing finite element data,
%    and writes that data out to three files that constitute an FEM model,
%    that is,
%    * a file of node coordinates;
%    * a file of elements defined by the nodes that form them;
%    * a file of node data.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character TEC_FILE_NAME(*), the name of the TEC file.
%
%    Output, integer DIM_NUM, the spatial dimension, inferred from the
%    names of the variables.
%
%    Output, integer NODE_NUM, the number of nodes, determined by the
%    "N=" argument.
%
%    Output, integer ELEMENT_NUM, the number of elements, inferred from the
%    "E=" argument.
%
%    Output, integer ELEMENT_ORDER, the order of the elements, inferred from
%    the "ZONETYPE=" argument.
%
%    Output, integer NODE_DATA_NUM, the number of data items per node,
%    inferred from the the number of node data items, minus those which are
%    inferred to be spatial coordinates.
%
%    Output, real NODE_COORD(DIM_NUM,NODE_NUM), the coordinates of nodes.
%
%    Output, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM);
%    the global index of local node I in element J.
%
%    Output, real NODE_DATA(NODE_DATA_NUM,NODE_NUM), the data values associated
%    with each node.
%
  dim_num = -1;
  node_num = -1;
  element_num = -1;
  element_order = -1;
  node_data_num = -1;
  node_coord = [];
  element_node = [];
  node_data = [];

  tec_file_unit = fopen ( tec_file_name );

  if ( tec_file_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEC_READ - Error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', tec_file_name );
    error ( 'TEC_READ - Error!' );
    return;
  end
%
%  Read and parse the TITLE line.
%  But it is optional, so you may have just read the VARIABLES line instead!
%
  line = ' ';

  while ( s_len_trim ( line ) == 0 )
    line = fgetl ( tec_file_unit );
  end
%
%  Read the VARIABLES line.
%
%  Because the TITLE line is apparently optional, we may have already
%  read the VARIABLES line!
%
  if ( s_begin ( line, 'TITLE=' ) )
    line = ' ';
    while ( s_len_trim ( line ) == 0 )
      line = fgetl ( tec_file_unit );
    end
  end

  if ( ~s_begin ( line, 'VARIABLES=' ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEC_READ - Fatal error!\n' );
    fprintf ( 1, '  The VARIABLES = line is missing in the file.\n' );
    error ( 'TEC_READ - Fatal error!' );
  end
%
%  Parse the VARIABLES line.
%  VARIABLES = name1 name2 name3...
%  The names may be quoted, and are separated by quotes, commas or spaces.
%
  [ variable_num, variable_name_length, variable_name ] ...
    = tec_variable_line_parse ( line );
%
%  Based on the variable names, determine the spatial dimension and the number
%  of node data items.
%
%  For now, we SIMPLY ASSUME that the spatial coordinates are listed first.
%  Hence, when we read the node data, we assume that the first DIM_NUM values
%  represent X, Y and possibly Z.
%
  dim_num = 0;
  node_data_num = variable_num;

  begin = 0;
  for variable = 1 : variable_num
    if ( variable_name_length(variable) == 1 )
      name = variable_name(begin+1);
      if ( ch_eqi ( name, 'X' ) || ...
          ch_eqi ( name, 'Y' ) || ...
          ch_eqi ( name, 'Z' ) )
        dim_num = dim_num + 1;
        node_data_num = node_data_num - 1;
      end
    end
    begin = begin + variable_name_length(variable);
  end
%
%  Read and parse the ZONE line.
%
  line = ' ';
  while ( s_len_trim ( line ) == 0 )
    line = fgetl ( tec_file_unit );
  end

  if ( ~s_begin ( line, 'ZONE' ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEC_READ - Fatal error!\n' );
    fprintf ( 1, '  The ZONE = line is missing in the file.\n' );
    error ( 'TEC_READ - Fatal error!' );
  end

  [ node_num, element_num, element_type ] = tec_zone_line_parse ( line );
%
%  Based on ELEMENT_TYPE, determine the element order.
%
  if ( s_eqi ( element_type, 'FETRIANGLE' ) )
    element_order = 3;
  elseif ( s_eqi ( element_type, 'FEQUADRILATERAL' ) )
    element_order = 4;
  elseif ( s_eqi ( element_type, 'FETETRAHEDRON' ) )
    element_order = 4;
  elseif ( s_eqi ( element_type, 'FEBRICK' ) )
    element_order = 8;
  else
    element_order = -1;
  end
%
%  Build up the format string for reading DIM_NUM + NODE_DATA_NUM reals.
%
  format = ' ';

  for i = 1 : dim_num + node_data_num
    format = strcat ( format, ' %f' );
  end

  node_coord = zeros ( dim_num, node_num );
  node_data = zeros ( node_data_num, node_num );
%
%  Now read the node coordinates and node data.
%
  for node = 1 : node_num

    line = ' ';
    while ( s_len_trim ( line ) == 0 )
      line = fgetl ( tec_file_unit );
    end

    [ x, count ] = sscanf ( line, format );

    if ( count == dim_num + node_data_num )
      node_coord(1:dim_num,node) = x(1:dim_num);

      node_data(1:node_data_num,node) = x(dim_num+1:dim_num+node_data_num);
    end

  end
%
%  Build up the format string for reading ELEMENT_ORDER integers.
%
  format = ' ';

  for element = 1 : element_order
    format = strcat ( format, ' %d' );
  end

  element_node = zeros ( element_order, element_num );
%
%  Now read the element data.
%
  for element = 1 : element_num

    line = ' ';
    while ( s_len_trim ( line ) == 0 )
      line = fgetl ( tec_file_unit );
    end

    [ x, count ] = sscanf ( line, format );

    if ( count == element_order )
      element_node(1:element_order,element) = x(1:element_order);
    end

  end

  fclose ( tec_file_unit );

  return
end