text:
classdef prtDecisionMap < prtDecision & prtActionBig
    % prtDecisionMap Maximum a-posteriori decision making
    %
    % prtDec = prtDecisionMap creates a prtDecisionBinaryMap
    % object, which can be used to perform Maximu a-posteriori decions.
    %
    % prtDecision objects are intended to be used either as members of
    % prtAlgorithm or prtClass objects.
    %
    % Example 1:
    %
    % ds = prtDataGenMary;                    % Load a data set
    % classifier = prtClassKnn;            % Create a clasifier
    % classifier = classifier.train(ds);   % Train the classifier
    % yOutClassifier = classifier.run(ds); % Run the classifier
    %
    % % Construct a prtAlgorithm object consisting of a prtClass object and
    % % a prtDecision object
    % algo = prtClassKnn + prtDecisionMap;
    %
    % algo = algo.train(ds);        % Train the algorithm
    % yOutAlgorithm = algo.run(ds); % Run the algorithm
    %
    % % Plot and compare the results
    % subplot(2,1,1); stem(yOutClassifier.getObservations); title('KNN Output');
    % subplot(2,1,2); stem(yOutAlgorithm.getObservations); title('KNN + Decision Output');
    %
    % Example 2:
    %
    % ds = prtDataGenMary;              % Load a data set
    % classifier = prtClassKnn;            % Create a clasifier
    % classifier = classifier.train(ds);   % Train the classifier
    %
    % % Plot the trained classifier
    % subplot(2,1,1); plot(classifier); title('KNN');
    %
    % % Set the classifiers internealDecider to be a prtDecsion object
    % classifier.internalDecider = prtDecisionMap;
    %
    % classifier = classifier.train(ds); % Train the classifier
    % subplot(2,1,2); plot(classifier); title('KNN + Decision');
    %
    % See also: prtDecisionBinary, prtDecisionBinarySpecifiedPd,
    % prtDecisionBinarySpecifiedPf, prtDecisionMap

     properties (SetAccess = private)
        name = 'MAP'  % MAP
        nameAbbreviation = 'MAP';  % MAP
     end

    properties (SetAccess = private, Hidden = true)
        runBinary = false;
        minPeDecision = [];
    end

    methods
        function self = prtDecisionMap(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end
    end
    methods (Access=protected,Hidden=true)

        function self = trainActionBig(self, dataSet)
            if isa(dataSet,'prtDataSetBigClass')
                self.classList = dataSet.uniqueClasses;
            else
                self.classList = 1:dataSet.nFeatures;
            end
            %binary classification
            if dataSet.nClasses == 2 && dataSet.nFeatures == 1
                warning('prtDecisionMap:binaryData','User specified MAP decision, but input data was binary, and classifier provided binary decision statistics.  MAP will default to prtDecisionBinaryMinPe, but there are subtle differences between these approaches.  This warning will only be shown once.  To turn it off permanently, use "warning off prtDecisionMap:binaryData" in your startup M-file');
                warning off prtDecisionMap:binaryData
                self.runBinary = true;
                self.minPeDecision = prtDecisionBinaryMinPe;
                self.minPeDecision = self.minPeDecision.trainBig(dataSet);
            end
        end

        function self = trainAction(self, dataSet)
            if isa(dataSet,'prtDataSetClass')
                self.classList = dataSet.uniqueClasses;
            else
                self.classList = 1:dataSet.nFeatures;
            end
            %binary classification
            if dataSet.nClasses == 2 && dataSet.nFeatures == 1
                warning('prtDecisionMap:binaryData','User specified MAP decision, but input data was binary, and classifier provided binary decision statistics.  MAP will default to prtDecisionBinaryMinPe, but there are subtle differences between these approaches.  This warning will only be shown once.  To turn it off permanently, use "warning off prtDecisionMap:binaryData" in your startup M-file');
                warning off prtDecisionMap:binaryData
                self.runBinary = true;
                self.minPeDecision = prtDecisionBinaryMinPe;
                self.minPeDecision = self.minPeDecision.train(dataSet);
            end
        end

        function dataSet = runAction(self,dataSet)
            yOut = dataSet.getObservations;

            %Under certain strange circumstances, e.g., in prtClassCascade,
            %the self.runBinary flag may be turned on, even when the actual
            %mode of operation should be "m-ary". It's not clear to me when
            %or why this is happening.  So I now check size(yOut,2) in
            %addition to self.runBinary...
            %       -Pete, 2011.11.21
            if size(yOut,2) == 1 && self.runBinary
                dataSet = self.minPeDecision.run(dataSet);
                return;
            else
                if size(yOut,2) > 1
                    [twiddle,index] = max(yOut,[],2); %#ok<ASGLU>
                else
                    error('prt:prtDecisionMap','Cannot run prtDecisionMap on algorithms with single-column output; use prtDecisionBinaryMinPe instead');
                end
            end
            if ~isempty(self.classList)
                if max(index(:)) <= length(self.classList)
                    classList = self.classList(index);
                else
                    classList = index;
                end
            else
                classList = index;
            end
            classList = classList(:);
            dataSet = dataSet.setObservations(classList);
        end

        function xOut = runActionFast(self,xIn,ds) %#ok<INUSD>
           [twiddle,index] = max(xIn,[],2); %#ok<ASGLU>
           xOut = self.classList(index);
           xOut = xOut(:);
        end
    end
end


text:
function y = tdis_prb(x,n)
% PURPOSE: calculates t-probabilities for elements in x-vector
%---------------------------------------------------
% USAGE: y = tdis_prb(x,n)
% where: x = vector containing computed t-values
%        n = degrees of freedom parameter
%---------------------------------------------------
% RETURNS:
%        y = a vector of marginal probability levels
% --------------------------------------------------
% SEE ALSO: fdis_prb(), chis_prb
%---------------------------------------------------

% written by:
% James P. LeSage, Dept of Economics
% University of Toledo
% 2801 W. Bancroft St,
% Toledo, OH 43606
% jpl@jpl.econ.utoledo.edu

if nargin ~= 2; error('Wrong # of arguments to tdis_prb'); end;
if n <=0; error('dof is negative or zero in tdis_prb'); end;

x2 = n./(n+x.^2);
one = find(x2 >= 1);
if length(one) > 0
    x2(one,1) = 1-1e-12;
end;
zip = find(x2 <= 0);
if length(zip) > 0
    x2(zip,1) = 1e-12;
end;

tmp = 1.0 - 0.5*betainc(x2,0.5*n,0.5);
y = 2*(1-tmp);

author: ambropo

repo: VAR-Toolbox

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


text:
% KM_DEMO_KCCA Demo file for kernel canonical correlation analysis
% algorithm.
%
% This script takes two multi-dimensional variables and uses kernel CCA to
% map them onto a single latent 1-D variable. The kernel matrices are
% decomposed using incomplete Cholesky decomposition in order to allow
% large data sets. This demo includes 3 flavors of the KCCA generalized
% eigenvalue problem, all yielding very similar results.
%
% Author: Steven Van Vaerenbergh (steven *at* gtas.dicom.unican.es), 2012.
%
% The algorithm in this file is based on the following publications:
% D. R. Hardoon, S. Szedmak and J. Shawe-Taylor, "Canonical Correlation
% Analysis: An Overview with Application to Learning Methods", Neural
% Computation, Volume 16 (12), Pages 2639--2664, 2004.
% F. R. Bach, M. I. Jordan, "Kernel Independent Component Analysis", Journal
% of Machine Learning Research, 3, 1-48, 2002.
%
% This file is part of the Kernel Methods Toolbox for MATLAB.
% https://github.com/steven2358/kmbox

close all; clear
rs = 1; % seed for random generator
rng('default')
rng(rs)

%% PARAMETERS
N = 1000;	% number of samples. method's complexity is O(NM^2)
Mmax = 50;  % max. M (number of components in incomplete Cholesky decomp.)
reg = 1E-5; % regularization
kerneltype = 'gauss';   % kernel type
kernelpar = 1;  % kernel parameter

%% PROGRAM
tic

% generate data
s = randn(N,1);	% latent signal
r1 = randn(N,1); r2 = randn(N,1);	% random (helper) variables

% option 1: two multi-dimensional variables that are mappable onto s
x1 = [tanh(r1-s)+0.1*r1 r1+3*s-1/10*(sin(3*s))];
x2 = [s - 2*(1-exp(-r2))./(1+exp(-r2)) r2.*s tanh(r2+s)];

% option 2: two invertible 1D nonlinear transformations
% x1 = tanh(0.8*s)+0.1*s;	% moderate saturation
% x2 = -1/10*(sin(s*3)+1.1*s*3);   % stairway

% clean up: remove mean
x1 = x1-repmat(mean(x1),N,1);
x2 = x2-repmat(mean(x2),N,1);

% normalize variance (to improve anisotropy) or any other preprocessing
x1 = x1*sqrt(diag(1./diag(x1'*x1)));
x2 = x2*sqrt(diag(1./diag(x2'*x2)));

% KCCA
[y1,y2,beta] = km_kcca(x1,x2,kerneltype,kernelpar,reg,1,'ICD',Mmax);

% scale the estimated signals to compare without the scalar ambiguity
scaling = sqrt(var(s))/sqrt(var(y1))*sign(s(1))*sign(y1(1));

% mean square errors
error1 = s-scaling*y1;
error2 = s-scaling*y2;
MSE1 = sum(error1.^2)/N;
MSE2 = sum(error2.^2)/N;

toc
%% OUTPUT

figure; hold all
plot(s)
plot(scaling*y1);
plot(scaling*y2);
legend('latent variable','projection 1','projection 2')

fprintf('2x %d data points\n',N)
fprintf('Canonical correlation: %f\n',beta)
fprintf('MSE1: %f\n',MSE1);
fprintf('MSE2: %f\n',MSE2);
fprintf('\n')

% figure;plot(sort(diag(real(betas))))  % check eigenvalues
% figure;plot(sort(diag(betas))); % imaginary part due to numerical error

author: epfl-lasa


text:
function [seriesNameArray] = dispFileSeries(inputFilePath,varargin)
	% [bfSeriesNo] = getSeriesNoFromName(inputFilePath,bfSeriesName,varargin)
	%
	% Returns series number for Bio-Formats series with name given by user.
	%
	% Biafra Ahanonu
	% started: 2022.07.15 [08:10:16]
	%
	% Inputs
	% 	inputFilePath - Str: path to Bio-Formats compatible file.
	%
	% Outputs
	% 	bfSeriesNo - Int: series number (1-based indexing) matching the input series name. NaN is output if no series is found.
	%
	% Options (input as Name-Value with Name = options.(Name))
	% 	% File ID connection from calling bfGetReader(inputFilePath), this is to save time on larger files by avoiding opening the connection again.
	% 	options.fileIdOpen = '';

	% Changelog
		%
	% TODO
		%

	% ========================
	% File ID connection from calling bfGetReader(inputFilePath), this is to save time on larger files by avoiding opening the connection again.
	options.fileIdOpen = [];
	% get options
	options = ciapkg.io.getOptions(options,varargin);
	% disp(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	% ========================

	try
		% By default output NaN if no series found.
		seriesNameArray = {};

		if isempty(options.fileIdOpen)
			disp(['Opening connection to Bio-Formats file:' inputFilePath])
			startReadTime = tic;
			fileIdOpen = bfGetReader(inputFilePath);
			toc(startReadTime)
		else
			fileIdOpen = options.fileIdOpen;
		end
		omeMeta = fileIdOpen.getMetadataStore();

		nSeries = fileIdOpen.getSeriesCount();
		seriesNameArray = cell([1 nSeries]);
		disp('===')
		disp(['Series in file: ' inputFilePath])
		for seriesNo = 1:nSeries
			% Convert 1-based to 0-based indexing.
			thisStr = char(omeMeta.getImageName(seriesNo-1));
			seriesNameArray{seriesNo} = thisStr;
			disp([num2str(seriesNo) ': "' thisStr '"'])
		end
		disp('===')
	catch err
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end

	function [outputs] = nestedfxn_exampleFxn(arg)
		% Always start nested functions with "nestedfxn_" prefix.
		% outputs = ;
	end
end
function [outputs] = localfxn_exampleFxn(arg)
	% Always start local functions with "localfxn_" prefix.
	% outputs = ;
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function printpdf(fname)
%PRINTPDF Prints the current figure into a pdf document
set(gca, 'LooseInset', get(gca, 'TightInset'));
fname = [regexprep(fname, '^(.*)\.pdf$', '$1'), '.eps'];
% fname = [regexprep(fname, '^(.*)\.pdf$', '$1'), '.pdf'];
print('-depsc', fname) ;
% print('-dpdf', fname);
if ~system(['epstopdf ', fname])
% if ~system(['/usr/bin/epstopdf ', fname])
system(['rm ', fname]);
end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function zi = lagrange_interp_nd_value ( m, n_1d, a, b, nd, zd, ni, xi )

%*****************************************************************************80
%
%% LAGRANGE_INTERP_ND_VALUE evaluates an ND Lagrange interpolant.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N_1D(M), the order of the 1D rule to be used
%    in each dimension.
%
%    Input, real A(M), B(M), the lower and upper interval endpoints in
%    each dimension.
%
%    Input, integer ND, the number of points in the product grid.
%
%    Input, real ZD(ND), the function evaluated at the points XD.
%
%    Input, integer NI, the number of points at which the interpolant
%    is to be evaluated.
%
%    Input, real XI(M,NI), the points at which the interpolant is to be
%    evaluated.
%
%    Output, real ZI(NI), the interpolant evaluated at the points XI.
%
  zi = zeros(ni,1);

  for j = 1 : ni

    w = ones ( nd, 1 );

    for i = 1 : m
      n = n_1d(i);
      x_1d(1:n) = cc_compute_points ( n );
      x_1d(1:n) = 0.5 * ( ( 1.0 - x_1d(1:n) ) * a(i) ...
                        + ( 1.0 + x_1d(1:n) ) * b(i) );
      value = lagrange_basis_1d ( n, x_1d, 1, xi(i,j) );
      w = r8vec_direct_product2 ( i, n, value, m, nd, w );
    end

    zi(j) = zd' * w;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function drawAffine(afnv, tsize, color, linewidth)

rect= round(aff2image(afnv', tsize));
inp	= reshape(rect,2,4);

topleft_r = inp(1,1);
topleft_c = inp(2,1);
botleft_r = inp(1,2);
botleft_c = inp(2,2);
topright_r = inp(1,3);
topright_c = inp(2,3);
botright_r = inp(1,4);
botright_c = inp(2,4);
p = line([topleft_c, topright_c], [topleft_r, topright_r]);
set(p, 'Color', color); set(p, 'LineWidth', linewidth); set(p, 'LineStyle', '-');
p = line([topright_c, botright_c], [topright_r, botright_r]);
set(p, 'Color', color); set(p, 'LineWidth', linewidth); set(p, 'LineStyle', '-');
p = line([botright_c, botleft_c], [botright_r, botleft_r]);
set(p, 'Color', color); set(p, 'LineWidth', linewidth); set(p, 'LineStyle', '-');
p = line([botleft_c, topleft_c], [botleft_r, topleft_r]);
set(p, 'Color', color); set(p, 'LineWidth', linewidth); set(p, 'LineStyle', '-');

author: yuxng


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MRiLab auto generated file: DO NOT EDIT!     %
% Generated by MRiLab "DoWriteXML2m" Generator %
% MRiLab Version 1.3                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [B1x, B1y, B1z, E1x, E1y, E1z, Pos]=Coil_9ChSurfChest
%====================================
p.Azimuth=pi/2;
p.CoilID=1;
p.CurrentDir=1;
p.Elevation=0;
p.PosX=-0.1;
p.PosY=-0.1;
p.PosZ=-0.1;
p.Radius=0.04;
p.Scale=0.1;
p.Segment=20;
[B1xt, B1yt, B1zt, E1xt, E1yt, E1zt, Post]=CoilCircle(p);
B1x(:,:,:,p.CoilID)=B1xt;
B1y(:,:,:,p.CoilID)=B1yt;
B1z(:,:,:,p.CoilID)=B1zt;
E1x(:,:,:,p.CoilID)=E1xt;
E1y(:,:,:,p.CoilID)=E1yt;
E1z(:,:,:,p.CoilID)=E1zt;
Pos(p.CoilID,:)=Post;
p=[];
%--------------------
p.Azimuth=pi/2;
p.CoilID=2;
p.CurrentDir=1;
p.Elevation=0;
p.PosX=-0.1;
p.PosY=-0.1;
p.PosZ=0;
p.Radius=0.04;
p.Scale=0.1;
p.Segment=20;
[B1xt, B1yt, B1zt, E1xt, E1yt, E1zt, Post]=CoilCircle(p);
B1x(:,:,:,p.CoilID)=B1xt;
B1y(:,:,:,p.CoilID)=B1yt;
B1z(:,:,:,p.CoilID)=B1zt;
E1x(:,:,:,p.CoilID)=E1xt;
E1y(:,:,:,p.CoilID)=E1yt;
E1z(:,:,:,p.CoilID)=E1zt;
Pos(p.CoilID,:)=Post;
p=[];
%--------------------
p.Azimuth=pi/2;
p.CoilID=3;
p.CurrentDir=1;
p.Elevation=0;
p.PosX=-0.1;
p.PosY=-0.1;
p.PosZ=0.1;
p.Radius=0.04;
p.Scale=0.1;
p.Segment=20;
[B1xt, B1yt, B1zt, E1xt, E1yt, E1zt, Post]=CoilCircle(p);
B1x(:,:,:,p.CoilID)=B1xt;
B1y(:,:,:,p.CoilID)=B1yt;
B1z(:,:,:,p.CoilID)=B1zt;
E1x(:,:,:,p.CoilID)=E1xt;
E1y(:,:,:,p.CoilID)=E1yt;
E1z(:,:,:,p.CoilID)=E1zt;
Pos(p.CoilID,:)=Post;
p=[];
%--------------------
p.Azimuth=pi/2;
p.CoilID=4;
p.CurrentDir=1;
p.Elevation=0;
p.PosX=0;
p.PosY=-0.1;
p.PosZ=-0.1;
p.Radius=0.04;
p.Scale=0.1;
p.Segment=20;
[B1xt, B1yt, B1zt, E1xt, E1yt, E1zt, Post]=CoilCircle(p);
B1x(:,:,:,p.CoilID)=B1xt;
B1y(:,:,:,p.CoilID)=B1yt;
B1z(:,:,:,p.CoilID)=B1zt;
E1x(:,:,:,p.CoilID)=E1xt;
E1y(:,:,:,p.CoilID)=E1yt;
E1z(:,:,:,p.CoilID)=E1zt;
Pos(p.CoilID,:)=Post;
p=[];
%--------------------
p.Azimuth=pi/2;
p.CoilID=5;
p.CurrentDir=1;
p.Elevation=0;
p.PosX=0;
p.PosY=-0.1;
p.PosZ=0;
p.Radius=0.04;
p.Scale=0.1;
p.Segment=20;
[B1xt, B1yt, B1zt, E1xt, E1yt, E1zt, Post]=CoilCircle(p);
B1x(:,:,:,p.CoilID)=B1xt;
B1y(:,:,:,p.CoilID)=B1yt;
B1z(:,:,:,p.CoilID)=B1zt;
E1x(:,:,:,p.CoilID)=E1xt;
E1y(:,:,:,p.CoilID)=E1yt;
E1z(:,:,:,p.CoilID)=E1zt;
Pos(p.CoilID,:)=Post;
p=[];
%--------------------
p.Azimuth=pi/2;
p.CoilID=6;
p.CurrentDir=1;
p.Elevation=0;
p.PosX=0;
p.PosY=-0.1;
p.PosZ=0.1;
p.Radius=0.04;
p.Scale=0.1;
p.Segment=20;
[B1xt, B1yt, B1zt, E1xt, E1yt, E1zt, Post]=CoilCircle(p);
B1x(:,:,:,p.CoilID)=B1xt;
B1y(:,:,:,p.CoilID)=B1yt;
B1z(:,:,:,p.CoilID)=B1zt;
E1x(:,:,:,p.CoilID)=E1xt;
E1y(:,:,:,p.CoilID)=E1yt;
E1z(:,:,:,p.CoilID)=E1zt;
Pos(p.CoilID,:)=Post;
p=[];
%--------------------
p.Azimuth=pi/2;
p.CoilID=7;
p.CurrentDir=1;
p.Elevation=0;
p.PosX=0.1;
p.PosY=-0.1;
p.PosZ=-0.1;
p.Radius=0.04;
p.Scale=0.1;
p.Segment=20;
[B1xt, B1yt, B1zt, E1xt, E1yt, E1zt, Post]=CoilCircle(p);
B1x(:,:,:,p.CoilID)=B1xt;
B1y(:,:,:,p.CoilID)=B1yt;
B1z(:,:,:,p.CoilID)=B1zt;
E1x(:,:,:,p.CoilID)=E1xt;
E1y(:,:,:,p.CoilID)=E1yt;
E1z(:,:,:,p.CoilID)=E1zt;
Pos(p.CoilID,:)=Post;
p=[];
%--------------------
p.Azimuth=pi/2;
p.CoilID=8;
p.CurrentDir=1;
p.Elevation=0;
p.PosX=0.1;
p.PosY=-0.1;
p.PosZ=0;
p.Radius=0.04;
p.Scale=0.1;
p.Segment=20;
[B1xt, B1yt, B1zt, E1xt, E1yt, E1zt, Post]=CoilCircle(p);
B1x(:,:,:,p.CoilID)=B1xt;
B1y(:,:,:,p.CoilID)=B1yt;
B1z(:,:,:,p.CoilID)=B1zt;
E1x(:,:,:,p.CoilID)=E1xt;
E1y(:,:,:,p.CoilID)=E1yt;
E1z(:,:,:,p.CoilID)=E1zt;
Pos(p.CoilID,:)=Post;
p=[];
%--------------------
p.Azimuth=pi/2;
p.CoilID=9;
p.CurrentDir=1;
p.Elevation=0;
p.PosX=0.1;
p.PosY=-0.1;
p.PosZ=0.1;
p.Radius=0.04;
p.Scale=0.1;
p.Segment=20;
[B1xt, B1yt, B1zt, E1xt, E1yt, E1zt, Post]=CoilCircle(p);
B1x(:,:,:,p.CoilID)=B1xt;
B1y(:,:,:,p.CoilID)=B1yt;
B1z(:,:,:,p.CoilID)=B1zt;
E1x(:,:,:,p.CoilID)=E1xt;
E1y(:,:,:,p.CoilID)=E1yt;
E1z(:,:,:,p.CoilID)=E1zt;
Pos(p.CoilID,:)=Post;
p=[];
%--------------------
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function test_conjgradmv()

N = 100;
x = rand(N,1);
A = sparse(double(rand(N) < 0.1));
K = A*A' + speye(N);

ind = 1:2:N;

xmv = x(ind);
Kmv = K(ind,ind);
ymv = Kmv*xmv;
y = zeros(N,1);
y(ind) = ymv;
I = false(N,1);
I(ind) = true;

x_pcg = pcg(Kmv,ymv,1e-6,100);
x_cg = conjgradmv(K,y,I,'tol',1e-6,'maxiter',100,'verbose',true);

norm(x_pcg-x_cg(ind))

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function test_deep_nmf(varargin)
%
% demonstration file for NMFLibrary.
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai on Apr. 05, 2017

    if nargin < 1
        clc;
        clear;
        close all;
        rng('default')

        m = 300;
        n = 500;
        V = rand(m,n);
        options = [];
        options.verbose = 1;
        options.max_epoch = 100;
        health_check_mode = false;
    else
        V = varargin{1};
        %rank = varargin{2};
        options = varargin{3};
        health_check_mode = true;
    end

    %% Initialize of rank to be factorized
    rank_layers = [49 25 16];

    %% Deep-Semi-NMF
    [w_deep_semi_nmf, infos_deep_semi_nmf] = deep_semi_nmf(V, rank_layers, options);

    %% Deep-ns-NMF
    options.theta = 0.5;
    %options.update_alg = 'mu';
    options.update_alg = 'apg';
    options.apg_maxiter = 10;
    [w_deep_ns_nmf, infos_deep_ns_nmf] = deep_ns_nmf(V, rank_layers, options);

    [w_deep_bi_nmf, infos_deep_bi_nmf] = deep_bidirectional_nmf(V, rank_layers, options);

    %% Plotting
    if ~health_check_mode
        display_graph('iter','cost', {'Deep-SemiNMF', 'Deep-nsNMF', 'Deep-Bidir-SemiNMF'}, {w_deep_semi_nmf, w_deep_ns_nmf, w_deep_bi_nmf}, {infos_deep_semi_nmf, infos_deep_ns_nmf, infos_deep_bi_nmf});
    end

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [A]=patchArea(F,V)

% function [A]=patchArea(F,V)
% ------------------------------------------------------------------------
% This simple function calculates the areas of the faces specified by F and
% V. The output is a vector A containing size(F,1) elements. The face areas
% are calculated via triangulation of the faces. If faces are already
% triangular triangulation is skipped are area calculation is direction
% performed.
%
%
%
% Kevin Mattheus Moerman
%
% 2011/04/12
% 2021/09/13 Updated to use more efficient patchEdgeCrossProduct method
% 2021/09/14 Renamed to patchArea
%------------------------------------------------------------------------

%%

C=patchEdgeCrossProduct(F,V);
A=sqrt(sum(C.^2,2));

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function value = sphere_unit_monomial_nd ( n, p )

%*****************************************************************************80
%
%% SPHERE_UNIT_MONOMIAL_ND integrates a monomial on a unit sphere in ND.
%
%  Integration region:
%
%    Points X(1:N) such that
%
%      Sum ( X(1:N)**2 ) == 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gerald Folland,
%    How to Integrate a Polynomial Over a Sphere,
%    American Mathematical Monthly,
%    Volume 108, May 2001, pages 446-448.
%
%  Parameters:
%
%    Input, integer N, the dimension of the space.
%
%    Input, integer P(N), the exponents of X(1) through X(N) in the monomial.
%    The exponents P(N) must be nonnegative.
%
%    Output, real SPHERE_UNIT_MONOMIAL_ND, the integral of
%    X1**P(1)*X2**P(2)*...*XN**P(N) over the unit sphere.
%
  for i = 1 : n
    if ( mod ( p(i), 2 ) == 1 )
      value = 0.0E+00;
      return
    end
  end

  temp = 0.0E+00;
  arg2 = 0.0E+00;

  for i = 1 : n
    arg1 = ( p(i) + 1 ) / 2.0E+00;
    temp = temp + gammaln ( arg1 );
    arg2 = arg2 + arg1;
  end
  temp = temp - gammaln ( arg2 );

  value = 2.0E+00 * exp ( temp );

  return
end


text:
function [f] = binomial(n,d)
%BINOMIAL  calculate the binomial coefficient
%          n and d may be complex and any equal size
%
%        n!
% f = --------
%     d!(n-d)!
%
%usage: b = binomial(n,d)
%
%tested on version 5.3.1
%
%see also: Gamma, Prod
%see also: mhelp binomial
%
%not correct for negative integer arguments!

%Paul Godfrey
%pgodfrey@conexant.com
%8-12-00

nsize=size(n); n=n(:);
dsize=size(d); d=d(:);

if min(nsize==dsize)==0
   error('Input argument size mismatch!')
end

if (0<n & 0<=d & round(n)==n & round(d)==d & n>=d & max(size(n))==1 & max(size(d))==1)
   f = prod((d+1):n)/prod(1:(n-d));
else
   z=gamma([n d n-d]+1);
   f=z(:,1)./(z(:,2).*z(:,3));
end

f=reshape(f,nsize);

return

%a demo of this routine is
dx=1/16;
dy=dx;
x=-4:dx:4;
y=-4:dy:4;
[X,Y]=meshgrid(x,y);
f=binomial(X,Y);
p=find(abs(f)>4);
%f(p)=sign(real(f(p)))*4;
f(p)=NaN;
mesh(x,y,real(f))
view([23 54]);
rotate3d

author: Sable

repo: mcbench-benchmarks


text:
%% This function will generate the ODE function based on calculation result of the sybolic expression
% Last Update: 2019/04/22
% Coded By: K

function Generate_ODE_RHS(Sindy_ODEs,var_num_state,var_num_control)
z=sym('z',[var_num_state,1]);
u=sym('u',[var_num_control,1]);
syms t
f= matlabFunction(Sindy_ODEs,'File','Sindy_ODE_RHS','Optimize',true,'Vars',{t,z,u});

author: dynamicslab

repo: SINDy-PI


text:
function value = monomial_int01 ( dim_num, expon )

%*****************************************************************************80
%
%% MONOMIAL_INT01 returns the integral of a monomial over the [0,1] hypercube.
%
%  Discussion:
%
%    This routine evaluates a monomial of the form
%
%      product ( 1 <= dim <= dim_num ) x(dim)^expon(dim)
%
%    where the exponents are nonnegative integers.  Note that
%    if the combination 0^0 is encountered, it should be treated
%    as 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer EXPON(DIM_NUM), the exponents.
%
%    Output, real VALUE, the value of the integral of the
%    monomial over the [0,1] hypercube.
%
  value = 1.0 / prod ( expon(1:dim_num) + 1 );

  return
end


text:
function [movieInfo] = getMovieInfo(inputMovie,varargin)
	% Wrapper to quickly get movie information. Options are the same as loadMovieList.
	% Biafra Ahanonu
	% started: 2020.09.29 [13:23:57]
	% inputs
		% inputMovie - [x y frames] or char string (path to movie).
	% outputs
		%

	% changelog
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
        % 2022.01.20 [22:50:25] - Directly call loadMovieList and bypass ciapkg.api.
	% TODO
		%

	% import ciapkg.api.* % import CIAtah functions in ciapkg package API.

	% ========================
	% DESCRIPTION
	% OPTIONS ARE THE SAME AS loadMovieList.
	% options.exampleOption = '';
	% get options
	% options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	% ========================

	try
		movieInfo = [];
		movieInfo = ciapkg.io.loadMovieList(inputMovie,'getMovieDims',1,'passArgs',varargin);
	catch err
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end
end


text:
function [U, Uc] = slrangespace(tar, varargin)
%SLRANGESPACE Determines the subspace of the range of X
%
% $ Syntax $
%   - U = slrangespace(tar)
%   - U = slrangespace(tar, ...)
%   - [U, Uc] = slrangespace(...)
%
% $ Arguments $
%   - tar:      the target, can be a sample matrix or covariance
%   - U:        the basis for the range space
%   - Uc:       the basis for the orthogonal complement of the range
%
% $ Description $
%   - U = slrangespace(tar) determines the range space of tar in default
%     settings. That is to set the dimension to be the rank of tar. Note
%     that tar can have two forms: a sample matrix or a covariance
%     given by the syntax {'cov', C}.
%
%   - U = slrangespace(tar, ...) determines the range space of X using
%     the specified dimension determination schemes. The arguments
%     input following X will be delivered to sldim_by_eigval for dimension
%     determination.
%
%   - [U, Uc] = slrangespace(...) also returns orthogonal complement of U.
%
% $ History $
%   - Created by Dahua Lin on Apr 25, 2005
%

%% parse and verify input arguments

if isnumeric(tar)
    target_type = 1;        % target is sample

    if ndims(tar) ~= 2
        error('sltoolbox:invalidarg', ...
            'When tar is sample matrix, it should be a 2D matrix');
    end

    X = tar;
    d = size(X, 1);

elseif iscell(tar)

    if length(tar) == 2 && ischar(tar{1}) ...
            && strcmpi(tar{1}, 'cov')

        target_type = 2;        % target is covariance

        C = tar{2};
        d = size(C, 1);

        if ~isequal(size(C), [d d])
            error('sltoolbox:invalidarg', ...
                'A covariance matrix should be square');
        end

    else
        error('sltoolbox:invalidarg', ...
            'The target should be a sample matrix or a covariance given by cell');
    end

else

    error('sltoolbox:invalidarg', ...
        'The target should be a sample matrix or a covariance given by cell');
end

if nargout == 0
    return;
elseif nargout == 1
    need_complement = false;
else
    need_complement = true;
end

%% compute

switch target_type

    case 1      % sample matrix

        if ~need_complement
            [U, D] = svd(X, 0);
        else
            [U, D] = svd(X);
        end

        evals = diag(D) .^ 2;
        clear D;

    case 2      % covariance

        [evals, U] = slsymeig(C);

end

evals = max(evals, 0);
k = sldim_by_eigval(evals, varargin{:});

%% output

if ~need_complement
    U = U(:, 1:k);
else
    if k == size(U, 2)
        Uc = zeros(d, 0);
    else
        Uc = U(:, k+1:end);
        U  = U(:, 1:k);
    end
end

author: lmthang


text:
function ll = gpsimLogLikelihood(model)

% GPSIMLOGLIKELIHOOD Compute the log likelihood of a GPSIM model.
% FORMAT
% DESC computes the log likelihood of the given Gaussian process
% for use in a single input motif protein network.
% ARG model : the model for which the log likelihood is computed.
% RETURN ll : the log likelihood of the data set.
%
% SEEALSO : gpsimCreate, gpsimLogLikeGradient, gpsimObjective
%
% COPYRIGHT : Neil D. Lawrence, 2006

% SHEFFIELDML

dim = size(model.y, 1);
ll = -dim*log(2*pi) - model.logDetK - model.m'*model.invK*model.m;
ll = ll*0.5;

% In case we need priors in.
if isfield(model, 'bprior'),
  ll = ll + kernPriorLogProb(model.kern);
  ll = ll + priorLogProb(model.bprior, model.B);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function Neff = getNeff(dens)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% getNeff(P)
%  returns an estimate of Neff, the 'effective' number of points in the NPDE
%   (= sum(1/weights^2) )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2003 Alexander Ihler; distributable under GPL -- see README.txt

  Neff = 1/sum(getWeights(dens).^2);

author: ShapeNet

repo: RenderForCNN


text:
function eeglabciplot(lower,upper,x,colour, alphaval)
% eeglabciplot(lower,upper)
% eeglabciplot(lower,upper,x)
% eeglabciplot(lower,upper,x,colour)
%
% Plots a shaded region on a graph between specified lower and upper confidence intervals (L and U).
% l and u must be vectors of the same length.
% Uses the 'fill' function, not 'area'. Therefore multiple shaded plots
% can be overlaid without a problem. Make them transparent for total visibility.
% x data can be specified, otherwise plots against index values.
% colour can be specified (eg 'k'). Defaults to blue.

% Author:   Raymond Reynolds 24/11/06
% Modified by Ramon Martinez Cancino

if length(lower)~=length(upper)
    error('lower and upper vectors must be same length')
end

if nargin<5
    alphaval = 1;
end
if nargin<4
    colour= 'r';
end

if nargin<3
    x=1:length(lower);
end

% convert to row vectors so fliplr can work
if find(size(x)==(max(size(x))))<2
x=x'; end
if find(size(lower)==(max(size(lower))))<2
lower=lower'; end
if find(size(upper)==(max(size(upper))))<2
upper=upper'; end

hdl_tmp = fill([x fliplr(x)],[upper fliplr(lower)],colour,'FaceAlpha', alphaval, 'EdgeColor', 'none');

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function M = slmetric_pw(X1, X2, mtype, varargin)
%SLMETRIC_PW Compute the metric between column vectors pairwisely
%
% [ Syntax ]
%   - M = slmetric_pw(X1, X2, mtype);
%   - M = slmetric_pw(X1, X2, mtype, ...);
%
% [ Arguments ]
%   - X1, X2:       the sample matrices
%   - mtype:        the string indicating the type of metric
%   - M:            the resulting metric matrix
%
% [ Description ]
%    - M = slmetric_pw(X1, X2, mtype) Computes the metrics between
%      column vectors of X1 and X2 pairwisely, using the metric
%      specified by mtype.
%
%      Both X1 and X2 are matrices with each column representing a
%      sample. X1 and X2 should have the same number of rows. Suppose
%      the size of X1 is d x n1, and the size of X2 is d x n2. Then
%      the output metric matrix M will be of size n1 x n2, in which
%      M(i, j) is the metric value between X1(:,i) and X2(:,j).
%
%    - M = slmetric_pw(X1, X2, mtype, ...) Some metric types requires
%      extra parameters, which should be specified in params.
%
%      The supported metrics of this function are listed as follows:
%      \{:
%        - eucdist:        Euclidean distance:
%                          $ ||x - y|| $
%
%        - sqdist:         Square of Euclidean distance:
%                          $ ||x - y||^2 $
%
%        - dotprod:        Canonical dot product:
%                          $ <x,y> = x^T * y $
%
%        - nrmcorr:        Normalized correlation (cosine angle):
%                          $ (x^T * y ) / (||x|| * ||y||) $
%
%        - corrdist:       Normalized Correlation distance
%                          $ 1 - nrmcorr(x, y) $
%
%        - angle:          Angle between two vectors (in radian)
%                          $ arccos (nrmcorr(x, y)) $
%        - quadfrm:        Quadratic form:
%                          $ x^T * Q * y $
%                         Q is specified in the 1st extra parameter
%
%        - quaddiff:       Quadratic form of difference:
%                          $ (x - y)^T * Q * (x - y) $
%                         Q is specified in the 1st extra parameter
%
%        - cityblk:        City block distance (abssum of difference)
%                          $ sum_i |x_i - y_i| $
%
%        - maxdiff:        Maximum absolute difference
%                          $ max_i |x_i - y_i| $
%
%        - mindiff:        Minimum absolute difference
%                          $ min_i |x_i - y_i| $
%
%        - minkowski:      Minkowski distance
%                          $ (\sum_i |x_i - y_i|^p)^(1/p) $
%                         The order p is specified in the 1st extra parameter
%
%        - wsqdist:        Weighted square of Euclidean distance
%                          $ \sum_i w_i (x_i - y_i)^2 $
%                         the weights w is specified in 1st extra parameter
%                         as a d x 1 column vector
%
%        - hamming:        Hamming distance with threshold t
%                          \{
%                              ht1 = x > t
%                              ht2 = y > t
%                              d = sum(ht1 ~= ht2)
%                          \}
%                         use threshold t as the first extra param.
%                         (by default, t is set to zero).
%
%        - hamming_nrm:    Normalized hamming distance, which equals the
%                          ratio of the elements that differ.
%                          \{
%                              ht1 = x > t
%                              ht2 = y > t
%                              d = sum(ht1 ~= ht2) / length(ht1)
%                          \}
%                          use threshold t as the first extra param.
%                         (by default, t is set to zero).
%
%        - intersect:      Histogram Intersection
%                           $ d = sum min(x, y) / min(sum(x), sum(y))$
%
%        - intersectdis:   Histogram intersection distance
%                           $ d = 1 - sum min(x, y) / min(sum(x), sum(y)) $
%
%        - chisq:          Chi-Square Distance
%                           $ d = sum (x(i) - y(i))^2/(2 * (x(i)+y(i))) $
%
%        - kldiv:          Kull-back Leibler divergence
%                           $ d = sum x(i) log (x(i) / y(i)) $
%
%        - jeffrey:        Jeffrey divergence
%                           $ d = KL(h1, (h1+h2)/2) + KL(h2, (h1+h2)/2) $
%      \:}
%
% [ Remarks ]
%   - Both X1 and X2 should be a matrix of numeric values, except
%     for case when metric type is 'hamming' or 'hamming_nrm'.
%     For hamming or hamming_nrm metric, the input matrix can be logical.
%
% [ Examples ]
%   - Compute different types of metrics in pairwise manner
%     \{
%         % prepare sample matrix
%         X1 = rand(10, 100);
%         X2 = rand(10, 150);
%
%         % compute the euclidean distances (L2)
%         % between the samples in X1 and X2
%         M = slmetric_pw(X1, X2, 'eucdist');
%
%         % compute the eucidean distances between the samples
%         % in X1 in a pairwise manner
%         M = slmetric_pw(X1, X1, 'eucdist');
%
%         % compute the city block distances (L1)
%         M = slmetric_pw(X1, X2, 'cityblk');
%
%         % compute the normalize correlations
%         M = slmetric_pw(X1, X2, 'nrmcorr');
%
%         % compute hamming distances
%         M = slmetric_pw(X1, X2, 'hamming', 0.5);
%         M2 = slmetric_pw((X1 > 0.5), (X2 > 0.5), 'hamming');
%         assert(isequal(M, M2));
%     \}
%
%   - Compute the parameterized metrics
%     \{
%         % compute weighted squared distances with user-supplied weights
%         weights = rand(10, 1);
%         M = slmetric_pw(X1, X2, 'wsqdist', weights);
%
%         % compute quadratic distances (x-y)^T * Q (x-y)
%         Q = rand(10, 10);
%         M = slmetric_pw(X1, X2, 'quaddiff', Q);
%
%         % compute Minkowski distance of order 3
%         M = slmetric_pw(X1, X2, 'minkowski', 3);
%     \}
%
% [ History ]
%   - Created by Dahua Lin on Dec 06th, 2005
%   - Modified by Dahua Lin on Apr 21st, 2005
%       - regularize the error reporting
%   - Modified by Dahua Lin on Sep 11st, 2005
%       - completely rewrite the core codes based on new mex computation
%         cores, and the runtime efficiency in both time and space is
%         significantly increased.
%   - Modified by Dahua Lin on Jul 02, 2007
%       - rewrite the core computation based on the bsxfun introduced in
%         MATLAB R2007a
%       - rewrite the core-mex for cityblk, maxdiff, mindiff
%       - introduce new metrics: corrdist, minkowski
%   - Modified by Dahua Lin on Jul 30, 2007
%       - Add the metric types for histograms, which are originally
%         implemented in slhistmetric_pw in sltoolbox v1.
%   - Modified by Dahua Lin on Aug 16, 2007
%       - revise some of the help contents
%

%% parse and verify input arguments
error(nargchk(3, inf, nargin));
assert(ischar(mtype), 'sltoolbox:slmetric_pw:invalidarg', ...
    'The metric type should be a string.');

if strcmp(mtype, 'hamming') || strcmp(mtype, 'hamming_nrm')
    assert((isnumeric(X1) || islogical(X1)) && ndims(X1) == 2 && ...
           (isnumeric(X2) || islogical(X2)) && ndims(X2) == 2, ...
        'sltoolbox:slmetric_pw:invalidarg', 'X1 and X2 should be numeric or logical matrices.');
else
    assert(isnumeric(X1) && ndims(X1) == 2 && isnumeric(X2) && ndims(X2) == 2, ...
        'sltoolbox:slmetric_pw:invalidarg', 'X1 and X2 should be numeric matrices.');
end

assert(isa(X2, class(X1)), ...
    'sltoolbox:slmetric_pw:invalidarg', 'X1 and X2 should be of the same class.');

if isempty(X1) || isempty(X2)
    M = [];
    return;
end

%% compute
switch mtype
    case {'eucdist', 'sqdist'}
        checkdim(X1, X2);
        M = bsxfun(@plus, sum(X1 .* X1, 1)', (-2) * X1' * X2);
        M = bsxfun(@plus, sum(X2 .* X2, 1), M);
        M(M < 0) = 0;
        if strcmp(mtype, 'eucdist')
            M = sqrt(M);
        end

    case 'dotprod'
        checkdim(X1, X2);
        M = X1' * X2;

    case {'nrmcorr', 'corrdist', 'angle'}
        checkdim(X1, X2);
        ns1 = sqrt(sum(X1 .* X1, 1));
        ns2 = sqrt(sum(X2 .* X2, 1));
        ns1(ns1 == 0) = 1;
        ns2(ns2 == 0) = 1;
        M = bsxfun(@times, X1' * X2, 1 ./ ns1');
        M = bsxfun(@times, M, 1 ./ ns2);
        switch mtype
            case 'corrdist'
                M = 1 - M;
            case 'angle'
                M = real(acos(M));
        end

    case 'quadfrm'
        Q = varargin{1};
        M = X1' * Q * X2;

    case 'quaddiff'
        checkdim(X1, X2);
        Q = varargin{1};
        M = X1' * (-(Q + Q')) * X2;
        M = bsxfun(@plus, M, sum(X1 .* (Q * X1), 1)');
        M = bsxfun(@plus, M, sum(X2 .* (Q * X2), 1));

    case 'cityblk'
        checkdim(X1, X2);
        M = pwmetrics_cimp(X1, X2, int32(1));

    case 'maxdiff'
        checkdim(X1, X2);
        M = pwmetrics_cimp(X1, X2, int32(3));

    case 'mindiff'
        checkdim(X1, X2);
        M = pwmetrics_cimp(X1, X2, int32(2));

    case 'minkowski'
        checkdim(X1, X2);
        pord = varargin{1};
        if ~isscalar(pord)
            error('sltoolbox:slmetric_pw:invalidparam', ...
                'the mikowski order should be a scalar');
        end
        pord = cast(pord, class(X1));
        M = pwmetrics_cimp(X1, X2, int32(4), pord);

    case 'wsqdist'
        d = checkdim(X1, X2);
        w = varargin{1};
        if ~isequal(size(w), [d, 1])
            error('sltoolbox:slmetric_pw:invalidparam', ...
                'the weights should be given as a d x 1 vector.');
        end
        wX2 = bsxfun(@times, X2, w);
        M = bsxfun(@plus, (-2) * X1' * wX2, sum(wX2 .* X2, 1));
        clear wX2;
        wX1 = bsxfun(@times, X1, w);
        M = bsxfun(@plus, M, sum(wX1 .* X1, 1)');

    case {'hamming', 'hamming_nrm'}
        checkdim(X1, X2);
        if islogical(X1) && islogical(X2)
            H1 = X1;
            H2 = X2;
        else
            if isempty(varargin)
                t = 0;
            else
                t = varargin{1};
                assert(isnumeric(t) && isscalar(t), ...
                    'sltoolbox:slmetric_pw:invalidparam', 't should be a numeric scalar.');
            end
            H1 = X1 > t;
            H2 = X2 > t;
        end
        M = pwhamming_cimp(H1, H2);
        if strcmp(mtype, 'hamming_nrm')
            M = M / size(H1, 1);
        end

    case 'intersect'
        checkdim(X1, X2);
        M = pwmetrics_cimp(X1, X2, int32(5));

    case 'intersectdis'
        checkdim(X1, X2);
        M = 1 - pwmetrics_cimp(X1, X2, int32(5));

    case 'chisq'
        checkdim(X1, X2);
        M = pwmetrics_cimp(X1, X2, int32(6));

    case 'kldiv'
        checkdim(X1, X2);
        M = pwmetrics_cimp(X1, X2, int32(7));

    case 'jeffrey'
        checkdim(X1, X2);
        M = pwmetrics_cimp(X1, X2, int32(8));

    otherwise
        error('sltoolbox:slmetric_pw:unknowntype', 'Unknown metric type %s', mtype);

end

%% Auxiliary function

function d = checkdim(X1, X2)

d = size(X1, 1);
if d ~= size(X2, 1)
    error('sltoolbox:slmetric_pw:sizmismatch', ...
        'X1 and X2 have different sample dimensions');
end


text:
function [data_train, labels_train, data_devel, labels_devel, raw_devel, PC, means_norm, stds_norm, devel_ids, devel_success] = ...
    Prepare_HOG_AU_data_dynamic(train_users, devel_users, au_train, rest_aus, UNBC_dir, features_dir)

%%
addpath(genpath('../data extraction/'));

% First extracting the labels
[ labels_train, valid_ids_train, filenames ] = extract_UNBC_labels(UNBC_dir, train_users, au_train);

[ labels_other, ~, ~ ] = extract_UNBC_labels(UNBC_dir, train_users, rest_aus);
labels_other = cat(1, labels_other{:});

% Reading in the HOG data (of only relevant frames)
[train_appearance_data, valid_ids_train_hog, vid_ids_train_string] = Read_HOG_files_dynamic(train_users, features_dir);

[train_geom_data] = Read_geom_files_dynamic(train_users,  features_dir);

% Subsample the data to make training quicker
labels_train = cat(1, labels_train{:});
valid_ids_train = logical(cat(1, valid_ids_train{:}));

if(numel(train_users) > 0)
    reduced_inds = false(size(labels_train,1),1);
    reduced_inds(labels_train > 0) = true;

    % make sure the same number of positive and negative samples is taken
    pos_count = sum(labels_train > 0);
    neg_count = sum(labels_train == 0);

    num_other = floor(pos_count / (size(labels_other, 2)));

    inds_all = 1:size(labels_train,1);

    for i=1:size(labels_other, 2)+1

        if(i > size(labels_other, 2))
            % fill the rest with a proportion of neutral
            inds_other = inds_all(sum(labels_other,2)==0 & ~labels_train );
                num_other_i = min(numel(inds_other), pos_count - sum(labels_train(reduced_inds,:)==0));
        else
                % take a proportion of each other AU
            inds_other = inds_all(labels_other(:, i) & ~labels_train );
            num_other_i = min(numel(inds_other), num_other);
        end
        inds_other_to_keep = inds_other(round(linspace(1, numel(inds_other), num_other_i)));
        reduced_inds(inds_other_to_keep) = true;

    end

    % Remove invalid ids based on CLM failing or AU not being labelled
    reduced_inds(~valid_ids_train) = false;
    reduced_inds(~valid_ids_train_hog) = false;

    labels_other = labels_other(reduced_inds, :);
    labels_train = labels_train(reduced_inds,:);
    train_appearance_data = train_appearance_data(reduced_inds,:);
    train_geom_data = train_geom_data(reduced_inds,:);
    vid_ids_train_string = vid_ids_train_string(reduced_inds,:);
end
%% Extract devel data

% First extracting the labels
[ labels_devel, valid_ids_devel, vid_ids_devel ] = extract_UNBC_labels(UNBC_dir, devel_users, au_train);

% Reading in the HOG data (of only relevant frames)
[devel_appearance_data, valid_ids_devel_hog, vid_ids_devel_string] = Read_HOG_files_dynamic(devel_users, features_dir);
devel_success = valid_ids_devel_hog;
devel_ids = vid_ids_devel_string;

[devel_geom_data] = Read_geom_files_dynamic(devel_users, features_dir);

labels_devel = cat(1, labels_devel{:});

% Peforming zone specific masking
if(au_train < 8 || au_train == 43 || au_train == 45) % upper face AUs ignore bottom face
    % normalise the data
    pca_file = '../../pca_generation/generic_face_upper.mat';
    load(pca_file);
elseif(au_train > 9) % lower face AUs ignore upper face and the sides
    % normalise the data
    pca_file = '../../pca_generation/generic_face_lower.mat';
    load(pca_file);
elseif(au_train == 9) % Central face model
    % normalise the data
    pca_file = '../../pca_generation/generic_face_rigid.mat';
    load(pca_file);
end

% Grab all data for validation as want good params for all the data
raw_devel = cat(2, devel_appearance_data, devel_geom_data);

devel_appearance_data = bsxfun(@times, bsxfun(@plus, devel_appearance_data, -means_norm), 1./stds_norm);

data_devel = devel_appearance_data * PC;

data_devel = cat(2, data_devel, devel_geom_data);

if(numel(train_users) > 0)
    train_appearance_data = bsxfun(@times, bsxfun(@plus, train_appearance_data, -means_norm), 1./stds_norm);

    data_train = train_appearance_data * PC;
    data_train = cat(2, data_train, train_geom_data);
else
    data_train = [];
end

geom_size = max(size(train_geom_data, 2), size(devel_geom_data, 2));

PC_n = zeros(size(PC)+geom_size);
PC_n(1:size(PC,1), 1:size(PC,2)) = PC;
PC_n(size(PC,1)+1:end, size(PC,2)+1:end) = eye(geom_size);
PC = PC_n;

means_norm = cat(2, means_norm, zeros(1, geom_size));
stds_norm = cat(2, stds_norm, ones(1, geom_size));

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function dat= proc_subtractMovingAverage(dat, ms, varargin)
%PROC_SUBTRACTMOVINGAVERAGE - Subtract moving average (high-pass) filter
%
%Usage:
% DAT= proc_subtractMovingAverage(DAT, MSEC, <METHOD='causal'>)
%
%Input:
% DAT    - data structure of continuous or epoched data
% MSEC   - length of interval in which the moving average is
%          to be calculated, unit [msec].
% METHOD - 'centered' or 'causal' (default).
%
%Output:
% DAT    - updated data structure

% Author(s): Benjamin Blankertz

misc_checkType(dat, 'STRUCT(x fs)');
misc_checkType(ms, '!DOUBLE[1]');

nSamples = round(ms*dat.fs/1000);
sx= size(dat.x);
xa= procutil_movingAverage(dat.x(:,:), nSamples, varargin{:});
dat.x= dat.x - reshape(xa, sx);

author: bbci


text:
% NTF | betaNTF | Simple beta-NTF implementation (Antoine Liutkus, 2012)
% process_video('NTF', 'betaNTF', 'dataset/demo.avi', 'output/demo_beta-NTF.avi');

% Compute a simple NTF model of 10 components
A = double(T);
r = 10;
[~,~,~,L] = betaNTF(A,r);
S = (A - L);
% For reconstruction
% for i = 1:r, B_hat(:,:,i) = W * diag(Q(i,:)) * H'; end

author: andrewssobral

repo: lrslibrary


text:
classdef DetectionEventHandle < handle
    properties
        detectedEvent;
        symbolHandle;
        textHandle;
    end

    properties (Access = private, Constant)
        visibleStringss = {'off','on'};
    end

    methods (Access = public)
        function obj = DetectionEventHandle(detectedEvent,symbolHandle, textHandle)
            if nargin > 0
                obj.detectedEvent = detectedEvent;
                obj.symbolHandle = symbolHandle;
                obj.textHandle = textHandle;
            end
        end

        function setVisible(obj,visible)
            obj.symbolHandle.Visible = obj.visibleStringss{visible+1};
            obj.textHandle.Visible = obj.visibleStringss{visible+1};
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
%This computes the hand-eye calibration using the method described in
%"Hand-Eye Calibration Using Dual Quaternions" from
%Konstantinos Daniilidis
%
%Input:
%Hmarker2world      a 4x4xNumber_of_Views Matrix of the form
%                   Hmarker2world(:,:,i) = [Ri_3x3 ti_3x1;[ 0 0 0 1]]
%                   with
%                   i = number of the view,
%                   Ri_3x3 the rotation matrix
%                   ti_3x1 the translation vector.
%                   Defining the transformation of the robot hand / marker
%                   to the robot base / external tracking device
%Hgrid2cam          a 4x4xNumber_of_Views Matrix (like above)
%                   Defining the transformation of the grid to the camera
%
%Output:
%Hcam2marker_       The transformation from the camera to the marker /
%                   robot arm
%err                The residuals from the least square processes
%
%Christian Wengert
%Computer Vision Laboratory
%ETH Zurich
%Sternwartstrasse 7
%CH-8092 Zurich
%www.vision.ee.ethz.ch/cwengert
%wengert@vision.ee.ethz.ch
function [Hcam2marker_, err] = hand_eye_dual_quaternion(Hmarker2world, Hgrid2cam)
%Our quaternions are like this (q1 q2 q3,s )
    %Get n
    n = size(Hmarker2world,3);
    %Make movements (a,B) which are interposition transformations
    %(marker2wordl and cam2grid)
    %transform A,B into dual quaternions
    for i=1:n-1
        A = inv(Hmarker2world(:,:,i+1))*Hmarker2world(:,:,i);
        B = Hgrid2cam(:,:,i+1)*inv(Hgrid2cam(:,:,i));
        [q,qprime] = getDualQuaternion(A(1:3,1:3),A(1:3,4));
        Qa(i).q = q;
        Qa(i).qprime = qprime;

        [q,qprime] = getDualQuaternion(B(1:3,1:3),B(1:3,4));
        Qb(i).q = q;
        Qb(i).qprime = qprime;
    end

    %The dual quaternion is (Q.q + epsilon*Q.prime)
    %a = Qa.q, a' = Qa.prime  idem for b
    S = [];
    for i=1:n-1
        S(:,:,i) = [Qa(i).q(1:3)-Qb(i).q(1:3)   crossprod(Qa(i).q(1:3)+Qb(i).q(1:3)) zeros(3,1) zeros(3,3);...
                    Qa(i).qprime(1:3)-Qb(i).qprime(1:3)   crossprod(Qa(i).qprime(1:3)+Qb(i).qprime(1:3)) Qa(i).q(1:3)-Qb(i).q(1:3)   crossprod(Qa(i).q(1:3)+Qb(i).q(1:3))];
    end

    %Construct T
    T = [];
    for i=1:n-1
        T = [T  S(:,:,i)'];
    end

    T = T';
    %SVD
    [U,S,V] = svd(T);

    %Solution, right null vectors of T
    v7 = V(:,7);
    v8 = V(:,8);

    u1 = v7(1:4);
    v1 = v7(5:8);

    u2 = v8(1:4);
    v2 = v8(5:8);
    %Now lambda1*v7+lambda2*v8 = [q;qprime]
    %
    %or other:
    %
    %lambda1^2*u1'*u1+2*lambda1*lambda2*u1'*u2+lambda2^2*u2'*u2 = 1
    %and
    %lambda1^2*u1'*v1 + lambda1*lambda2*(u1'*v2+u2'*v1)+lambda2*u2'*v1 = 0
    %Setting lambda1/lambda2 = s
    %lambda1^2/lambda2^2*u1'*v1 + lambda1*lambda2/lambda2^2*(u1'*v2+u2'*v1)+lambda2^2/lambda2^2*u2'*v1 = 0
    %s^2*u1'*v1 + s*(u1'*v2+u2'*v1)+u2'*v1 = 0
    %s^2*u1'*v1 + s*(u1'*v2+u2'*v1)+u2'*v1 = 0
    a = u1'*v1;
    b = (u1'*v2+u2'*v1);
%     c = u2'*v1;
    c = u2'*v2 ;
    s = roots([a b c]);

    %insert into equation
    val1 = s(1)^2*u1'*u1+2*s(1)*u1'*u2+u2'*u2;
    val2 = s(2)^2*u1'*u1+2*s(2)*u1'*u2+u2'*u2;
    %Take bigger value
    if(val1>val2)
        s = s(1);
        val = val1;
    else
        s = s(2);
        val = val2;
    end
    %Get lambdas
    lambda2 = sqrt(1/val);
    lambda1 = s*lambda2;

    %This algorithm gives quaternion with the form of (s, q1 q2
    %q3)->contrary to the notation we used above (q1 q2 q3,s )
    %Therefore we must rearrange the elements!
    qfinal = lambda1*v7+lambda2*v8;
    q = [qfinal(2:4);qfinal(1)];
    qprime = [qfinal(6:8);qfinal(5)];

    %Extract transformation
    R = q2dcm(q);
    t = 2*qmult(qprime,qconj(q));
    t = t(1:3);

    %Assign output arguments
    Hcam2marker_ = [R -R*t;[0 0 0 1]]^-1;
    err=[];

%Creates a dual quaternion from a rotation matrix and a translation vector
function [q,qprime] = getDualQuaternion(R,t)
    %Conversion from R,t to the screw representation [d,theta,l,m]

    r = rodrigues(R);
    theta = norm(r);
    l = r/norm(theta);
    %Pitch d
    d = l'*t;
    %Make point c
    c = .5*(t-d*l)+cot(theta/2)*cross(l,t);
    %moment vector
    m = cross(c,l);
    %Rotation quaternion
    %(q1 q2 q3,s )
    q = [sin(theta/2)*l; cos(theta/2)];
    %Get dual
    qprime = [.5*(q(4)*t+cross(t,q(1:3)));-.5*q(1:3)'*t];

author: christianwengert


text:
clear all; close all;
I=imread('snowflakes.png');
se=strel('disk', 5);
J=imopen(I, se);
figure;
subplot(121);  imshow(I);
subplot(122);  imshow(J, []);

author: HuangCongQing


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Ex0601: Perform T Test')

SpharmMatDir = '../';
CodeDir = [SpharmMatDir 'code'];
addpath(CodeDir);

DataFolder = 'data';

%% (1) Input data directory
InDataDir = [SpharmMatDir DataFolder '/Ex0601/hip06_reg'];

%% (2) List input file names
inFiles = dir([InDataDir '/*_reg.mat']); inNames={};
for i=1:length(inFiles)
    inNames{end+1} = [InDataDir '/' inFiles(i).name];
end

%% (3) Display input objects (optional)
    %Available values for Space- 'object';'param';'both'
dispConfs.Space = 'object';
    % Available values for Mesh- 'orig';'quad32';'quad64';'quad128';'quad256';'quad512';'icosa1';'icosa2'; ...
    %    'icosa3';'icosa4';'icosa5';'icosa6'
dispConfs.Mesh = 'quad64';
    % Available values for Shape- 'solid';'mesh';'both'
dispConfs.Shade = 'both';
    % Available values for Overlay- 'none';'adc_paramap'
dispConfs.Overlay = 'adc_paramap';
    % Available values for Export- 'screen';'png';'both'
dispConfs.Export = 'png';
dispConfs.Degree = [];
dispConfs.Template = '';

SpharmMatUtilDisplayObjs(dispConfs, inNames, CodeDir);

%% (4) Perform T Test
confs.Atlas = [SpharmMatDir DataFolder '/Ex0601/hip07_atlas/atlas.mat'];
confs.Smoothing_FWHM =5;
    % Available values for EqualVariance- 'Yes';'No'
confs.EqualVariance = 'Yes';
    % Available values for Signal-
    % 'vl_defm_org';'vl_defm_nrm';'vl_defm_pca';'vl_defm_fld'
confs.Signal = 'vl_defm_nrm';
    % Available values for SampleMesh-
    % 'icosa1';'icosa2';'icosa3';'icosa4';'icosa5';'icosa6'
confs.SampleMesh = 'icosa3';
confs.OutputNamePrefix = 't_map';
confs.GroupIDs = 'Ctrl,PT';
confs.OutDirectory = [SpharmMatDir DataFolder '/Ex0601/hip08_stat'];

if ~exist(confs.OutDirectory,'dir')
    mkdir(confs.OutDirectory);
end

inFiles = dir([InDataDir '/a*_reg.mat']); inNames1={};
for i=1:length(inFiles)
    inNames1{end+1} = [InDataDir '/' inFiles(i).name];
end

inFiles = dir([InDataDir '/b*_reg.mat']); inNames2={};
for i=1:length(inFiles)
    inNames2{end+1} = [InDataDir '/' inFiles(i).name];
end

SpharmMatStatAnalysis(confs, inNames1, inNames2, 't_map', CodeDir);

%% (5) Display output objects (optional)
dispConfs.Threshold_p_value = 0.05;
    % Available values for Overlay- 'p-value';'t-map'
dispConfs.Overlay = 'p-value';
    % Available values for Colormap-
    % 'jet';'hot';'summer';'cool';'autumn';'winter';'spring'
dispConfs.Colormap = 'jet';

inFile = dir([confs.OutDirectory '/t_map_*.mat']); inName={};
for i=1:length(inFile)
    inName{end+1} = [confs.OutDirectory '/' inFile(i).name];
end

h = figure('Name', 'Render Statistical Results','NumberTitle', 'off');
cameratoolbar(h, 'Show');
SpharmMatDisplayStat(dispConfs, inName, 'res_t_map', h, CodeDir);

dispConfs.Overlay = 't-map';

h = figure('Name', 'Render Statistical Results','NumberTitle', 'off');
cameratoolbar(h, 'Show');
SpharmMatDisplayStat(dispConfs, inName, 'res_t_map', h, CodeDir);

rmpath(CodeDir);

clear;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function vbrfa2011_artificial_experiment(seed)

if nargin < 1
  seed = 10;
end
rand('state', seed);
randn('state', seed);

% Observations:
%
% - When increasing M, Laplace -> independent t
%
% - Variance of corruption makes no big difference between Laplace and
% independent t

% From Bishop, standard deviations 5,4,3,2,1,1,1,1,... :

M = 10;
N = 100;
D = 4;

% Covariance matrix singular values:

R = orth(randn(M));
eig = ([D:-1:1,zeros(1,M-D)]).^2;
Cov = R*diag(eig)*R';
% $$$ eig = (1+[D:-1:1,zeros(1,M-D)]).^2;
% $$$ Cov = R*diag(eig-1)*R';

mu = randn(M,1);

% Noiseless data
Y = mvnrnd(mu, Cov, N)';

% Noisy data
Yn = Y + 1*randn(M,N);

% Corrupted data
I = rand(M,N) < 0.02;
Yno = Yn + I.*unifrnd(-30,30,M,N);

%
% Construct the model
%

Dh = M-1;

% $$$ tsplot(Yno, 'k')
% $$$ addtsplot(Y, 'b');
% $$$ return

for model=1:4

  % It seems that (at least in some cases) it would be better to combine the
  % isotropic noise with X instead W for better and more stable results..?

  noise_iso = noise_module_isotropic(M, 1, ...
                                     'init', struct('a_tau',1, ...
                                                    'b_tau',1));

  % ARD for W
  W_module = factor_module_ard(Dh+1, M, ...
                               'update_alpha', 1, ...
                               'noise_module', noise_iso);

  % IID for X
  mu = [1; zeros(Dh,1)];
  Cov = diag([1e-6; ones(Dh,1)]);
  X_module = factor_module_iid(Dh+1, N, ...
                               'prior', struct('mu', mu, ...
                                               'CovX', Cov), ...
                               'noise_module', []);

  switch model

   case 1 % GAUSSIAN
    noise = noise_module_fixed(M, N, 1);
    id = 'gaussian';

   case 2 % MULTIVARIATE T
    noise = noise_module_multivariate_t(M, N, ...
                                        'update_nu', [1], ...
                                        'init', struct('nu', 0.1));
    id = 'multi-t';

   case 3 % INDEPENDENT T
    noise = noise_module_independent_t(M, N, ...
                                       'update_nu', [1], ...
                                       'nu', 'pooled', ...
                                       'init', struct('nu', 0.1));
    id = 'ind-t';

   case 4 % LAPLACE
    noise = noise_module_laplace(M, N);
    id = 'laplace';

   otherwise
    error('Unknown model requested')

  end

  %
  % Run the method
  %

  % General function call for running the method
  Q = vbfa(Dh+1, ...
           Yno, ...
           W_module, ...
           X_module, ...
           noise, ...
           'maxiter', 50, ...
           'rotate', [1:10, 10:5:100], ...
           'rotation_maxiter', 20, ...
           'update_x', 1, ...
           'update_w', 1, ...
           'update_noise', 1);

  Yh = Q.W'*Q.X;

  rmse1_train(model) = rmse(Yh, Yno)
  rmse2_test(model) = rmse(Yh(~I), Y(~I))
  rmse3_outliers(model) = rmse(Yh(I), Y(I))

% $$$   tsplot(Yno(1:10,:), 'k')
% $$$   addtsplot(Y(1:10,:), 'b');
% $$$   addtsplot(Yh(1:10,:), 'r');
end

rmse1_train
rmse2_test
rmse3_outliers

%Q.W*diag(Q.rhow)*Q.W' + sum(Q.CovW,3)


text:
function [y, extra, invhess] = glmevfwd(net, x, t, x_test, invhess)
%GLMEVFWD Forward propagation with evidence for GLM
%
%	Description
%	Y = GLMEVFWD(NET, X, T, X_TEST) takes a network data structure  NET
%	together with the input X and target T training data and input test
%	data X_TEST. It returns the normal forward propagation through the
%	network Y together with a matrix EXTRA which consists of error bars
%	(variance) for a regression problem or moderated outputs for a
%	classification problem.
%
%	The optional argument (and return value)  INVHESS is the inverse of
%	the network Hessian computed on the training data inputs and targets.
%	Passing it in avoids recomputing it, which can be a significant
%	saving for large training sets.
%
%	See also
%	FEVBAYES
%

%	Copyright (c) Ian T Nabney (1996-2001)

[y, a] = glmfwd(net, x_test);
if nargin == 4
  [extra, invhess] = fevbayes(net, y, a, x, t, x_test);
else
  [extra, invhess] = fevbayes(net, y, a, x, t, x_test, invhess);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [] = dualcolorplot(x,y,F,OPT)
%DUALCOLORPLOT linear plot with dividing line.
% DUALCOLORPLOT(X,Y,F) plots vector Y versus vector X, and places a
% horizontal line on the plot at y-value F.  The line divides the plot into
% red (top) and blue (bottom).  The user is able to click on the dividing
% line then drag it up and down to change the proportion of red to blue.
%
% DUALCOLORPLOT(X,Y,F,'interp') forces linear interpolation between top and
% bottom colors so that they meet at the horizontal line.  This is slightly
% slower than the default for very large vectors.
%
% A context menu is provided for each of the lines in order to alter their
% appearance.  Simply right-click on any of the lines to see the
% options.  A context menu is also provided for the axes to turn on or off
% automatic scaling.  This can be useful if in zoom mode.
%
% NOTE:
% Calling DUALCOLORPLOT will set the following properties of the current
%         figure:  windowbuttondownfcn
%                  windowbuttonupfcn,
%                  windowbuttonmotionfcn
%
% and the uicontextmenu for the current axes.
%
% THIS FUNCTION REQUIRES LINECMENU BE ON THE PATH.
%
%
% EXAMPLES:
%           x = -2*pi:.01:2*pi;
%           f = @(x) sin(x).*exp(x/5);
%           dualcolorplot(x,f(x),0)  % Right click on lines and axes!
%
%
%           % Show the difference interpolation makes by exaggeration...
%           x = 0:.25:8*pi;  % A thinly populated plot
%           y = sin(x);
%           subplot(1,2,1)
%           dualcolorplot(x,y,0)
%           title('No Interpolation.')
%           subplot(1,2,2)
%           dualcolorplot(x,y,0,'interp')
%           title('With Interpolation')
%
%
% See also plot
%
% Author:  Matt Fig
% Data:  3/30/2011

NARG = nargin;

if NARG==3
    flg = 0;
elseif NARG == 4
   if ~ischar(OPT) || ~strncmpi(OPT,'interp',6)
      error('Unknown Fourth option.  See help.')
   end
   flg = 1;  % User wants to interpolate.
else
    error('dualcolorplot requires 3 or 4 input arguments.  See help.')
end

fh = get(0,'currentfigure');

if isempty(fh)
    fh = figure('windowbuttondownfcn',{@fh_wndbtnfcn},...
                'windowbuttonupfcn',{@fh_wndbtnfcn},...
                'windowbuttonmotionfcn',{@fh_wndmtnfcn});
    AX = axes;
else
    set(fh,'windowbuttondownfcn',{@fh_wndbtnfcn},...
           'windowbuttonupfcn',{@fh_wndbtnfcn},...
           'windowbuttonmotionfcn',{@fh_wndmtnfcn});
    AX = get(fh,'currentaxes');

    if isempty(AX)
        AX = axes;
    end
end

L = [];  % This hold the handles to the lines.
V = 1:length(x);  % For indexing... avoid FIND.
mxx = max(x(:));  mnx = min(x(:));  % For the axis limits.
mxy = max(y(:));  mny = min(y(:));
D = (mxy-mny)/10;  % Resize buffer.
fbd = false;  % figure button up/down.
axzm =  true;  % On and off resize axes...  for zoom state.
LUICM = linecmenu;  % Call the LINECMENU function.
AUICM = uicontextmenu('userdata',axzm);
AUIM = uimenu(AUICM,'label','Set Resize OFF');
set(AUIM,'callback',{@ax_cntxtcb});  % Callback for uimenu.
set_colors();  % This is the function which does the work.

    function [] = set_colors(varargin)
    % Sets the appropriate colors, and divides x and y.
        idxlt = y<=F;  % Find the lower portion of the data.
        idxgt = V(~idxlt);
        idxlt = V(idxlt);
        y1 = y;  x1 = x;  % The lower portion of the plot.
        y2 = y;  x2 = x;  % The upper portion of the plot.
        y1(idxgt) = nan;  % NAN won't get plotted.
        y2(idxlt) = nan;

        if flg % User wants interpolated values.
            NN = isnan(y1);
            A = strfind(NN,[1 0]);  % Where interpolations will go.
            B = strfind(NN,[0 1]);
            y1(A) = F;
            y1(B+1) = F;  % At interpolation point, y is the horiz line.
            y2(B) = F;
            y2(A+1) = F;
            % Interpolations for the x variables.
            idx1 = (x(A+1)-x(A))./(y(A+1)-y(A)).*(F-y(A))+x(A);
            idx2 = (x(B+1)-x(B))./(y(B+1)-y(B)).*(F-y(B))+x(B);
            x1(A) = idx1;
            x1(B+1) = idx2;
            x2(B) = idx2;
            x2(A+1) = idx1;
        end

        if isempty(L)  % First time through
            L(1) = line(x2,y2,'color','r');
            L(2) = line(x1,y1,'color','b');
            L(3) = line([mnx mxx],[F F],'color','g',...
                        'tag','HorizLine');
                    % Store data in Horz line userdata.
            set(L(3),'userdata',{L,flg,x,y,D,mnx,mxx,mny,mxy})
            set(AX,'uicontextmenu',AUICM);
            axis([mnx mxx min(mny,F)-D max(mxy,F)+D]) % Adjust axis
            set(L,'uicontextmenu',LUICM)
        else  % Called from callback.
            UD = get(varargin{1},'userdata');
            set(UD{1}(1),'xdata',x2,'ydata',y2);
            set(UD{1}(2),'xdata',x1,'ydata',y1);
            UD{3} = x;
            UD{4} = y;
            set(varargin{1},'userdata',UD)
        end
    end

    function [] = fh_wndbtnfcn(varargin)
    % Windowbutton-up/down-function for figure.
        fbd = ~fbd; % Button up or down?
    end

    function [] = fh_wndmtnfcn(varargin)
    % Windowbuttonmotionfcn for figure.
        CL = get(fh,'currentobject');
        CO = findobj(fh,'tag','HorizLine');  % The horizontal line(s).
        if ~isempty(CL) && any(CL==CO)  % Is current one of the horiz?
            UD = get(CL,'userdata') ;
            flg = UD{2};  % Retrieve the flag for interpolation.

            if fbd  % Button is held down.
                AXD = get(gca,{'currentpo','uicontextm'});% AXES DATA
                set(CL,'ydata',AXD{1}(:,2));

                if get(AXD{2},'userdata')
                    axis([UD{6} UD{7} min(UD{8},AXD{1}(1,2))-...
                         UD{5} max(UD{9},AXD{1}(1,2))+UD{5}])
                end

                F = AXD{1}(1,2);
                x = UD{3};
                y = UD{4};
                set_colors(CL);  % Call the drawing function.
            end

        end
    end

    function [] = ax_cntxtcb(varargin)
    % Callback for the contextmenus to the axes.
        P = get(varargin{1},'parent');
        axzm = ~get(P,'userdata');  % The on/off state is stored here.

        if axzm  % Change the label for the next call.
            set(AUIM,'label','Set Resize OFF')
        else
            set(AUIM,'label','Set Resize ON')
        end

        set(P,'userdata',axzm)  % Set the new value.
    end
end


text:
function s = resampstr(p,m,n)
%RESAMPSTR Stratified resampling
%
%  Description
%    S = RESAMPSTR(P) returns a new set of indices according to
%      the probabilities P. P is array of probabilities, which are
%      not necessarily normalized, though they must be
%      non-negative, and not all zero. The size of S is the size of P.
%
%    S = RESAMPSTR(P,M,N) returns an M by N matrix.
%
%    Default is to use no-sort resampling. For sorted resampling use
%      [PS,PI]=SORT(P);
%      S=PI(RESAMPSTR(PS));
%    Sorted re-sampling is slower but has slightly smaller
%    variance. Stratified resampling is unbiased, almost as fast as
%    deterministic resampling (RESAMPDET), and has only a slightly
%    larger variance.
%
%    In stratified resampling indices are sampled using random
%    numbers u_j~U[(j-1)/n,j/n], where n is length of P. Compare
%    this to simple random resampling where u_j~U[0,1].
%
%  Reference
%    Kitagawa, G., Monte Carlo Filter and Smoother for Non-Gaussian
%    Nonlinear State Space Models, Journal of Computational and
%    Graphical Statistics, 5(1):1-25, 1996.
%
%  See also
%    RESAMPSIM, RESAMPRES, RESAMPDET
%
% Copyright (c) 2003-2004,2010 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin<2
  [m,n]=size(p);
elseif nargin==2
  n=m;
end
mn=m.*n;
pn=p./sum(p(:)).*mn;
s=zeros(m,n);
r=rand(1,mn);
k=0;
c=0;
for i=1:numel(pn)
  c=c+pn(i);
  if c>=1
    a=floor(c);
    c=c-a;
    s(k+[1:a])=i;
    k=k+a;
  end
  if k<mn && c>=r(k+1)
    c=c-1;
    k=k+1;
    s(k)=i;
  end
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function cS = cordierite

cs_Crd = crystalSymmetry('mmm',[0.568 1 0.549]);
N = Miller({1,0,0},{0,1,0},{0,0,1},{1,1,0},{1,1,1},{1,3,0},{1,1,2},{0,1,1},cs_Crd);
dist = [2.1, 1.05, 3.5, 2.25, 5.3, 3.65, 8.7, 4.45];
cS = crystalShape(N./dist);

author: mtex-toolbox

repo: mtex

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


text:
function timer_cputime_test03 ( )

%*****************************************************************************80
%
%% TEST03 times the unvectorized EXP routine.
%
%  Discussion:
%
%    For the MATLAB implementation of this test, the unvectorized
%    test limits had to be DRASTICALLY reduced.  But then I noticed
%    that some kind of MATLAB cleverness causes all the loops
%    after the first cycle to be computed in almost no time.
%    This test will have to be rethought.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 December 2005
%
%  Author:
%
%    John Burkardt
%
  n_log_min = 12;
  n_log_max = 16;
  n_min = 2^n_log_min;
  n_max = 2^n_log_max;
  n_rep = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST03\n' );
  fprintf ( 1, '  Time the unvectorized loops:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    for i = 1 : n\n' );
  fprintf ( 1, '      y(i) =        x(i)\n' );
  fprintf ( 1, '      y(i) = PI *   x(i)\n' );
  fprintf ( 1, '      y(i) = sqrt ( x(i) )\n' );
  fprintf ( 1, '      y(i) = exp  ( x(i) )\n' );
  fprintf ( 1, '    end\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data vectors will be of minimum size %d\n', n_min );
  fprintf ( 1, '  Data vectors will be of maximum size %d\n', n_max );
  fprintf ( 1, '  Number of repetitions of the operation: %d\n', n_rep );

  for func = 1 : 4

    for i_rep = 1 : n_rep

      for n_log = n_log_min : n_log_max

        n = 2^n_log;

        x = rand(n,1);

        time1 = cputime;

        if ( func == 1 )
          for i = 1: n
            y(i) = x(i);
          end
        elseif ( func == 2 )
          for i = 1: n
            y(i) = pi * x(i);
          end
        elseif ( func == 3 )
          for i = 1: n
            y(i) = sqrt ( x(i) );
          end
        elseif ( func == 4 )
          for i = 1: n
            y(i) = exp ( x(i) );
          end
        end

        time2 = cputime;

        delta(n_log,i_rep) = time2 - time1;

      end

    end

    fprintf ( 1, '\n' );
    fprintf ( 1, 'Test03 Results:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'Vector Size  Rep #1        Rep #2        ' );
    fprintf ( 1, 'Rep #3        Rep #4        Rep #5\n' );
    fprintf ( 1, '\n' );

    for n_log = n_log_min : n_log_max
      n = 2^n_log;
      fprintf ( 1, '%d', n );
      for j = 1 : n_rep
        fprintf ( 1, '%14f', delta(n_log,j) );
      end
      fprintf ( 1, '\n' );
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function F = in_fread_manscan(sFile, sfid, iEpoch, SamplesBounds)
% IN_FREAD_MANSCAN:  Read a block of recordings from a MANSCAN file
%
% USAGE:  F = in_fread_manscan(sFile, sfid, iEpoch, SamplesBounds)  : Read all channels
%         F = in_fread_manscan(sFile, sfid)                         : Read all channels, all the times

% @=============================================================================
% 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, 2012-2014

if (nargin < 3) || isempty(iEpoch)
    iEpoch = 1;
end
if (nargin < 4) || isempty(SamplesBounds)
    if ~isempty(sFile.epochs)
        SamplesBounds = round(sFile.epochs(iEpoch).times .* sFile.prop.sfreq);
    else
        SamplesBounds = round(sFile.prop.times .* sFile.prop.sfreq);
    end
end

% Read data block
nChannels  = length(sFile.header.epoch(1).ChannelOrder);
nTimes = SamplesBounds(2) - SamplesBounds(1) + 1;
% Epoch offset
epochOffset = sFile.header.epoch(iEpoch).StartData1;
% Time offset
timeOffset = 2 * SamplesBounds(1) * nChannels;
% Total offset
totalOffset = epochOffset + timeOffset;

% Set position in file
fseek(sfid, totalOffset, 'bof');
% Read value
F = fread(sfid, [nChannels,nTimes], 'int16');

% Apply gains
F = bst_bsxfun(@rdivide, double(F), double(sFile.header.Gains));

author: brainstorm-tools


text:
function [res] = runsvm(Ks,lk)
% Copyright 2012 Nino Shervashidze, Karsten Borgwardt
% runsvm(Ks,lk)
% K = 1 x h cell array of kernelmatrices (n*n)
% lk = vector of class labels (n*1)
% cv = number of folds in cross-validation

% independent scheme
% best c

addpath('~/code/libsvm');
n=length(lk) % size of the dataset
% randomly permute labels: r will be also used for permuting the kernel matrices
r = randperm(n);
lk = lk(r);

% specify range of c-values
cvalues = (10 .^ [-7:2:7]) / size(lk,1);

cv = 10;
p80 = ceil(n * (1-2/cv));
p90 = ceil(n * (1-1/cv));
fs = n - p90; % fold size

% output variables
res.optkernel=zeros(cv,1);
res.optc=zeros(cv,1);
res.accuracy=zeros(cv,1);

% cross-validation loop
opth=zeros(1,cv);
for k = 1:cv
  imresult=[];

  height = length(Ks);
  for h=1:height
    K = Ks{h}(r,r);
    K_current = K([k*fs+1:size(K,2),1:(k-1)*fs,(k-1)*fs+1:k*fs],[k*fs+1:size(K,2),1:(k-1)*fs,(k-1)*fs+1:k*fs]);
    lk_current = lk([k*fs+1:size(K,2),1:(k-1)*fs,(k-1)*fs+1:k*fs]);
    K_current = makepos(K_current);
    K1 = [(1:size(K_current,1))', normalizekm(K_current)];

    for i = 1:size(cvalues,2)
      % train on 80%, predict on 10% (from 81% to 90%)
      size(lk_current(1:p80));
      size(K1(1:p80,1:p80+1));
      model = svmtrain(lk_current(1:p80,1), K1(1:p80,1:p80+1), strcat(['-t 4  -c ' num2str(cvalues(i))]));
      [predict_label, accuracy, dec_values] = svmpredict(lk_current(p80+1:p90,1),K1(p80+1:p90,1:p80+1), model);
      imresult(h,i)= accuracy(1);
    end
  end

  % determine optimal h and c
  [junk,position]= max(imresult(:));
  [optimalh, indoptimalc]=ind2sub(size(imresult),position);

  opth(k)=optimalh;
  res.optc(k)= cvalues(indoptimalc);
  res.optkernel(k)=optimalh;
  % train on 90% with optimal c, predict on 10% (from 91% to 100%)
  K = Ks{optimalh}(r,r);
  K_current = K([k*fs+1:size(K,2),1:(k-1)*fs,(k-1)*fs+1:k*fs],[k*fs+1:size(K,2),1:(k-1)*fs,(k-1)*fs+1:k*fs]);
  lk_current = lk([k*fs+1:size(K,2),1:(k-1)*fs,(k-1)*fs+1:k*fs]);
  K_current = makepos(K_current);
  K1 = [(1:size(K_current,1))', normalizekm(K_current)];

  model = svmtrain(lk_current(1:p90,1), K1(1:p90,1:p90+1),strcat(['-t 4  -c ' num2str(cvalues(indoptimalc))]) );
  [predict_label, accuracy, dec_values] = svmpredict(lk_current(p90+1:size(K,1),1), K1(p90+1:size(K,1),1:p90+1), model);
  res.accuracy(k)=accuracy(1)
end
res.mean_acc =  mean(res.accuracy)
res.std_acc = std(res.accuracy)
end

function result = makepos(K)
pd = 0;
addc = 10e-7;
while (pd ==  0)

addc = addc * 10
try
if (isinf(addc) == 1)
pd = 1;
else
chol(normalizekm(K + eye(size(K,1),size(K,1)) * addc));
pd = 1;
end
catch

end

end
if (isinf(addc)==0)
result = K + eye(size(K,1),size(K,1)) * addc;
else
result = eye(size(K,1));
end
end

author: muhanzhang

repo: DGCNN

save_path: github-repos/MATLAB/muhanzhang-DGCNN


text:
function c = phaseunlockreal(c,a,M,varargin)
%PHASEUNLOCKREAL  Undo phase lock of Gabor coefficients
%   Usage:  c=phaseunlockreal(c,a,M);
%
%   `phaseunlockreal(c,a,M)` removes phase locking from the Gabor coefficients *c*.
%   The coefficient must have been obtained from a |dgtreal| with parameter *a*
%   and using the 'timeinv' flag.
%
%   Phase locking the coefficients modifies them so as if they were obtained
%   from a frequency-invariant Gabor system. A filter bank produces phase locked
%   coefficients.
%
%   See also: dgt, phaselockreal
%
%   References: puc95

%   AUTHOR:    Peter Balazs, Peter L. Søndergaard, Zdenek Prusa
%   TESTING:   OK
%   REFERENCE: OK

if nargin<3
  error('%s: Too few input parameters.',upper(mfilename));
end;

if  ~isscalar(a) || ~isnumeric(a) || rem(a,1)~=0
  error('a must be integer');
end;

if  ~isscalar(M) || ~isnumeric(M) || rem(M,1)~=0
  error('M must be integer');
end;

definput.flags.accuracy={'normal', 'precise'};
flags=ltfatarghelper({},definput,varargin);

M2=size(c,1);
M2user = floor(M/2) + 1;

if M2~=M2user
    error('%s: Size of s does not comply with M.',upper(mfilename));
end

N=size(c,2);

if flags.do_normal
    TimeInd = (0:(N-1))*a;
    FreqInd = (0:(M2-1));

    phase = FreqInd'*TimeInd;
    phase = mod(phase,M);
    phase = exp(-2*1i*pi*phase/M);

    % Handle multisignals
    c=bsxfun(@times,c,phase);
elseif flags.do_precise

    frames = ifftreal(c,M);
    for n=0:N-1
        frames(:,n+1,:) = circshift(frames(:,n+1,:),n*a);
    end
    c = fftreal(frames);

end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function w = moments_normal_01 ( m )

%*****************************************************************************80
%
%% MOMENTS_NORMAL_01 returns moments of the standard Normal distribution.
%
%  Discussion:
%
%    pdf(x) = exp ( -x^2/2 ) / sqrt ( pi * 2 )
%    mu(k) = integral ( -oo < x < +oo ) x^k pdf(x) dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of moments desired.
%
%    Output, real W(M), the weighted integrals of X^0 through X^(M-1).
%
  w = zeros ( m, 1 );

  w(1) = 1.0;

  for k = 3 : 2: m
    w(k) = r8_factorial2 ( k - 2 );
  end

  for k = 2 : 2 : m
    w(k) = 0.0;
  end

  return
end


text:
close all
clear

%------------------------------------------------------------------------------------------
nb = 1000;%Total number of initial states
nmpcd=zeros(1,nb);%Store computing time for one time-step for NMPC-DCBF
impcd=zeros(1,nb);%Store computing time for one time-step for iMPC-DCBF
nmpcinf=zeros(1,1);%Store infeasible rate for one time-step for NMPC-DCBF
impcinf=zeros(1,1);%Store infeasible rate for one time-step for iMPC-DCBF
i=24;%i is number of horozon
ts = 1;%Total time steps
load(gamma1_0p4gamma2_0p4\InitialStateData');%Load the generated 1000 initial states
[a, b, c, d]=compare(i,ts,nb,InitialMat);%Involves NMPC-DCBF and iMPC-DCBF
nmpcd(1,:)=a;
impcd(1,:)=b;
nmpcinf(1,1)=c;
impcinf(1,1)=d;
nmpcplot1=nmpcd(1,:);
nmpcplot1(nmpcplot1==-1)=[];%Eliminate infeasible one time-step trajectories
impcplot1=impcd(1,:);
impcplot1(impcplot1==-1)=[];%Eliminate infeasible one time-step trajectories

save('timecom24','nmpcplot1','impcplot1');
save('feasibility24','nmpcinf','impcinf');

distnmpc = fitdist(nmpcplot1','Normal');
distimpc = fitdist(impcplot1','Normal');
mu1=distnmpc.mu;%Get mean of sample of computing time for NMPC-DCBF
mu2=distimpc.mu;%Get mean of sample of computing time for iMPC-DCBF
sigma1=distnmpc.sigma;%Get variance of sample of computing time for NMPC-DCBF
sigma2=distimpc.sigma;%Get variance of sample of computing time for iMPC-DCBF

%Initialize atmosphere parameters
function [tnmpc, timpc, ratiotnmpc, ratiotimpc]=compare(N11,ttsim,samplen,InitialMat)
tnmpc=[];
timpc=[];
for i=1:samplen
N1=N11;
tsim=ttsim;
xini1=InitialMat(1,i);
yini1=InitialMat(2,i);
thetaini1=InitialMat(3,i);
vini1=InitialMat(4,i);
x01=[xini1;yini1;thetaini1;vini1];
x02=[xini1 yini1 thetaini1 vini1];
t1 = nmpcdcbf(x01, N1, tsim);
t2 = impcdcbf(x02, N1, tsim);
tindex=[N11 i];
disp(tindex);
tnmpc=[tnmpc t1];%Computing time for NMPC-DCBF
timpc=[timpc t2];%Computing time for iMPC-DCBF
end
nnmpc1 = length(tnmpc);
nimpc1 = length(timpc);
tnmpcs = tnmpc;
timpcs = timpc;
tnmpcs(tnmpcs==-1)=[];
timpcs(timpcs==-1)=[];
nnmpc2 = length(tnmpcs);
nimpc2 = length(timpcs);
ratiotnmpc = (nnmpc1-nnmpc2)/nnmpc1;%Infeasible rate for NMPC-DCBF
ratiotimpc = (nimpc1-nimpc2)/nimpc1;%Infeasible rate for iMPC-DCBF
end

function t1 = nmpcdcbf(x00, N1, ttsim)
%% General Flags
run_nmpc_dcbf_one = true;
run_nmpc_dcbf_multiple = true;

%% Setup and Parameters
x0 = x00;
dt = 0.1;
time_total = ttsim *dt;
P = zeros(4,4);%Weight matrix P
P(1,1) = 10; P(2,2) = 10; P(3,3) = 10;P(4,4) = 10;
Q = zeros(4,4);%Weight matrix Q
Q(1,1) = 10; Q(2,2) = 10; Q(3,3) = 10;Q(4,4) = 10;
R = zeros(4,4);%Weight matrix R
R(1,1) = 1; R(2,2) = 1;R(3,3) = 1000;R(4,4) = 1000;
%Variables range as below
xmin = [-10; -10; -10;-10];
xmax = [10; 10; 10;10];
umin = [-7; -5;-inf;-inf];
umax = [7; 5;inf;inf];

%% Discrete-time unicycle model
system.dt = dt;
system.A = [1 0 0 0;
    0 1 0 0;
    0 0 1 0;
    0 0 0 1];
system.B = [x0(4,1)*cos(x0(3,1))*dt;
    x0(4,1)*sin(x0(3,1))*dt;
    0;
    0];
system.C =[0 0 0 0;
    0 0 0 0;
    1*dt 0 0 0;
    0 1*dt 0 0];
system.xl = xmin;
system.xu = xmax;
system.ul = umin;
system.uu = umax;

%% NMPC-DCBF parameters
params_nmpc_dcbf.Q = Q;
params_nmpc_dcbf.R = R;
params_nmpc_dcbf.P = P;
params_nmpc_dcbf.gamma1 = 0.4;
params_nmpc_dcbf.gamma2 = 0.4;

%% Obstacle
obs.pos1 = [0; 0];
obs.r1 = 1;

%% Simulate NMPC-DCBF with other N
params_nmpc_dcbf.N = N1;
if run_nmpc_dcbf_one
    fprintf('Run NMPC-DCBF-24\n');
    controller_nmpc_dcbf_multiple24 = NMPCDCBF2(x0, system, params_nmpc_dcbf);%Highest order mcbf=2
    controller_nmpc_dcbf_multiple24.obs = obs;
    controller_nmpc_dcbf_multiple24.sim(time_total);
end
t1=controller_nmpc_dcbf_multiple24.tt;
end

function t2=impcdcbf(x00, N1, ttsim)
t2=0;
xo = 0;
yo = 0;
r = 1;
N = N1;%Number of horizon
K1 = 1000;%Maximum iteration times, jmax
K = 1000;%Maximum iteration times, jmax
dt = 0.1;
gamma1 = 4;
gamma2 = 4;
nsim = ttsim;

tic;

% Constraints

umin = [-7; -5; -inf; -inf;];
umax = [7; 5; Inf; Inf];
xmin = [-10; -10; -10; -10];
xmax = [ 10;  10;  10;  10];

% Objective function
Q = diag([10 10 10 10]);
QN = Q;
R = 1 * eye(4);
R(3,3) = 1000;
R(4,4) = 1000;
% Initial and reference states
x0 = x00;
xr = [3; 0.01; 0; 0];
ur = [0; 0; 1; 1];

% Dynamic system initialization
BB = [0 0 0 0;0 0 0 0;dt 0 0 0;0 dt 0 0];

% Convex MPC frame
[nx, nu] = size(BB);
% Initialize states set (x00,x01,,,x0N）
x_0 = [];
x0 = x0';
u_0 = zeros(nu, N);
u0 = zeros(nu, 1);% we may need to change this for a better warm-up start
impc = zeros(nx, nsim + 1);
impc(:, 1) = x0;
x0new = x0;
for i=1 : (N+1)
    x_0 = [x_0 x0new];
    x0new = [dt*x0new(4)*cos(x0new(3))+x0new(1);dt*x0new(4)*sin(x0new(3))+x0new(2);dt*u0(1)+x0new(3);dt*u0(2)+x0new(4)];
end

abc = tangentline(x_0, r);
Acbf = cbfmat(abc, nx, nu, dt, gamma1, x0);%First order CBF constraints
lcbf = -abc(3, :);
lcbf(1) = [];
lcbf = lcbf';
ucbf = inf * ones(N, 1);
A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0);%Second order CBF constraints
lcbf2 = getlcbf2(abc, dt, gamma1, gamma2);
ucbf2 = inf * ones(N-1, 1);

% Cast MPC problem to a QP: x = (x(0),x(1),...,x(N),u(0),...,u(N-1))
% - quadratic objective
P = blkdiag( kron(speye(N), Q), QN, kron(speye(N), R) );
% - linear objective
q = [repmat(-Q*xr, N, 1); -QN*xr; repmat(-R*ur, N, 1)];
% - linear dynamics
Ax = getax(x_0, dt, N, nx);
Bu = kron([sparse(1, N); speye(N)], BB);
Cx = getcx(x_0, u_0, dt, N, nx);
Aeq = [Ax, Bu];
leq = [-x0; -Cx];
ueq = leq;
% - input and state constraints
Aineq = speye((N+1)*nx + N*nu);
lineq = [repmat(xmin, N+1, 1); repmat(umin, N, 1)];
uineq = [repmat(xmax, N+1, 1); repmat(umax, N, 1)];

% - OSQP constraints
A = [Aeq; Aineq; Acbf; A2cbf];
l = [leq; lineq; lcbf; lcbf2];
u = [ueq; uineq; ucbf; ucbf2];

% Create an OSQP object
prob = osqp;
% Setup workspace
prob.setup(P, q, A, l, u, 'warm_start', true);
% Solve
res = prob.solve();

% Check solver status
if ~strcmp(res.info.status, 'solved')
%    error('OSQP did not solve the problem!')
        t2=-1;
        return
end
ctrlx = res.x(1:(N+1)*nx);
ctrlu = res.x((N+1)*nx+1:(N+1)*nx+N*nu);

for j = 1 : (K1-1)%Iterations for the first time-step
    storagex = ctrlx;%store the state variables and inputs in order to compare them with the next iteration
    storageu = ctrlu;
    x_0 = trans(x0, ctrlx, nx, N+1);
    u_0 = transu(ctrlu, nu, N);
    abc = tangentline(x_0, r);
    Acbf = cbfmat(abc, nx, nu, dt, gamma1, x0);
    lcbf = -abc(3, :);
    lcbf(1) = [];
    lcbf = lcbf';
    ucbf = inf * ones(N, 1);
    A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0);
    lcbf2 = getlcbf2(abc, dt, gamma1, gamma2);
    ucbf2 = inf * ones(N-1, 1);
    % Cast MPC problem to a QP: x = (x(0),x(1),...,x(N),u(0),...,u(N-1))
    % - quadratic objective
    P = blkdiag( kron(speye(N), Q), QN, kron(speye(N), R) );
    % - linear objective
    q = [repmat(-Q*xr, N, 1); -QN*xr; repmat(-R*ur, N, 1)];
    % - linear dynamics
    Ax = getax(x_0, dt, N, nx);
    Bu = kron([sparse(1, N); speye(N)], BB);
    Cx = getcx(x_0, u_0, dt, N, nx);
    Aeq = [Ax, Bu];
    leq = [-x0; -Cx];
    ueq = leq;
    % - input and state constraints
    Aineq = speye((N+1)*nx + N*nu);
    lineq = [repmat(xmin, N+1, 1); repmat(umin, N, 1)];
    uineq = [repmat(xmax, N+1, 1); repmat(umax, N, 1)];

    % - OSQP constraints
    A = [Aeq; Aineq; Acbf; A2cbf];
    l = [leq; lineq; lcbf; lcbf2];
    u = [ueq; uineq; ucbf; ucbf2];
    % Create an OSQP object
    prob = osqp;
    % Setup workspace
    prob.setup(P, q, A, l, u, 'warm_start', true);
    % Solve
    res = prob.solve();

    % Check solver status
    if ~strcmp(res.info.status, 'solved')
%        error('OSQP did not solve the problem!')
        t2=-1;
        return
    end
    ctrlx = res.x(1:(N+1)*nx);
    ctrlu = res.x((N+1)*nx+1:(N+1)*nx+N*nu);
    testx = (storagex - ctrlx)'*(storagex - ctrlx);
    testu = (storageu - ctrlu)'*(storageu - ctrlu);
    test = (testx)/(storagex'*storagex);
    if (test)^(0.5)<=10^(-2)&& (testx/((N+1)*nx))^(0.5)<=10^(-4)%Convergence criterion
      break
    end
end % Move for the first step
ctrl = ctrlu(1:nu);
x0 = [dt*x0(4)*cos(x0(3))+x0(1);dt*x0(4)*sin(x0(3))+x0(2);dt*ctrl(1)+x0(3);dt*ctrl(2)+x0(4)];
ctrlu = rewrite(ctrlu, nu, N);
t2=t2+toc;
return

x_0 = newinit(x0, ctrlu, nx, nu, N, dt);
u_0 = transu(ctrlu, nu, N);
impc(:, 2) = x0;
storagex = ctrlx;
storageu = ctrlu;

for i = 1 : (nsim-1)%Iterations for the left time steps

   for j = 1 : K
    abc = tangentline(x_0, r);
    Acbf = cbfmat(abc, nx, nu, dt, gamma1, x0);
    lcbf = -abc(3, :);
    lcbf(1) = [];
    lcbf = lcbf';
    ucbf = inf * ones(N, 1);
    A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0);
    lcbf2 = getlcbf2(abc, dt, gamma1, gamma2);
    ucbf2 = inf * ones(N-1, 1);
    % Cast MPC problem to a QP: x = (x(0),x(1),...,x(N),u(0),...,u(N-1))
    % - quadratic objective
    P = blkdiag( kron(speye(N), Q), QN, kron(speye(N), R) );
    % - linear objective
    q = [repmat(-Q*xr, N, 1); -QN*xr; repmat(-R*ur, N, 1)];
    % - linear dynamics
    Ax = getax(x_0, dt, N, nx);
    Bu = kron([sparse(1, N); speye(N)], BB);
    Cx = getcx(x_0, u_0, dt, N, nx);
    Aeq = [Ax, Bu];
    leq = [-x0; -Cx];
    ueq = leq;
    % - input and state constraints
    Aineq = speye((N+1)*nx + N*nu);
    lineq = [repmat(xmin, N+1, 1); repmat(umin, N, 1)];
    uineq = [repmat(xmax, N+1, 1); repmat(umax, N, 1)];

    % - OSQP constraints
    A = [Aeq; Aineq; Acbf; A2cbf];
    l = [leq; lineq; lcbf; lcbf2];
    u = [ueq; uineq; ucbf; ucbf2];
    % Create an OSQP object
    prob = osqp;
    % Setup workspace
    prob.setup(P, q, A, l, u, 'warm_start', true);
    % Solve
    res = prob.solve();

    % Check solver status
    if ~strcmp(res.info.status, 'solved')
%        error('OSQP did not solve the problem!')
        t2=-1;
        return
    end
    ctrlx = res.x(1:(N+1)*nx);
    ctrlu = res.x((N+1)*nx+1:(N+1)*nx+N*nu);
    x0 = ctrlx(1:nx);
    testx0 = [testx0 x0];
    x_0 = trans(x0, ctrlx, nx, N+1);
    u_0 = transu(ctrlu, nu, N);
    testx = (storagex - ctrlx)'*(storagex - ctrlx);
    testu = (storageu - ctrlu)'*(storageu - ctrlu);
    test = (testx)/(storagex'*storagex);
    if (test)^(0.5)<=10^(-2)&& (testx/((N+1)*nx))^(0.5)<=10^(-4)%Convergence criterion
      break
    end
    storagex = ctrlx;
    storageu = ctrlu;
   end
   ctrl = ctrlu(1:nu);
   x0 = [dt*x0(4)*cos(x0(3))+x0(1);dt*x0(4)*sin(x0(3))+x0(2);dt*ctrl(1)+x0(3);dt*ctrl(2)+x0(4)];
   ctrlu = rewrite(ctrlu, nu, N);
   x_0 = newinit(x0, ctrlu, nx, nu, N, dt);
   u_0 = transu(ctrlu, nu, N);
   impc(:, i+2) = x0;
end
end

% Linerize the CBF constraints (get a, b, c for lines)
function abc = tangentline(xy, r)% x and y from initialize states set, abc are coeeficients for linear equation a*x+b*y+c=0
[xx, ~] = size(xy);%xx=2,yy=N+1
xy(xx,:) = []; % this part should be changed for other case
xy((xx-1),:) = []; % this part should be changed other case
[xx, yy] = size(xy);%xx=2,yy=N+1
xyjiao = zeros(xx, yy);%intersection points
for i = 1 : xx
    for j = 1 : yy
        xyjiao(i, j) = r * xy(i, j) * (1 / (xy(:, j)' * xy(:, j)))^(0.5);%calculate coordinates of intersection points
    end
end
cc = -r^2 * ones(1, yy);
abc = [xyjiao; cc];
end

% Get CBF constraints matrix
function Acbf = cbfmat(abc, nx, nu, dt, gamma, x0)
[~, yy] = size(abc);
Acbfx = zeros((yy-1), yy*nx);
Acbfu = zeros((yy-1), (yy-1)*nu);
for i = 1 : (yy-1)
    Acbfx(i, (i*nx)+1) = abc(1, (i+1));
    Acbfx(i, (i*nx)+2) = abc(2, (i+1));
end
for i = 1 : (yy-1)
    Acbfu(i, ((i-1)*nu+3)) = - (1 - dt * gamma)^(i) * (abc(1, 1) * x0(1, 1) + abc(2, 1) * x0(2, 1) + abc(3, 1));
end
Acbf = [Acbfx Acbfu];
end

% Transfer vector x into matrix
function res = trans(x0, vector, nxx, nyy)%nxx=nx,nyy=N+1
 res = zeros(nxx, nyy);
 res(:,1) = x0;
 for i = 1 : (nyy -1)
     res(:,i+1)= vector(((i)*nxx+1):(i+1)*nxx);
 end
end

% Transfer vector u into matrix
function resu = transu(vector, nxx, nyy)%nxx=nu,nyy=N
 resu = zeros(nxx, nyy);
 for i = 1 : (nyy)
     resu(:,i)= vector(((i-1)*nxx+1):(i)*nxx);
 end
end

% Rewrite u vector
function reu = rewrite(vector, nu, N)
append = vector((N-1)*nu+1:N*nu);
vector(1:nu) = [];
reu = [vector;append];
end

% Get new x_0
function x_0 = newinit(x0, ctrlu, nx, nu, N, dt)
 x_0 = zeros(nx, N+1);
 x_0(:, 1) = x0;
 for i=1 : N
    u0 = ctrlu((i-1)*nu+1:i*nu);
    x0 = [dt*x0(4)*cos(x0(3))+x0(1);dt*x0(4)*sin(x0(3))+x0(2);dt*u0(1)+x0(3);dt*u0(2)+x0(4)];
    x_0(:, i + 1) = x0;
 end
end

% Get AA matrix
function AA = getaa(x0, dt)
 AA = [1 0 -x0(4)*sin(x0(3))*dt cos(x0(3))*dt;0 1 x0(4)*cos(x0(3))*dt sin(x0(3))*dt;0 0 1 0;0 0 0 1];
end
% Get CC matrix
function CC = getcc(x0, x1, u0, dt)
 CC = [x0(4)*sin(x0(3))*x0(3)*dt-x0(4)*cos(x0(3))*dt+x1(1)-x0(1);-x0(4)*cos(x0(3))*x0(3)*dt-x0(4)*sin(x0(3))*dt+x1(2)-x0(2);-u0(1)*dt+x1(3)-x0(3);-u0(2)*dt+x1(4)-x0(4)];
end
% Get Ax matrix
function Ax = getax(x_0, dt, N, nx)
x0 = x_0(:,1);
AA = getaa(x0, dt);
Ax = kron(speye(N+1), -speye(nx)) + kron(sparse(diag(ones(N, 1), -1)), AA);
for i = 1 : (N-1)
   x0 = x_0(:,i+1);
   AA = getaa(x0, dt);
   Ax(nx*(i+1)+1:nx*(i+1)+nx,nx*i+1:nx*i+nx) = AA;
end
end
% Get Cx matrix
function Cx = getcx(x_0, u_0, dt, N, nx)
Cx = zeros(N*nx, 1);
for i = 1 : N
    u0 = u_0(:,i);
    x0 = x_0(:,i);
    x1 = x_0(:,i+1);
    CC = getcc(x0, x1, u0, dt);
    Cx((i-1)*nx+1:(i-1)*nx+nx) = CC;
end
end
% Get A2cbf
function A2cbf = cbfmat2(abc, nx, nu, dt, gamma1, gamma2, x0)
[~, yy] = size(abc);
Acbfx2 = zeros((yy-2), yy*nx);
Acbfx22 = zeros((yy-2), yy*nx);
Acbfu2 = zeros((yy-2), (yy-1)*nu);
for i = 1 : (yy-2)
    Acbfx2(i, (i*nx)+1) = (gamma1-1/dt)*abc(1, (i+1));
    Acbfx2(i, (i*nx)+2) = (gamma1-1/dt)*abc(2, (i+1));
    Acbfx2(i, ((i+1)*nx)+1) = (1/dt)*abc(1, (i+2));
    Acbfx2(i, ((i+1)*nx)+2) = (1/dt)*abc(2, (i+2));
end
for i = 1 : (yy-2)
    Acbfx22(i, (1*nx)+1) = -(1-dt*gamma2)^(i)/dt*abc(1, (1+1));
    Acbfx22(i, (1*nx)+2) = -(1-dt*gamma2)^(i)/dt*abc(2, (1+1));
end
Acbfx2 = Acbfx2 + Acbfx22;
for i = 1 : (yy-2)
    Acbfu2(i, ((i-1)*nu+4)) = - (1 - dt * gamma2)^(i)*(gamma1-1/dt)*(abc(1, 1) * x0(1, 1) + abc(2, 1) * x0(2, 1) + abc(3, 1));
end
A2cbf = [Acbfx2 Acbfu2];
end
% Get lcbf2
function lcbf2 = getlcbf2(abc, dt, gamma1, gamma2)
[~, yy] = size(abc);
lcbf2 = zeros((yy-2),1);
for i = 1 : (yy-2)
    lcbf2(i, 1) = -abc(3, (i+2))/dt-(gamma1-1/dt)*abc(3, (i+1))+(1 - dt * gamma2)^(i)*abc(3, (1+1))/dt;
end
end

author: HybridRobotics

repo: NMPC-DCLF-DCBF