text:
function [mu, s, a, gf] = fit_logn(X, Y, varargin)
%fits lognormal distribution to histogram data
%FIT_LOGN fits lognormal distribution to histogram data (value and count),
%initial guess for mean and standard deviation can be given as optional
%arguments. Return values are mean, standard deviation, scaling factor (for
%non-unit distributions) and goodness of fit (should be close to 1)
%
%Example:
%	[MU, sigma, a, gf] = fitlogn(xdata, ydata [, mean guess [, std guess]])
%
%See also: logn2mean, mean2logn, lognfit, fit_norm, lognpdf, lognstat

p = inputParser;
p.addRequired('X',       @(X) isnumeric(X) & length(X) > 2 );
p.addRequired('Y',       @(X) isnumeric(X) & length(X) > 2 );
p.addOptional('M',   [], @(X) isnumeric(X) & numel(X) == 1 );
p.addOptional('SIG', [], @(X) isnumeric(X) & numel(X) == 1 & X > 0);
p.parse(X, Y, varargin{:});

warning('off', 'curvefit:fit:complexXusingOnlyReal');
warning('off', 'MATLAB:singularMatrix');

X = X(:);
Y = Y(:);
if any(X <= 0)
	Y(X <= 0) = [];
	X(X <= 0) = [];
end

Xn = log(X);
Yn = Y.*X;
Xs = (Xn(1:end-1)+diff(Xn)/2);
Ys = Yn(1:end-1);

if ~isempty(p.Results.M) && ~isempty(p.Results.SIG)
	[MU, S] = mean2logn(p.Results.M, p.Results.SIG);
elseif ~isempty(p.Results.M) && isempty(p.Results.SIG)
	[MU, S] = mean2logn(p.Results.M, 0.1*p.Results.M^2);
else
	MU     = sum( Ys.*Xs )/sum( Ys );
	S      = sqrt(abs(sum( Ys.*(Xs-MU).^2 )/sum( Ys )));
end

% Initial Fit
[mu, s, a] = fit_norm(Xn, Yn, MU, S);
gf         = gof( Y, lognpdf(X, mu, s)*a );

% Reset Starting Parameters
if gf < 0.9999
	nsteps = min([50 max([3 ceil(10^(1-gf))])]);
	s_vec  = linspace(max([0 (2-nsteps)*S]), nsteps*S, nsteps);
	s_vec(s_vec == 0) = eps;
	cfs    = cell(size(s_vec));
	gofs   = zeros(size(s_vec));
	for j = 1:length(s_vec)
		[cfs{j}(2), cfs{j}(3), cfs{j}(1)] = fit_norm(Xn, Yn, MU, s_vec(j));
		gofs(j) = gof( Y, lognpdf(X,cfs{j}(2),cfs{j}(3))*cfs{j}(1) );
	end
	[~, i]         = max(gofs);
	[mu, s, a, gf] = fit_norm(Xn, Yn, cfs{i}(2), cfs{i}(3));
end

warning('on', 'curvefit:fit:complexXusingOnlyReal');
warning('on', 'MATLAB:singularMatrix');
end

function G = gof(data, fitdata, varargin)
%calculates goodness-of-fit value R^2
%GOF calculates goodness-of-fit value R^2 from two data sets with an
%optional weight input
%
%Example:
%    gof(original values, fitted values [, weight])
%
% See also: fit
%

	idx		=	~isnan(data) & ~isnan(fitdata);
	data	=	data(idx);
	fitdata	=	fitdata(idx);
	if nargin > 2 && isnumeric(varargin{1})
		weights		=	varargin{1}(idx);
	else
		weights		=	ones(size(data));
	end

	G = 1 - sum(weights.*(data-fitdata).^2)/sum(weights.*(data - mean(data)).^2);

end

% Copyright 2009-2013 Alexandra Heidsieck <aheidsieck@tum.de>,
%                     IMETUM, Technische Universitaet Muenchen


text:
function [f_unk_1]=rlse(yk_1,hk_1,G_un,rk)

% Estimate sk_1: A matrix with dimension rxr
Sk_1=inv(G_un'*inv(rk)*G_un);

% Recusive solution for unknown excitation
f_unk_1=Sk_1*G_un'*inv(rk)*(yk_1-hk_1);

author: Sable

repo: mcbench-benchmarks


text:
function fmarginal = add_ev_to_dmarginal(fmarginal, evidence, ns)
% ADD_EV_TO_DMARGINAL 'pump up' observed nodes back to their original size.
% fmarginal = add_ev_to_dmarginal(fmarginal, evidence, ns)
%
% We introduce 0s into the array in positions which are incompatible with the evidence.

dom = fmarginal.domain;
odom = dom(~isemptycell(evidence(dom)));
vals = cat(1, evidence{odom});
index = mk_multi_index(length(dom), find_equiv_posns(odom, dom), vals);
T = 0*myones(ns(dom));
ens = ns(:)';
ens(odom) = 1;
T(index{:}) = myreshape(fmarginal.T, ens(dom));
fmarginal.T = T;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% Calculate gridness score for an autocorrelogram
%
% Calculates a gridness score by expanding a circle around the centre field and
% calculating a correlation value of the expanded circle with it's rotated versions.
% The expansion is done up until the smallest side of the autocorrelogram.
% Can also calculate grid statistics.
%
% Gridness score value by itslef is calculated as a maximum over sliding mean
% of expanded circles. The widtg of the sliding window is given by variable
% numGridnessRadii. This is done in order to keep the values the same with
% historical development of gridness score.
%
% NB! The function assumes that the provided input is indeed an autocorrelgoram.
% This means that the center pixel must have the maximum value. This must be true
% for both square and rectangular autocorrelogram. This will be important if you
% manually create square matrix from a rectangular one. Imagine you have a rectangular
% matrix aCorr of size 309x305, aCorr(155, 153) == 1. If you make it rectangular
% like  this `aSquare = aCorr(1:305, :);`, then the centre pixel will not have
% the highest values, i.e. `aSquare(153, 153) ~= 1`. You should do it like this
% `aSquare = aCorr(3:309-2, :)`, then `aSquare(153, 153)` yields `1`.
%
%  USAGE
%   [score, <stats>] = analyses.gridnessScore(aCorr, <options>)
%   aCorr       A 2D autocorrelogram. aCorr can be square or rectangular. See the note (NB) above!
%   <options>   Optional list of property-value pairs (see table below)
%
%   ==============================================================================================
%    Properties         Values
%   ----------------------------------------------------------------------------------------------
%    'minOrientation'   Value of minimal difference of inner fields orientation (in degrees). If
%                       there are fields that differ in orientation for less than minOrientation,
%                       then only the closest to the centre field are left. Default value is 15.
%   'debug'             True or False. If set to True, the function produces some debug output.
%   ==============================================================================================
%   score       Gridness score. Ranges from -2 to 2. 2 is more a theoretical bound for a perfect grid.
%               More practical value is around 1.3.
%   stats       If this variable is requested, then it is a structure with the following statistics:
%       spacing         3-element vector with distances from the centre field to neighbour fields.
%       orientation     3-element vector with orientations between the centre field and neighbour fields.
%       ellipse         Ellipse fitted to the grid. Contains the centre, radii and orientation in
%                       radians, stored as [Cx, Cy, Rx, Ry, theta].
%       ellipseTheta    Radius of the ellipse in degrees wrapped in range [0..180].
%
function [gscore, varargout] = gridnessScore(aCorr, varargin)
    nout = max(nargout, 1) - 1;
    inp = inputParser;
    defaultMinOrientation = 15;
    gridStat.orientation = [];
    gridStat.spacing = [];
    gridStat.ellipse = [];
    gridStat.ellipseTheta = nan;

    % input argument check functions
    checkDScalar = @(x) helpers.isdscalar(x, '>0');

    % fill input parser object
    addRequired(inp, 'aCorr');
    addParameter(inp, 'minOrientation', defaultMinOrientation, checkDScalar);
    addParameter(inp, 'debug', false);

    parse(inp, aCorr, varargin{:});

    % get parsed arguments
    minOrientation = inp.Results.minOrientation;
    isDebug = inp.Results.debug;

    halfSize = ceil(size(aCorr)/2);
    half_height = halfSize(1);
    half_width = halfSize(2);
    aCorrRad = min(halfSize);
    aCorrSize = size(aCorr);

    if aCorrSize(1) == 1 || aCorrSize(2) == 1
        gscore = nan;
        if nout > 0
            varargout{1} = gridStat;
            varargout{2} = 0;
            varargout{3} = nan(6, 2);
            varargout{4} = 0;
        end
        return;
    end

    % contourc is efficient if aCorr is normalized
    maxValue = max(max(aCorr));
    if maxValue ~= 1
        aCorr = aCorr / maxValue;
    end

    cFieldRadius = findCentreRadius(aCorr, half_width, half_height);
    if isDebug
        fprintf('Center radius is %f\n', cFieldRadius);
    end
    if cFieldRadius == 0 || cFieldRadius == 1 || cFieldRadius == -1 || cFieldRadius >=min(halfSize)
        gscore = nan;
        if nout > 0
            varargout{1} = gridStat;
            varargout{2} = 0;
            varargout{3} = nan(6, 2);
            varargout{4} = 0;
        end
        return;
    end

    % Meshgrid for expanding circle
    [rr, cc] = meshgrid(1:size(aCorr, 2), 1:size(aCorr, 1));

    % Define iteration radius step size for the gridness score
%     if cFieldRadius>=aCorrRad  % modified by Weijian Zong,20201012
%     cFieldRadius=aCorrRad;  % modified by Weijian Zong,20201012
%     else% modified by Weijian Zong,20201012
%     end % modified by Weijian Zong,20201012
    radSteps = cFieldRadius:aCorrRad;
    radSteps(1) = [];
    numSteps = length(radSteps);

    GNS = zeros(numSteps, 2);
    rotCorr = zeros(1, 5);
    rotAngles_deg = 30*(1:5);

    % aCorr is rotated outside the loop for speed
    rotatedCorr = cell(1, length(rotAngles_deg));
    for i = 1:length(rotAngles_deg)
        rotatedCorr{i} = imrotate(aCorr, rotAngles_deg(i), 'bilinear', 'crop');
    end

    mainCircle = sqrt((cc - half_height).^2 + (rr - half_width).^2);
    innerCircle = mainCircle > cFieldRadius;

    % Define expanding ring of autocorrellogram and do x30 correlations
    for i = 1:numSteps
        ind = (innerCircle & (mainCircle < radSteps(i)));
        tempCorr = reshape(aCorr(ind), 1, [])';
        for j = 1:5
            rotatedCircle = reshape(rotatedCorr{j}(ind), 1, [])';
            rotCorr(j) =  corr(tempCorr, rotatedCircle);
            if isDebug
                fprintf('Step %u, angle %u, corr value %f\n', i-1, rotAngles_deg(j), rotCorr(j));
            end
        end
        GNS(i, 1) = min(rotCorr([2, 4])) - max(rotCorr([1, 3, 5]));
        GNS(i, 2) = radSteps(i);
    end

    % Find the biggest gridness score and radius
    [~, gscoreLoc] = max(GNS(:, 1));
    % See function help about numGridnessRadii
    numGridnessRadii = 3;
    numStep = numSteps - numGridnessRadii;
    if numStep < 1
        numStep = 1;
    end

    if numStep == 1
        gscore = nanmean(GNS(:, 1));
    else
        meanGridnessArray = zeros(numStep, 1);
        for ii = 1:numStep
            meanGridnessArray(ii) = nanmean(GNS(ii:ii + numGridnessRadii-1, 1));
        end

        [gscore, gInd] = max(meanGridnessArray);
        gscoreLoc = gInd + (numGridnessRadii-1)/2;
    end

    varargout{4} = radSteps(gscoreLoc);

    % Return if we do not need to calculate grid statistics
    if nout < 1
        return;
    end

    %% Calculate gridness score statistics
    bestCorr = (mainCircle < radSteps(gscoreLoc) * 1.25) .* aCorr;
    regionalMaxMap = imregionalmax(bestCorr, 4);
    se = strel('square', 3);
    im2 = imdilate(regionalMaxMap, se); % dilate map to eliminate fragmentation
    cc = bwconncomp(im2, 8);
    stats = regionprops(cc, 'Centroid');

    if length(stats) < 5
        warning('BNT:numFields', 'Not enough inner fields has been found. Can''t calculate grid properties');

        varargout{1} = gridStat;
        varargout{2} = cFieldRadius;
        varargout{3} = nan(6, 2);
        varargout{4} = radSteps(gscoreLoc);
        return;
    end

    allCoords = [stats(:).Centroid];
    centresOfMass(:, 1) = allCoords(1:2:end);
    centresOfMass(:, 2) = allCoords(2:2:end);

    % Calculate orientation for each field relative to the centre field
    orientation = (atan2(centresOfMass(:, 2) - half_height, centresOfMass(:, 1) - half_width)); % atan2(Y, X)
    peaksToCentre = sqDistance(centresOfMass', [half_width half_height]');
    zeroInd = find(orientation == 0, 1);
    orientation(zeroInd) = []; % remove zero value, so that we do not have a side effect with minOrientation
    stats(zeroInd) = [];
    peaksToCentre(zeroInd) = [];
    centresOfMass(zeroInd, :) = [];

    % filter fields that have similar orientation
    orientDistSq = CircStat2012a.circ_dist2(orientation);
    closeFields = abs(orientDistSq) < deg2rad(minOrientation);
    [rows, cols] = size(closeFields);
    closeFields(1:(rows+1):rows*cols) = 0; % assign zero to diagonal elements
    closeFields(tril(true(rows))) = 0; % assign zero to lower triangular of a matrix. Matrix is
                                       % symmetric and we do not need these values.
    [rows, cols] = find(closeFields); % find non-empty elements, they correspond to indices of close fields
    if ~isempty(rows)
        indToDelete = zeros(1, length(rows));
        for i = 1:length(rows)
            % fieldPeaks = [fields([rows(i) cols(i)]).peakX; fields([rows(i) cols(i)]).peakY];
            % peaksToCentre = sqDistance(fieldPeaks, [half_width; half_height]);
            if peaksToCentre(rows(i)) > peaksToCentre(cols(i))
                indToDelete(i) = rows(i);
            else
                indToDelete(i) = cols(i);
            end
        end
        indToDelete = unique(indToDelete);
        stats(indToDelete) = [];
        peaksToCentre(indToDelete) = [];

        if length(stats) < 4
            warning('BNT:numFields', 'Not enough inner fields has been found. Can''t calculate grid properties');

            varargout{1} = gridStat;
            varargout{2} = cFieldRadius;
            varargout{3} = nan(6, 2);
            varargout{4} = radSteps(gscoreLoc);
            return;
        end

        allCoords = [stats(:).Centroid];
        clear centresOfMass;
        centresOfMass(:, 1) = allCoords(1:2:end);
        centresOfMass(:, 2) = allCoords(2:2:end);
    end

    % % get fields peak coordinates
    % fieldPeaks = zeros(length(stats), 2);
    % for i = 1:length(stats)
    %     [~, maxInd] = max(bestCorr(stats(i).PixelIdxList));
    %     fieldPeaks(i, :) = stats(i).PixelList(maxInd, :);
    % end
    % % fieldPeaks = [fields(:).peakX; fields(:).peakY]'; %
%     peaksToCentre = sqDistance(centresOfMass', [half_width half_height]');
    [~, sortInd] = sort(peaksToCentre);
    stats = stats(sortInd);
    centresOfMass = centresOfMass(sortInd, :);

    % leave only 6 closest neighbours (if available)
    if length(stats) > 5
%         stats = stats(1:6);
        centresOfMass = centresOfMass(1:6, :);
    else
%         stats = stats(1:end);
        centresOfMass = centresOfMass(1:end, :);
    end

    % centresOfMass = [fields(:).x; fields(:).y]';
    % Calculate orientation for each field relative to the centre field
    orientation = rad2deg(atan2(centresOfMass(:, 2) - half_height, centresOfMass(:, 1) - half_width)); % atan2(Y, X)

    % Calculate distances between centre of masses for each field and the centre field
    spacing = sqrt((centresOfMass(:, 1) - half_width).^2 + (centresOfMass(:, 2) - half_height).^2);

%     % Plot grid polygon points
%     figure, plot.colorMap(bestCorr), hold on;
%     plot(centresOfMass(:, 1), centresOfMass(:, 2), '+k', 'markersize', 8);

    ell = general.fitEllipse(centresOfMass(:, 1), centresOfMass(:, 2));
    ellipseTheta = rad2deg(general.wrap(ell(end)) + pi);
%     drawEllipse(ell, 'linewidth', 2, 'color', [1 1 1]);

    % Determine axes orientation, spacing and deviation
    [~, bBC] = sort(abs(orientation));
    [~, bBC2] = sort(abs(orientation - orientation(bBC(1))));

    % leave only three values, because autocorrelogram is symmetric
    orientation = orientation(bBC2(1:3));
    spacing = spacing(bBC2(1:3));
    [orientation, orientSortInd] = sort(orientation);

    spacing = spacing(orientSortInd);

    gridStat.orientation = orientation;
    gridStat.spacing = spacing;
    gridStat.ellipse = ell;
    gridStat.ellipseTheta = ellipseTheta;

    varargout{1} = gridStat;
    varargout{2} = cFieldRadius;
    varargout{3} = centresOfMass;
    varargout{4} = radSteps(gscoreLoc);
end

function D = sqDistance(X, Y)
    D = bsxfun(@plus,dot(X,X,1)',dot(Y,Y,1))-2*(X'*Y);
end

function cFieldRadius = findCentreRadius(aCorr, half_width, half_height)
    cFieldRadius = 0;

    % Search for fields only around the centre
    peakCoords = [half_width, half_height];
    [~, fields] = SpatialTuning_BNT.placefieldAdaptive(aCorr, 'minPeak', 0, 'minBins', 2, ...
        'peakCoords', peakCoords);
    if isempty(fields)
        return;
    end

    peakLoc = [fields(:).peakX; fields(:).peakY]';

    % get all distances and check two minimums of them
    allDistances = sqDistance(peakLoc', [half_width half_height]'); % point should be in format [x y]
    [~, sortIndices] = sort(allDistances);
    if length(sortIndices) >= 2
        % this is a bit leagacy code. By using peakCoords with placefieldAdaptive, we
        % should always get just a single field. Keeping it as I did not test the version without it properly.
        twoMinIndices = sortIndices(1:2);

        if abs(allDistances(twoMinIndices(1)) - allDistances(twoMinIndices(2))) < 2
            % two fields with close middle points. Let's select one with minimum square
            [~, minInd] = min([fields(twoMinIndices).area]);
            closestFieldInd = twoMinIndices(minInd);
        else
            closestFieldInd = twoMinIndices(1); % get the first minimum
        end
    else
        closestFieldInd = sortIndices(1);
    end
    cFieldRadius = floor(sqrt(fields(closestFieldInd).area / pi));
end

author: kavli-ntnu


text:
function fgrabm_init
%FGRABM_INIT FrameGrabM Initialization routine
%   FGRABM_INIT is the first activity that must be performed to start up
%   the underlying video capture library.  This produces a global variable
%   FGRABM that will be used by other FrameGrabM functions.  Most other
%   FrameGrabM functions will call FGRABM_INIT automatically if needed.
%
%   Version 0.8 - 06 March 2012
global FGRABM

% Declare our Java classpath:
import com.academiken.framegrabm.util.*;

% Create our persistent structure:
if isempty(FGRABM)
    FGRABM = struct();
end

% Instanciate the Java object:
if ~isfield(FGRABM, 'class')
    FGRABM.class = GrabCore;
end

% Set the default device:
if ~isfield(FGRABM, 'defaultDevice')
    FGRABM.defaultDevice = 1;
end

% Prune off formats field so formats can be refreshed:
if isfield(FGRABM, 'formats')
    FGRABM = rmfield(FGRABM, 'formats');
end

% Initialize:
FGRABM.class.init()


text:
% Weighted Prediction Error for dereverberation
% reference: Nakatani, Tomohiro, et al. "Speech dereverberation based on
% variance-normalized delayed linear prediction." IEEE TASLP 18.7 (2010)
% ZitengWANG@201903

clear all
addpath('..\STFT\')
addpath('..\Simulation\')
addpath('..\Simulation\RIR-Generator\')

%% simulation start
flatStart = 1;
prefix = '';   % for saving file

speechDir = '..\Simulation\Data\';
speechFile = 'fajw0_sa1.wav';
saveDir = 'GeneratedData\';
if ~exist(saveDir)
    mkdir(saveDir)
end
speech = audioread([speechDir speechFile]);

% configuration
cfg = [];
cfg.fs = 16000;                     % sampling rate
cfg.room = [6 5 3];                 % room dimension (m)
cfg.T60 = 0;                      % reverberation time (s)

cfg.Nch = 6;
cfg.micCenter = [2 3 1.5];          % array center (m)
cfg.micCoordinate =[0.0425,0.0,0.0;
        0.02125,0.03680608,0.0;
        -0.02125,0.03680608,0.0;
        -0.0425,0.0,0.0;
        -0.02125,-0.03680608,0.0;
        0.02125,-0.03680608,0.0;];  % microphone array coordinates

cfg.az = 180;
cfg.el = 0;
cfg.dist = 3;

cfg.SNR = 20;
if cfg.SNR ~= inf               % inf means no noise
    cfg.noiseType = 'white';    % choice {'white' 'diffuse' 'recorded'}
    if strcmp(cfg.noiseType, 'recorded')
        % check first the noise is longer than speech!
        cfg.noiseFile = '';
    end
end

cfg.SIR = inf;

% setup room and collect data
setup_room
setup_noise

%% offline processing start
% the following parts are specific to the algorithm
Nfft = 512;

Y = stft_multi_2(y, Nfft);
[Nframe, Nbin, Nch] = size(Y);

iterMax = 5;
wpe.delay = 3;
wpe.taps = 10;
XeEst = Y;
for bin = 1:Nbin
    % get audio context first
    YbarFrm = cell(Nframe, 1);
    for frm = (wpe.taps + wpe.delay):Nframe
        Ybar = squeeze(Y(frm - wpe.delay,bin,:));
        for tap=1:wpe.taps-1
            Ybar = [Ybar; squeeze(Y(frm - wpe.delay - tap,bin ,:))];
        end
        YbarFrm{frm} = Ybar;
    end

    % iterate
    X = squeeze(Y(:,bin,:));
    for iter = 1:iterMax
        % calculate mean signal power
        Xpow = mean(abs(X).^2, 2);
        XpowMax = max(Xpow);
        % calculate correlation matrix and correlation vector
        R = 0;
        P = 0;
        for frm = (wpe.taps + wpe.delay):Nframe
            Ybar = YbarFrm{frm};
            YbarTmp = Ybar / max(Xpow(frm), 1e-10*XpowMax);
            R = R + YbarTmp * Ybar' ;
            P = P + YbarTmp * squeeze(Y(frm,bin,:))';
        end
        hWPE = R \ P;

        % apply the filter
        for frm = (wpe.taps + wpe.delay):Nframe
            X(frm, :) = squeeze(Y(frm,bin,:)) - hWPE'*YbarFrm{frm};
        end
    end
    XeEst(:,bin,:) = permute(X,[1,3,2]);
end

xeEst = istft_multi_2(XeEst, length(speech));
audiowrite([saveDir prefix 'WPE' postfix '.wav'], xeEst, fs);

author: ZitengWang

repo: MASP

save_path: github-repos/MATLAB/ZitengWang-MASP


text:
function [ lo, hi ] = p02_box ( m )

%*****************************************************************************80
%
%% P02_BOX returns a bounding box for problem 02.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Output, real LO(M), HI(M), the low and high corners of the box.
%
  center = [ 0.0, 0.0 ];
  r1 = 1.0;
  r2 = 0.4;

  lo(1:m) = [ center(1) - r1, center(2) - r1 ];
  hi(1:m) = [ center(1) + r1, center(2) + r1 ];

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x, infos] = spg_nmf(V, rank, in_options)
% Stochastic projected gradient for non-negative matrix factorization (SMU-NMF) algorithm.
%
% Inputs:
%       matrix      V
%       rank        rank
%       options     options
% Output:
%       w           solution of w
%       infos       information
%
%
% This file is part of NMFLibrary.
%
% Created by H.Kasai on Mar. 28, 2017
%
% Change log:
%
%       Mar. 14, 2018 (Hiroyuki Kasai): Fixed algorithm.
%
%       May. 20, 2019 (Hiroyuki Kasai): Added initialization module.
%
%       Jul. 12, 2022 (Hiroyuki Kasai): Modified code structures.
%

    % set dimensions and samples
    [m, n] = size(V);

    % set local options
    local_options = [];
    local_options.W_sub_mode          = 'STD';
    local_options.H_sub_mode          = 'LS';
    local_options.accel               = false;
    local_options.ls                  = false;
    local_options.precon              = false;
    local_options.h_repeat            = 1;
    local_options.rep_mode            = 'fix';
    local_options.robust              = false;

    % check input options
    if ~exist('in_options', 'var') || isempty(in_options)
        in_options = struct();
    end
    % merge options
    options = mergeOptions(get_nmf_default_options(), local_options);
    options = mergeOptions(options, in_options);

    % initialize factors
    init_options = options;
    [init_factors, ~] = generate_init_factors(V, rank, init_options);
    Wt = init_factors.W;
    H = init_factors.H;
    R = init_factors.R;

    % determine sub_mode
    if options.accel
        options.H_sub_mode = 'ACC';
    else
        if options.ls
            options.H_sub_mode = 'LS';
        else
            options.H_sub_mode = 'STD';
        end
        options.h_repeat = 1;
    end

    % permute samples
    if options.permute_on
        perm_idx = randperm(n);
    else
        perm_idx = 1:n;
    end
    V = V(:,perm_idx);
    H = H(:,perm_idx);

    % initialize
    method_name = sprintf('SPG-NMF (%s,%s)', options.W_sub_mode, options.H_sub_mode);
    epoch = 0;
    grad_calc_count = 0;
    nesterov_alpha = 1;

    if options.verbose > 0
        fprintf('# %s: started ...\n', method_name);
    end

    if strcmp(options.rep_mode, 'adaptive')
        rhoh = 1+(m+m*rank)/(1*(rank+1));
        alpha = 2;
        delta = 0.01;
    end

    % store initial info
    clear infos;
    [infos, f_val, optgap] = store_nmf_info(V, Wt, H, R, options, [], epoch, grad_calc_count, 0);

    if options.verbose > 1
        fprintf('SPG-NMF (%s,%s): Epoch = 0000, cost = %.16e, optgap = %.4e\n', options.W_sub_mode, options.H_sub_mode, f_val, optgap);
    end

    % set start time
    start_time = tic();

    % main outer loop
    while true

        % check stop condition
        [stop_flag, reason, max_reached_flag] = check_stop_condition(epoch, infos, options);
        if stop_flag
            display_stop_reason(epoch, infos, options, method_name, reason, max_reached_flag);
            break;
        end

        cnt = 0;

        % main inner loop
        for t = 1 : options.batch_size : n - 1
            cnt = cnt + 1;

            % retrieve vt and ht
            vt = V(:,t:t+options.batch_size-1);
            ht = H(:,t:t+options.batch_size-1);

            % uddate ht
            Wtv = Wt.' * vt;
            WtW = Wt.' * Wt;
            if strcmp(options.H_sub_mode, 'ACC')
                if strcmp(options.rep_mode, 'adaptive')
                    gamma = 1;
                    eps0 = 1;
                    j = 1;
                    rhoh_alpha = rhoh*alpha;
                    %while j <= floor(1+rhoh*alpha) &&  gamma >= delta*eps0
                    while j <= rhoh_alpha && gamma >= delta*eps0
                        ht0 = ht;
                        ht = ht .* (Wtv) ./ (WtW * ht);
                        ht = ht + (ht<eps) .* eps;
                        if j == 1
                            eps0 = norm(ht0-ht);
                        end
                        gamma = norm(ht0-ht);
                        j = j+1;
                    end
                else
                    for ii=1:options.h_repeat
                        ht = ht .* (Wtv) ./ (WtW * ht);
                        ht = ht + (ht<eps) .* eps;
                    end
                end
            elseif strcmp(options.H_sub_mode, 'LS')
                ht = calc_nls_nmf(vt, Wt, 1e-16);
                ht = ht + (ht<eps) * eps;
            else
                ht = ht .* (Wtv) ./ (WtW * ht);
                ht = ht + (ht<eps) .* eps;
            end

            % update W
            if strcmp(options.W_sub_mode, 'Precon')
                Wt = projection_precon_mnls(vt, ht', Wt);
            elseif strcmp(options.W_sub_mode, 'Nesterov')
                options.nesterov_maxit = 1;
                options.nesterov_func_type = 'smooth';
                %options.nesterov_func_type = 'stochastic';
                %options.nesterov_func_type = 'strong_alpha_beta';

                [Wt, apg_iter, nesterov_alpha] = nesterov_mnls(vt, ht', Wt, nesterov_alpha, options.nesterov_maxit, options.nesterov_func_type);
                %nesterov_alpha
            else
                Wt = projection_mnls(vt, ht', Wt);
            end
            Wt = Wt + (Wt<eps) .* eps;

            % store new h
            H(:,t:t+options.batch_size-1) = ht;

            grad_calc_count = grad_calc_count + m * options.batch_size;
        end

        % measure elapsed time
        elapsed_time = toc(start_time);

        % update epoch
        epoch = epoch + 1;

        % store info
        infos = store_nmf_info(V, Wt, H, R, options, infos, epoch, grad_calc_count, elapsed_time);

        % display info
        display_info(method_name, epoch, infos, options);

    end

    x.W = Wt;
    x.H(:,perm_idx) = H;

end

author: hiroyuki-kasai

repo: NMFLibrary


text:
function [varargout]=meshView(varargin)

% function [varargout]=meshView(varargin)
% ------------------------------------------------------------------------
% 2018/01/23 Updated to allow for subfigure plotting
% 2018/04/16 Added control of direction of cutting (X,Y, or Z) and also
% what side is viewed/cut away.
% 2019/09/07 Fixed bug in relation to single element/face input
% 2019/09/07 Updated handling of cutting direction
% ------------------------------------------------------------------------
%% Parse input

switch nargin
    case 1
        meshStruct=varargin{1};
        optionStruct=[];
    case 2
        meshStruct=varargin{1};
        optionStruct=varargin{2};
    otherwise
        error('Wrong number of input arguments');
end

defaultOptionStruct.hFig=[];
defaultOptionStruct.numSliceSteps=25; %Number of animation steps
defaultOptionStruct.cMap=[];
defaultOptionStruct.cutDir=2;
defaultOptionStruct.cutSide=1;
defaultOptionStruct.faceAlpha1=0.2;
defaultOptionStruct.faceAlpha2=1;
defaultOptionStruct.lightWeightPlot=1;
defaultOptionStruct.edgeWidth=1;
defaultOptionStruct.edgeColor='k';
[optionStruct]=structComplete(optionStruct,defaultOptionStruct,0);

%% Get control parameters

figHandles=optionStruct.hFig;
if ~isempty(figHandles)
    hFig=figHandles(1);
    if numel(figHandles)==2
        hSub=figHandles(2);
    else
        hSub=[];
    end
else
    hFig=[];
    hSub=[];
end

numSliceSteps=optionStruct.numSliceSteps;
cMap=optionStruct.cMap;
cutDir=optionStruct.cutDir;
cutSide=optionStruct.cutSide;
faceAlpha1=optionStruct.faceAlpha1;
faceAlpha2=optionStruct.faceAlpha2;
lightWeightPlot=optionStruct.lightWeightPlot;
edgeWidth=optionStruct.edgeWidth;
edgeColor=optionStruct.edgeColor;
fontSize=15;

%% Access mesh data

%                 nodes: [415×3 double]
%         facesBoundary: [320×3 double]
%        boundaryMarker: [320×1 double]
%                 faces: [8144×3 double]
%              elements: [2036×4 double]
%     elementMaterialID: [2036×1 double]
%        faceMaterialID: [8144×1 double]
%        loadNameStruct: [1×1 struct]

E=meshStruct.elements;
V=meshStruct.nodes;

if isfield(meshStruct,'facesBoundary')
    Fb=meshStruct.facesBoundary;
    Cb=meshStruct.boundaryMarker;
else
    F=meshStruct.faces;
    indBoundary=tesBoundary(F);
    Fb=F(indBoundary,:);
    Cb=ones(size(Fb,1),1);
end

if isfield(meshStruct,'elementData')
    CE=meshStruct.elementData;
    numColormapLevels=250;
elseif isfield(meshStruct,'elementMaterialID')
    CE=meshStruct.elementMaterialID;
    numColormapLevels=numel(unique(CE(:)));
else
    CE=ones(size(E,1),1);
    numColormapLevels=250;
end

if isempty(cMap)
    cMap=gjet(numColormapLevels);
end

%%
% prepare figure
if isempty(hFig)
    hFig=cFigure; %Store figure handle
    title('Cut view of mesh','FontSize',fontSize);
end

figure(hFig);
if ~isempty(hSub)
    subplot(hSub);
end
hold on;

if ~isempty(Fb)
    gpatch(Fb,V,0.5*ones(1,3),'none',faceAlpha1);
end

hp=gpatch(Fb,V,Cb,edgeColor,faceAlpha2,edgeWidth); %Graphics object to vary property of during animation

camlight headlight;
axisGeom(gca,fontSize);

if abs(max(CE(:))-min(CE(:)))>eps(1)
    caxis([min(CE(:)) max(CE(:))]);
end

colormap(gca,cMap);
if ~isfield(meshStruct,'elementData')
    icolorbar;
else
    colorbar;
end

%%
% Set up animation

[VE]=patchCentre(E,V);
XE=VE(:,cutDir);
animStruct.Time=linspace(0,1,numSliceSteps); %Time vector
cutLevel=linspace(min(XE)-max(eps(XE)),max(XE)+max(eps(XE)),numSliceSteps); %Property to set

for q=1:1:numSliceSteps %Step through time
    cutLevelNow=cutLevel(q); %The current cut level

    if cutSide==1
        logicCutView=XE>=cutLevelNow;
    elseif cutSide==-1
        logicCutView=XE<=cutLevelNow;
    end

    [Fs,Cs]=element2patch(E(logicCutView,:),CE(logicCutView));

    if lightWeightPlot==1
        [indBoundary]=tesBoundary(Fs);
        Fs=Fs(indBoundary,:);
        Cs=Cs(indBoundary,:);
    end

    %Set entries in animation structure
    animStruct.Handles{q}=[hp hp]; %Handles of objects to animate
    animStruct.Props{q}={'Faces','CData'}; %Properties of objects to animate
    animStruct.Set{q}={Fs,Cs}; %Property values for to set in order to animate
end

%Add animation layer
anim8(hFig,animStruct);
set(hFig.UserData.anim8.sliderHandles{1},'Value',round(numSliceSteps/2)); %Set to middle
gdrawnow;

%% Collect output

switch nargout
    case 1
        varargout{1}=hFig;
    case 2
        varargout{1}=hFig;
        varargout{2}=hp;
end

%%
% _*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:
%This is an example of using the matlab scripts for drawing MEMS fabrication process cross-section. Only support deposit, lithography, and etch(isotropic and anisotropic is not differentied here.). Materials can be used is Si, SiNi,SiO2, metal, PR(photoresistor).

% basic syntax:
% 1. deposit ('materialName', layerThickness), where materialName is one of Si, SiNi, SiO2, metal, PR. layerThickness is an integer, normally it is between 1 and 20.
% 2. liga(mask), where mask is an series of char combined by ' ' and '-', this function lithography photoresistor, so it should used after an layer of PR is deposited.
% 3. etch('materialName'): material which is appointed will be etched from top to bottom until meet an fully cover other material layer.
% 4. removepr('PR'): remove all the photoresistor.
% 5. drawlayers(); draw current result.
% DATE: 30/8/2007
% muwn.gu@gmail.com
%========donot modify the head
clear all;
clear global;
loadlayerprofile('MEMS');
global layerMatrix;
global Layerprofile;
figure(1);
%===================

%*******************replace below to your own process
mask1 = ' --     ---     --      -  ';	% metal 1
mask2 = '    -----------            ';   % SiO2 beam
mask3 = '  ---------------          ';
mask4 = '                       --- ';  % SiO2
mask5 = '                      ---  ';	% metal 2

deposit ('Si', 10);
deposit('SiO2', 20);
deposit('SiNi', 5);

% electrode
deposit('metal', 5); % metal 1 thick
deposit('PR', 10);
liga(mask1);
etch ('metal');
removepr('PR');
drawlayers();

% beams
deposit('SiO2', 5);
deposit('PR', 10);
liga(mask2);
etch('SiO2');
removepr('PR');
drawlayers();

deposit('PoSi', 10);
deposit('PR', 10);
liga(mask3);
etch('PoSi');
removepr('PR');
drawlayers();
etch('SiO2');

% electrode
deposit('metal', 2); % metal 2 thin
deposit('PR', 10);
liga(mask4);
etch ('metal');
removepr('PR');
drawlayers();
%**********************************************

%=================donot modify below===================
axis ([15 80 0 Layerprofile.currentY])


text:
function Datacube_conved = Convolution(Datacube,W,ConvFilter)
% Given W, performing convolution on Datacube with zero padding

[Width,Height,Channel] = size(Datacube);
num_filters = size(W,2);
Datacube_conved = cell(num_filters,1);
magSize = (ConvFilter.PatchSize-1)/2;
magChannel = (ConvFilter.Channel-1)/2;
Tempcube_conved = zeros(Width+ConvFilter.PatchSize-1, Height+ConvFilter.PatchSize-1, Channel+ConvFilter.Channel-1);
Tempcube_conved((magSize+1):end-magSize,(magSize+1):end-magSize,(magChannel+1):end-magChannel) = Datacube;
X = im2colstep(Tempcube_conved,[ConvFilter.PatchSize,ConvFilter.PatchSize,ConvFilter.Channel]);
mu = mean(X,2);
% mu = mean(X,1);
X = bsxfun(@minus, X, mu);
for i=1:num_filters
    Datacube_conved{i} = reshape(W(:,i)'*X,[Width,Height,Channel]);
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function x = stft_iw(Sx, opt)
% ISTFT  Inverse short-time Fourier transform
%
% Very closely based on Steven Schimmel's stft.m and istft.m from
% his SPHSC 503: Speech Signal Processing course at Univ. Washington.

if nargin<2, opt = struct(); end

if ~isfield(opt, 'window'), opt.window = round(size(Sx,2)/16); end
if length(opt.window) == 1, opt.window = hamming(opt.window); end
opt.overlap = length(opt.window)-1;

window = opt.window / norm(opt.window, 2); % Unit norm

Nwin = length(window);
n = size(Sx, 2);

% regenerate the full spectrum 0...2pi (minus zero Hz value)
Sx = [Sx; conj(Sx(floor((Nwin+1)/2):-1:2,:))];

% take the inverse fft over the columns
xbuf = real(ifft(Sx,[],1));

% apply the window to the columns
xbuf = xbuf .* repmat(window(:),1,size(xbuf,2));

% overlap-add the columns
x = unbuffer(xbuf,Nwin,opt.overlap);

%%% subfunction
function y = unbuffer(x,w,o)
% UNBUFFER  undo the effect of 'buffering' by overlap-add (see BUFFER)
%    A = UNBUFFER(B,WINDOWLEN,OVERLAP) returns the signal A that is
%    the unbuffered version of B.

y    = [];
skip = w - o;
N    = ceil(w/skip);
L    = (size(x,2) - 1) * skip + size(x,1);

% zero pad columns to make length nearest integer multiple of skip
if size(x,1)<skip*N, x(skip*N,end) = 0; end;

% selectively reshape columns of input into 1-d signals
for i = 1:N
    t = reshape(x(:,i:N:end),1,[]);
    l = length(t);
    y(i,l+(i-1)*skip) = 0;
    y(i,[1:l]+(i-1)*skip) = t;
end;

% overlap-add
y = sum(y,1);
y = y(1:L);

author: HeLiangHIT


text:
function test_failed = test_blocprocoffline()

% Scanario 1) Block reading from a wav and block writing to a wav
%
%
inName = 'test_in.wav';
outName = 'test_out.wav';
f = 2*rand(44100,1)*0.9-1;
wavwsave(f,44100,inName);

fs=block(inName,'offline','outfile',outName);

flag = 1;
while flag
   [fb,flag] = blockread();
   blockwrite(fb/2);
end

% Scanario 2) blockwrite from vector to wav
%
%

f = gspi;
f2 = 2*rand(numel(f),1)*0.9-1;

fs = block([f,f2],'fs',44100,'offline','outfile',outName);

flag = 1;
while flag
   [fb,flag] = blockread(44100);
   blockwrite(fb/2);
end

delete(inName);
delete(outName);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function test_ft_analysispipeline

% MEM 2gb
% WALLTIME 02:53:11
% DEPENDENCY ft_analysispipeline

% the style of this test script is also used in test_ft_datatype and test_bug2185

dirlist = {
  dccnpath('/home/common/matlab/fieldtrip/data/test/latest')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20131231')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20130630')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20121231')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20120630')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20111231')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20110630')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20101231')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20100630')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20091231')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20090630')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20081231')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20080630')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20071231')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20070630')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20061231')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20060630')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20051231')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20050630')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20040623')
  dccnpath('/home/common/matlab/fieldtrip/data/test/20031128')
  };

for j=1:length(dirlist)
  filelist = hcp_filelist(dirlist{j});

  [dummy, dummy, x] = cellfun(@fileparts, filelist, 'uniformoutput', false);
  sel = strcmp(x, '.mat');
  filelist = filelist(sel);
  clear p f x

  for i=1:length(filelist)

    % skip the large files
    d = dir(filelist{i});
    if d.bytes>50000000
      continue
    end

    try
      fprintf('processing data structure from %s\n', filelist{i});
      var = loadvar(filelist{i});
      disp(var)
    catch
      % some of the mat files are corrupt, this should not spoil the test
      disp(lasterr);
      continue
    end

    cfg = [];
    cfg.showcallinfo = 'no';
    ft_analysispipeline(cfg, var);
    set(gcf, 'Name', shortname(filelist{i}));
    drawnow
    close all

  end % for filelist

end % for dirlist
end % main function

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function str = shortname(str)
len = 50;
if length(str)>len
  begchar = length(str)-len+4;
  endchar = length(str);
  str = ['...' str(begchar:endchar)];
end
end % function shortname

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [sts] = ft_spiketriggeredspectrum_fft(cfg, data, spike)

% FT_SPIKETRIGGEREDSPECTRUM_FFT computes the Fourier spectrum (amplitude and phase)
% of the LFP around the % spikes.  A phase of zero corresponds to the spike being on
% the peak of the LFP oscillation. A phase of 180 degree corresponds to the spike being
% in the through of the oscillation. A phase of 45 degrees corresponds to the spike
% being just after the peak in the LFP.
%
% If the triggered spike leads a spike in another channel, then the angle of the Fourier
% spectrum of that other channel will be negative. Earlier phases are in clockwise
% direction.
%
% Use as
%   [sts] = ft_spiketriggeredspectrum_convol(cfg,data,spike)
% or
%   [sts] = ft_spiketriggeredspectrum_convol(cfg,data)
% where the spike data can either be contained in the DATA input or in the SPIKE input.
%
% The input DATA should be organised as the raw datatype, obtained from FT_PREPROCESSING
% or FT_APPENDSPIKE.
%
% The (optional) input SPIKE should be organised as the spike or the raw datatype,
% obtained from FT_SPIKE_MAKETRIALS or FT_PREPROCESSING (in that case, conversion is done
% within the function)
%
% Important is that data.time and spike.trialtime should be referenced relative to the
% same trial trigger times.
%
% The configuration should be according to
%   cfg.timwin       = [begin end], time around each spike (default = [-0.1 0.1])
%   cfg.foilim       = [begin end], frequency band of interest (default = [0 150])
%   cfg.taper        = 'dpss', 'hanning' or many others, see WINDOW (default = 'hanning')
%   cfg.tapsmofrq    = number, the amount of spectral smoothing through
%                      multi-tapering. Note that 4 Hz smoothing means plus-minus 4 Hz,
%                      i.e. a 8 Hz smoothing box. Note: multitapering rotates phases (no
%                      problem for consistency)
%   cfg.spikechannel = string, name of spike channels to trigger on cfg.channel = Nx1
%                      cell-array with selection of channels (default = 'all'),
%                      see FT_CHANNELSELECTION for details
%   cfg.feedback     = 'no', 'text', 'textbar', 'gui' (default = 'no')
%
% The output STS data structure can be input to FT_SPIKETRIGGEREDSPECTRUM_STAT
%
% This function uses a NaN-aware spectral estimation technique, which will default to the
% standard Matlab FFT routine if no NaNs are present. The fft_along_rows subfunction below
% demonstrates the expected function behavior.
%
% See FT_SPIKETRIGGEREDINTERPOLATION to remove segments of LFP around spikes.
% See FT_SPIKETRIGGEREDSPECTRUM_CONVOL for an alternative implementation based
% on convolution

% Copyright (C) 2008, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble provenance data spike

% check input data structure
data = ft_checkdata(data,'datatype', 'raw', 'feedback', 'yes');
if nargin==3
  spike = ft_checkdata(spike, 'datatype', {'spike'}, 'feedback', 'yes');
end

% these were supported in the past, but are not any more (for consistency with other spike functions)
cfg = ft_checkconfig(cfg, 'forbidden', {'inputfile','outputfile'});

%get the options
cfg.timwin       = ft_getopt(cfg, 'timwin',[-0.1 0.1]);
cfg.spikechannel = ft_getopt(cfg,'spikechannel', 'all');
cfg.channel      = ft_getopt(cfg,'channel', 'all');
cfg.feedback     = ft_getopt(cfg,'feedback', 'no');
cfg.tapsmofrq    = ft_getopt(cfg,'tapsmofrq', 4);
cfg.taper        = ft_getopt(cfg,'taper', 'hanning');
cfg.foilim       = ft_getopt(cfg,'foilim', [0 150]);

% ensure that the options are valid
cfg = ft_checkopt(cfg,'timwin','doublevector');
cfg = ft_checkopt(cfg,'spikechannel',{'cell', 'char', 'double', 'empty'});
cfg = ft_checkopt(cfg,'channel', {'cell', 'char', 'double'});
cfg = ft_checkopt(cfg,'feedback', 'char', {'yes', 'no'});
cfg = ft_checkopt(cfg,'taper', 'char');
cfg = ft_checkopt(cfg,'tapsmofrq', 'doublescalar');
cfg = ft_checkopt(cfg,'foilim', 'doublevector');

if strcmp(cfg.taper, 'sine')
  error('sorry, sine taper is not yet implemented');
end

% get the spikechannels
if nargin==2

  % autodetect the spikechannels and EEG channels
  [spikechannel, eegchannel] = detectspikechan(data);

  % make the final selection of spike channels and check
  if strcmp(cfg.spikechannel, 'all'),
    cfg.spikechannel = spikechannel;
  else
    cfg.spikechannel = ft_channelselection(cfg.spikechannel, data.label);
    if ~all(ismember(cfg.spikechannel,spikechannel)),
      error('some selected spike channels appear eeg channels');
    end
  end

  % make the final selection of EEG channels and check
  if strcmp(cfg.channel,'all')
    cfg.channel = eegchannel;
  else
    cfg.channel      = ft_channelselection(cfg.channel, data.label);
    if ~all(ismember(cfg.channel,eegchannel)),
      warning('some of the selected eeg channels appear spike channels');
    end
  end

  % select the data and convert to a spike structure
  tmpcfg = [];
  tmpcfg.channel = cfg.spikechannel;
  data_spk = ft_selectdata(tmpcfg, data);
  tmpcfg.channel = cfg.channel;
  data     = ft_selectdata(tmpcfg, data); % leave only LFP
  spike    = ft_checkdata(data_spk,'datatype', 'spike');
  clear data_spk % remove the continuous data
else
  cfg.spikechannel = ft_channelselection(cfg.spikechannel, spike.label);
  cfg.channel      = ft_channelselection(cfg.channel, data.label);
end

% determine the channel indices and number of chans
chansel          = match_str(data.label, cfg.channel); % selected channels
nchansel         = length(cfg.channel);                % number of channels
spikesel         = match_str(spike.label, cfg.spikechannel);
nspikesel        = length(spikesel); % number of spike channels
if nspikesel==0, error('no spike channel selected'); end

% construct the taper
if ~isfield(data, 'fsample'), data.fsample = 1/mean(diff(data.time{1})); end
begpad = round(cfg.timwin(1)*data.fsample);
endpad = round(cfg.timwin(2)*data.fsample);
numsmp = endpad - begpad + 1;
if ~strcmp(cfg.taper,'dpss')
  taper  = window(cfg.taper, numsmp);
  taper  = taper./norm(taper);
else
  % not implemented yet: keep tapers, or selecting only a subset of them.
  taper  = dpss(numsmp, cfg.tapsmofrq);
  taper  = taper(:,1:end-1);            % we get 2*NW-1 tapers
  taper  = sum(taper,2)./size(taper,2); % using the linearity of multitapering
end
taper  = sparse(diag(taper));

% preallocate the output structures for different units / trials
ntrial      = length(data.trial);
spectrum    = cell(nspikesel,ntrial);
spiketime   = cell(nspikesel,ntrial);
spiketrial  = cell(nspikesel,ntrial);

% select the frequencies
freqaxis = linspace(0, data.fsample, numsmp);
fbeg = nearest(freqaxis, cfg.foilim(1));
fend = nearest(freqaxis, cfg.foilim(2));

% update the configuration to account for rounding off differences
cfg.foilim(1) = freqaxis(fbeg);
cfg.foilim(2) = freqaxis(fend);

% make a representation of the spike, this is used for the phase rotation
spike_repr = zeros(1,numsmp);
time       = linspace(cfg.timwin(1),cfg.timwin(2), numsmp);
spike_repr(1-begpad) = 1;
spike_fft = specest_nanfft(spike_repr, time);
spike_fft = spike_fft(fbeg:fend);
spike_fft = spike_fft./abs(spike_fft);
rephase   = sparse(diag(conj(spike_fft)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute the spectra
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ft_progress('init', 'text',     'Please wait...');
for iUnit  = 1:nspikesel
  for iTrial = 1:ntrial

    % select the spikes that fell in the trial and convert to samples
    timeBins = [data.time{iTrial}  data.time{iTrial}(end)+1/data.fsample] - (0.5/data.fsample);
    hasTrial = spike.trial{spikesel(iUnit)} == iTrial; % find the spikes that are in the trial
    ts       = spike.time{spikesel(iUnit)}(hasTrial); % get the spike times for these spikes
    ts       = ts(ts>=timeBins(1) & ts<=timeBins(end)); % only select those spikes that fall in the trial window
    [ignore,spikesmp] = histc(ts,timeBins);
    if ~isempty(ts)
      ts(spikesmp==0 | spikesmp==length(timeBins)) = [];
    end

    spikesmp(spikesmp==0 | spikesmp==length(timeBins)) = [];

    % store in the output cell arrays as column vectors
    spiketime{iUnit, iTrial}  = ts(:);
    tr = iTrial*ones(size(spikesmp));
    spiketrial{iUnit, iTrial} = tr(:);

    % preallocate the spectrum
    spectrum{iUnit, iTrial} = zeros(length(spikesmp), nchansel, fend-fbeg+1);

    % compute the spiketriggered spectrum
    ft_progress(iTrial/ntrial, 'spectrally decomposing data for trial %d of %d, %d spikes for unit %d', iTrial, ntrial, length(spikesmp), iUnit);
    for j=1:length(spikesmp)

      % selected samples
      begsmp = spikesmp(j) + begpad;
      endsmp = spikesmp(j) + endpad;

      % handle spikes near the borders of the trials
      if (begsmp<1)
        segment = nan(nchansel, numsmp);
      elseif endsmp>size(data.trial{iTrial},2)
        segment = nan(nchansel, numsmp);
      else
        segment = data.trial{iTrial}(chansel,begsmp:endsmp);
      end

      % substract the DC component from every segment, to avoid any leakage of the taper
      segmentMean = repmat(nanmean(segment,2),1,numsmp); % nChan x Numsmp
      segment     = segment - segmentMean; % LFP has average of zero now (no DC)

      % taper the data segment around the spike and compute the fft
      segment_fft = specest_nanfft(segment * taper, time);

      % select the desired output frquencies and normalize
      segment_fft = segment_fft(:,fbeg:fend) ./ sqrt(numsmp/2);

      % rotate the estimated phase at each frequency to correct for the segment t=0 not being at the first sample
      segment_fft = segment_fft * rephase;

      % store the result for this spike in this trial
      spectrum{iUnit, iTrial}(j,:,:) = segment_fft;

    end % for each spike in this trial
  end % for each trial
end
ft_progress('close');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% collect the results in a structure that is a spike structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sts.lfplabel       = data.label(chansel);
sts.freq           = freqaxis(fbeg:fend);
sts.dimord         = 'rpt_chan_freq';
for iUnit = 1:nspikesel
  sts.fourierspctrm{iUnit}  = cat(1, spectrum{iUnit,:});
  spectrum(iUnit,:) = {[]}; % free from the memory
  sts.time{iUnit}           = cat(1,spiketime{iUnit,:});
  sts.trial{iUnit}          = cat(1,spiketrial{iUnit,:});
end
sts.dimord    = '{chan}_spike_lfpchan_freq';
sts.trialtime = spike.trialtime;
sts.label     = spike.label(spikesel);

% do the general cleanup and bookkeeping at the end of the function

ft_postamble previous   data
ft_postamble provenance sts
ft_postamble history    sts

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [spikelabel, eeglabel] = detectspikechan(data)

maxRate = 1000; % default on what we still consider a neuronal signal

% autodetect the spike channels
ntrial = length(data.trial);
nchans  = length(data.label);
spikechan = zeros(nchans,1);
for i=1:ntrial
  for j=1:nchans
    hasAllInts    = all(isnan(data.trial{i}(j,:)) | data.trial{i}(j,:) == round(data.trial{i}(j,:)));
    hasAllPosInts = all(isnan(data.trial{i}(j,:)) | data.trial{i}(j,:)>=0);
    fr            = nansum(data.trial{i}(j,:),2) ./ (data.time{i}(end)-data.time{i}(1));
    spikechan(j) = spikechan(j) + double(hasAllInts & hasAllPosInts & fr<=maxRate);
  end
end
spikechan = (spikechan==ntrial);

spikelabel = data.label(spikechan);
eeglabel   = data.label(~spikechan);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION for demonstration purpose
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = fft_along_rows(x)
y = fft(x, [], 2); % use normal Matlab function to compute the fft along 2nd dimension


text:
function y = fshift(x,s)
% FSHIFT Fractional circular shift
%   Syntax:
%
%       >> y = fshift(x,s)
%
%   FSHIFT circularly shifts the elements of vector x by a (possibly
%   non-integer) number of elements s. FSHIFT works by applying a linear
%   phase in the spectrum domain and is equivalent to CIRCSHIFT for integer
%   values of argument s (to machine precision).

% (c) 2005 Francois Bouffard
%          fbouffar@gel.ulaval.ca

needtr = 0; if size(x,1) == 1; x = x(:); needtr = 1; end;
N = size(x,1);
r = floor(N/2)+1; f = ((1:N)-r)/(N/2);
p = exp(-j*s*pi*f)';
y = ifft(fft(x).*ifftshift(p)); if isreal(x); y = real(y); end;
if needtr; y = y.'; end;


text:
function hpdi = hpdi(x, p)
% HPDI - Estimates the Bayesian HPD intervals
%
%   Y = HPDI(X,P) returns a Highest Posterior Density (HPD) interval
%   for each column of X. P must be a scalar. Y is a 2 row matrix
%   where ith column is HPDI for ith column of X.

%   References:
%      [1] Chen, M.-H., Shao, Q.-M., and Ibrahim, J. Q., (2000).
%          Monte Carlo Methods in Bayesian Computation. Springer-Verlag.
%
% Copyright (C) 2001 Aki Vehtari

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

if nargin < 2
  error('Not enough arguments')
end

m=size(x,2);
pts=linspace(0.1,99.9-p,20);
pt1=prctile(x,pts);
pt2=prctile(x,p+pts);
cis=abs(pt2-pt1);
[foo,hpdpi]=min(cis);
if m==1
  hpdi=[pt1(hpdpi); pt2(hpdpi)];
else
  hpdpi=sub2ind(size(pt1),hpdpi,1:m);
  hpdi=[pt1(hpdpi); pt2(hpdpi)];
end

author: gpstuff-dev

repo: gpstuff

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


text:
function [y] = local2globalVel2(V, lon,lat)

% SYNTAX:
%   [y] = local2globalVel2(V, X);
%
% INPUT:
%   V = local position vector(s)
%   lon = longitude of orifin vector in radians
%   lat = latirude of orifin vector in radians
%
% OUTPUT:
%   y = global position vector(s)
%
% DESCRIPTION:
%   Rototation from local-level reference frame to Earth-fixed reference frame

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by:
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

%initialize new position vector
y = zeros(size(V));

for i = 1 : size(V,2)
    %rotation matrix from global to local reference system
    R = [-sin(lon) cos(lon) 0;
         -sin(lat)*cos(lon) -sin(lat)*sin(lon) cos(lat);
         +cos(lat)*cos(lon) +cos(lat)*sin(lon) sin(lat)];

    %rototraslation
    y(:,i) = R\V(:,i);
end

author: goGPS-Project


text:
function R = randrot(n, N)
% Generates uniformly random rotation matrices.
%
% function R = randrot(n, N)
%
% R is a n-by-n-by-N matrix such that each slice R(:, :, i) is an
% orthogonal matrix of size n of determinant +1 (i.e., a matrix in SO(n)).
% By default, N = 1.
% Complexity: N times O(n^3).
% Theory in Diaconis and Shahshahani 1987 for the uniformity on O(n);
% With details in Mezzadri 2007,
% "How to generate random matrices from the classical compact groups."
% To ensure matrices in SO(n), we permute the two first columns when
% the determinant is -1.
%
% See also: randskew

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Sept. 25, 2012.
% Contributors:
% Change log:

    if nargin < 2
        N = 1;
    end

    if n == 1
        R = ones(1, 1, N);
        return;
    end

    R = zeros(n, n, N);

    for i = 1 : N

        % Generated as such, Q is uniformly distributed over O(n), the set
        % of orthogonal matrices.
        A = randn(n);
        [Q, RR] = qr(A);
        Q = Q * diag(sign(diag(RR))); %% Mezzadri 2007

        % If Q is in O(n) but not in SO(n), we permute the two first
        % columns of Q such that det(new Q) = -det(Q), hence the new Q will
        % be in SO(n), uniformly distributed.
        if det(Q) < 0
            Q(:, [1 2]) = Q(:, [2 1]);
        end

        R(:, :, i) = Q;

    end

end

author: andrewssobral

repo: lrslibrary


text:
function robThreeDynamicsDisplay(model, spaceNum)

% ROBTHREEDYNAMICSDISPLAY Display the robot dynamics model.

% FGPLVM

if nargin > 1
  spacing = repmat(32, 1, spaceNum);
else
  spaceNum = 0;
  spacing = [];
end
spacing = char(spacing);
fprintf(spacing);
fprintf('Third Tailored dynamics model for robot:\n')
fprintf(spacing);
fprintf('  Lambda value: %2.4f\n', model.lambda)
fprintf(spacing);
fprintf('  Sigma value: %2.4f\n', model.sigma2)

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% robust_kmeans() - an extension of Matlab kmeans() that removes outlier
%        components from all clusters.
%        This is a helper function called from pop_clust().

function  [IDX,C,sumd,D,outliers] = robust_kmeans(data,N,STD,MAXiter,method)
% data - pre-clustering data matrix.
% N - number of wanted clusters.

if nargin < 5
    method = 'kmeans';
end;

flag  = 1;
not_outliers = 1:size(data,1);
old_outliers = [];
if strcmpi(method, 'kmeans')
    [IDX,C,sumd,D] = kmeans(data,N,'replicates',30,'emptyaction','drop'); % Cluster using K-means algorithm
else
    [IDX,C,sumd,D] = kmeanscluster(data,N); % Cluster using K-means algorithm
end;
if STD >= 2 % STD for returned outlier
    rSTD = STD -1;
else
    rSTD = STD;
end
loop = 0;

while flag
     loop =  loop + 1;
	std_all = [];
    ref_D = 0;
	for k = 1:N
        tmp = ['cls' num2str(k) ' = find(IDX=='  num2str(k) ')''; ' ]; %find the component indices belonging to each cluster (cls1 = ...).
        eval(tmp);
        tmp = ['std' num2str(k) ' = std(D(cls'  num2str(k) ' ,' num2str(k) ')); ' ]; %compute the std of each cluster
        eval(tmp);
        std_all = [std_all ['std' num2str(k)  '  ']];
        tmp = [ 'ref_D = ' num2str(ref_D) ' + mean(D(cls'  num2str(k) ' ,' num2str(k) '));' ];
        eval(tmp);
	end
	std_all = [ '[ ' std_all ' ]' ];
    std_all = eval(std_all);

	% Find the outliers
    % Outlier definition - its distance from its cluster center is bigger
    % than STD times the std of the cluster, as long as the distance is bigger
    % than the mean distance times STD (avoid problems where all points turn to be outliers).
	outliers = [];
    ref_D = ref_D/N;
	for k = 1:N
        tmp = ['cls' num2str(k) '(find(D(find(IDX=='  num2str(k) ')'' , ' num2str(k) ') > ' num2str(STD)  '*std' num2str(k) ')); ' ];
        optionalO = eval(tmp);
        Oind = find(D(optionalO,k) >  ref_D*STD);
        outliers = [outliers optionalO(Oind)];
	end
    if isempty(outliers) | (loop == MAXiter)
        flag = 0;
    end
    l = length(old_outliers);
    returned_outliers = [];

    for k = 1:l
        tmp = sum((C-ones(N,1)*data(old_outliers(k),:)).^2,2)'; % Find the distance of each former outlier to the current cluster
        if isempty(find(tmp <= std_all*rSTD))  %Check if the outlier is still an outlier (far from each cluster center more than STD-1 times its std).
            returned_outliers = [returned_outliers old_outliers(k)];
        end;
    end
    outliers = not_outliers(outliers);
    outliers = [outliers returned_outliers ];
	tmp = ones(1,size(data,1));
	tmp(outliers) = 0;
	not_outliers = (find(tmp==1));

    if strcmpi(method, 'kmeans')
        [IDX,C,sumd,D] = kmeans(data(not_outliers,:),N,'replicates',30,'emptyaction','drop');
    else
        [IDX,C,sumd,D] = kmeanscluster(data(not_outliers,:),N);
    end;
    old_outliers = outliers;
    old_IDX = zeros(size(data,1),1);
    old_IDX(sort(not_outliers)) = IDX;

end

IDX = old_IDX;

author: PatternRecognition

repo: OpenBMI


text:
classdef prtPreProcEnergyNormalizeRows < prtPreProc
    % prtPreProcEnergyNormalizeRows Normalize the rows of the data to have unit
    % energy
    %

    properties (SetAccess=private)

        name = 'Energy Normalize Rows'  %  MinMax Rows
        nameAbbreviation = 'ENR'  % MMR
    end

    properties
        %no properties
        energyOffset = 0;
    end

    methods

        function self = prtPreProcEnergyNormalizeRows(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end
    end

    methods (Access = protected, Hidden = true)

        function self = trainAction(self,DataSet) %#ok<INUSD>
            %do nothing
        end

        function DataSet = runAction(self,DataSet) %#ok<MANU>

            theData = DataSet.getObservations;
            theData = bsxfun(@rdivide,theData,self.energyOffset + sqrt(sum(theData.^2,2)));
            DataSet = DataSet.setObservations(theData);
        end

        function xOut = runActionFast(self,xIn,ds) %#ok<INUSD>
            xOut = bsxfun(@rdivide,xIn,self.energyOffset + sqrt(sum(xIn.^2,2)));
        end
    end

    methods (Hidden)
        function str = exportSimpleText(self)
            titleText = sprintf('%% prtPreProcEnergyNormalizeRows\n');
            energyOffsetText = prtUtilMatrixToText(full(self.energyOffset),'varName','energyOffset');
            str = sprintf('%s%s',titleText,energyOffsetText); % No parameters
        end
    end

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function plotAcrobot(t,z,u,p)

[energy.total, energy.potential, energy.kinetic] = acrobotEnergy(z,p);

subplot(2,3,1);
plot(t,z(1,:))
xlabel('t')
ylabel('q1')
title('link one angle')
subplot(2,3,2);
plot(t,z(2,:))
xlabel('t')
ylabel('q2')
title('link two angle')
subplot(2,3,4);
plot(t,z(3,:))
xlabel('t')
ylabel('dq1')
title('link one rate')
subplot(2,3,5);
plot(t,z(4,:))
xlabel('t')
ylabel('dq2')
title('link two rate')

subplot(2,3,3); hold on
plot(t,energy.total,'k')
plot(t,energy.potential,'r')
plot(t,energy.kinetic,'b')
legend('total','potential','kinetic')
xlabel('t')
ylabel('e');
title('mechanical energy')

subplot(2,3,6)
plot(t,u)
xlabel('t')
ylabel('u')
title('torque between links')

end

author: MatthewPeterKelly

repo: OptimTraj


text:
%% Copyright (C) 2014-2016 Colin B. Macdonald
%%
%% This file is part of OctSymPy
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym isallconstant (@var{x})
%% Whether all elements of a symbolic array are constant.
%%
%% Example:
%% @example
%% @group
%% A = [1 2 sym(pi); sym(4) 5 6]
%%   @result{} A = (sym 2×3 matrix)
%%       ⎡1  2  π⎤
%%       ⎢       ⎥
%%       ⎣4  5  6⎦
%%
%% isallconstant (A)
%%   @result{} ans = 1
%% @end group
%%
%% @group
%% A(1) = sym('x')
%%   @result{} A = (sym 2×3 matrix)
%%       ⎡x  2  π⎤
%%       ⎢       ⎥
%%       ⎣4  5  6⎦
%%
%% isallconstant (A)
%%   @result{} ans = 0
%% @end group
%% @end example
%%
%% @seealso{@@sym/isconstant, @@sym/symvar, findsymbols}
%% @end defmethod

function z = isallconstant(x)

  z = isempty (findsymbols (x));

end

%!assert (isallconstant([sym(1) 2 3]))

%!test
%! syms x
%! assert (~isallconstant([sym(1) x 3]))

%!test
%! syms x
%! assert (~isallconstant([sym(1) x; sym(2) 3]))

repo: octsympy


text:
function gen_caffe_validation_imageset

% only validation images
% don't worry about padding. Use small batch size for p5 feat extraction

conf = config;

images_v = load(fullfile(conf.root_path, ...
    'validation_images_crop15_square256.mat'), 'validation_images');
images_v = images_v.validation_images;

images = cat(4, images_v.img);

validation_filenames = {images_v.filename};

assert(length(images) == length(validation_filenames));

% images must be (height x width x channels x num)
assert(ndims(images) == 4);
% assume square images
assert(size(images, 1) == size(images, 2));
% images must have 3 channels
assert(size(images, 3) == 3);
% Convert to BGR
images = images(:,:,[3 2 1],:);
% Switch width and height
images = permute(images, [2 1 3 4]);

savefilename = conf.validation_set_path;
imageset_to_leveldb(images, savefilename);
fprintf('Done writing level db\n');

save(fullfile(conf.cache_path, 'validation_filenames.mat'), 'validation_filenames');

author: rksltnl


text:
function [func, grad, Hess] = FuncGradHessSub(x, y, F, R, kin, rho)
% `SubGradHess` provides the function value, the gradient, and the Hessian
% of the subproblem of DCA and BDCA to be used in `fminunc`.
%
% USAGE:
%
%    [func, grad, Hess] = FuncGradHessSub(x, y, F, R, kin, rho)
%
% INPUTS:
%    x,y:     points
%    F:       Forward stoichiometric matrix
%    kin:     kinetics parameter in :math:`R^{2n}`
%    rho:     strongly comvex modulus
%
%
% OUTPUTS:
%    f:       function value
%    grad:    gradient
%    H:       Hessian

m            = size(F,1);
FR           = [F,R];
RF           = [R,F];
FR_plus_RF   = FR+RF;
exp_x        = exp(kin+FR'*x);
exp_y        = exp(kin+FR'*y);
ux           = FR*exp_x;
wx           = RF*exp_x;
uy_plus_wy   = FR_plus_RF*exp_y;
FRdiag_x     = FR*diag(exp_x);
dux          = FRdiag_x*FR';
dwx          = FRdiag_x*RF';
FRdiag_y     = FR*diag(exp_y);
duy_plus_dwy = FRdiag_y*FR_plus_RF';
gx           = 2*(norm(ux)^2+norm(wx)^2)+rho/2*norm(x)^2;
dgx          = 4*(dux*ux+dwx*wx)+rho*x;
dhy          = 2*duy_plus_dwy*uy_plus_wy+rho*y;

func = gx-dhy'*x;
grad = dgx-dhy;
H    = dux*dux.'+dwx*dwx.';
for j = 1:m
    H = H+ux(j)*FR*diag(exp_x.*FR(j,:)')*FR'+ ...
                                     wx(j)*FR*diag(exp_x.*RF(j,:)')*FR';
end
Hess = 4*H+rho*eye(m);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%% End of FuncGradHessSub.m %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: opencobra

repo: cobratoolbox


text:
function chronux
% This library performs time-frequency analysis (mostly using the
% multi-taper spectral estimation method) of univariate and multivariate
% data, both for continuous processes such as LFP/EEG and for point
% processes such as spike times. Point process can either be stored as
% times or as a binned process of counts. The routines in this library
% are named differently for the three cases. For calculations
% that can be performed for each of the three data types, we use suffixes
% c, pb, or pt to refer to continuous, point process binned counts, or
% point process times. For example, the spectrum calculation is performed
% mtspectrumc for continuous processes, mtspectrumpb for a binned point
% process, and mtspectrumpt for a point process consisting of times. There
% are also routines for calculating hybrid quantities involving one continuous
% and one point process. These are suffixed in a similar manner. For
% example, coherencycpb calculates the coherency between a binned point process
% and a continuous process.
%
% Certain variables are used repeatedly in this library.
%
% DATA
% data in most cases can be univariate or multivariate, and either point process,
% or continuous.
%
%      Continuous data: Continuous data is assumed to be a matrix with
%                       dimensions samples x channels/trials.
%
%      Point Process: A single time series of spike times can be in the form of
%                     a column vector.
%                     Multichannel/trial spike time data is not amenable to this
%                     storage format, since there are generally different
%                     number of spikes in each channel/trial. Instead,
%                     multichannel/trial spike data is stored in a structure
%                     array. A structure is a matlab data object with various
%                     fields. These fields contain the elements
%                       e.g. The command data=struct('times',[]); creates an empty
%                            structure with field 'times'. Similarly, the command
%                            data=struct('times',[1 2 3]); creates the structure with
%                            the field 'times' containing integers 1, 2, and 3.
%
%                     We can also have a structure array (or an array of structures)
%                     defined for example, by
%                     data(1)=struct('times',rand(1,100)); and
%                     data(2)=struct('times',rand(1,200));
%                     This is a 2 dimensional structure array where the
%                     first field is a 100 dimensional random vector, and
%                     the second field is a 200 dimensional random vector.
%                     This format allows storage of multichannel point
%                     process times in a single variable data.
%
%                     The above holds for point processes stored as times.
%                     If instead, the point processes are binned, then one
%                     can use a matrix to represent them
%
%
%      Summary: data - array of continuous data with dimensions time x channels
%                      structural array of spike times with dimensions
%                               equal to the number of channels
%                      1d array of spike times as a column vector
%                      array of binned spike counts with dimensions time x channels
%
% PARAMETERS:
% These are various parameters used in the spectral calculations. Since
% these parameters are used by most routines in Chronux, they are stored in
% a single structure params. The fields of params are
%
% tapers : precalculated tapers from dpss or in the one of the following
%          forms:
%          (1) A numeric vector [TW K] where TW is the
%              time-bandwidth product and K is the number of
%              tapers to be used (less than or equal to
%              2TW-1).
%          (2) A numeric vector [W T p] where W is the
%              bandwidth, T is the duration of the data and p
%              is an integer such that 2TW-p tapers are used. In
%              this form there is no default i.e. to specify
%              the bandwidth, you have to specify T and p as
%              well. Note that the units of W and T have to be
%              consistent: if W is in Hz, T must be in seconds
%              and vice versa. Note that these units must also
%              be consistent with the units of params.Fs: W can
%              be in Hz if and only if params.Fs is in Hz.
%              The default is to use form 1 with TW=3 and K=5
%
%
% pad:   (padding factor for the FFT) - optional (can take values -1,0,1,2...).
%         -1 corresponds to no padding, 0 corresponds to padding
%         to the next highest power of 2 etc.
%			  e.g. For N = 500, if PAD = -1, we do not pad; if PAD = 0, we pad the FFT
%			       to 512 points, if pad=1, we pad to 1024 points etc.
%			       Defaults to 0.
%
% Fs:sampling frequency.optional (default 1)
%
%
% fpass: frequencies in an fft calculation can range from 0 to Fs/2 where
%        Fs is the sampling frequency. Sometimes it may be useful to
%        compute fourier transforms (and resulting quantities like the
%        spectrum over a smaller range of frequencies). This is specified
%        by fpass, which can be in the form [fmin fmax] where fmin >=0 and
%        fmax<=Fs/2. optional (default [0 Fs/2])
%
% err=[errtype p] calculates theoretical error bars (confidence levels)
%                 when errtype=1 and jackknife error bars when errchk=2. In each case, the
%                 error is calculated at a p value specified by p. -
%                 optional (default 0)
%
% trialave: trialave controls whether or not to average over channels/trials for
%           multichannel/trial analyses. trialave=0 (default) implies no trial
%           averaging, trialave=1 implies that the quantity of interest is averaged
%           over channels/trials. optional (default 0)
%
% Other parameters are discussed in individual routines as and when they
% are used.

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%% Copyright (C) 2014-2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym erf (@var{x})
%% Symbolic erf function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = erf (x)
%%   @result{} y = (sym) erf(x)
%% @end group
%% @end example
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @end defmethod

function y = erf(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('erf', x);
end

%!error erf (sym(1), 2)
%!assert (isequaln (erf (sym(nan)), sym(nan)))

%!shared x, d
%! d = 1;
%! x = sym('1');

%!test
%! f1 = erf(x);
%! f2 = erf(d);
%! assert( abs(double(f1) - f2) < 1e-15 )

%!test
%! D = [d d; d d];
%! A = [x x; x x];
%! f1 = erf(A);
%! f2 = erf(D);
%! assert( all(all( abs(double(f1) - f2) < 1e-15 )))

%!test
%! % round trip
%! y = sym('y');
%! A = erf (d);
%! f = erf (y);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B, -eps)


text:
function element_size = p08_element_size ( )

%*****************************************************************************80
%
%% P08_ELEMENT_SIZE returns a typical element size for problem 08.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, real ( kind = 8 ) ELEMENT_SIZE, a typical element size.
%
  element_size = 0.005;

  return
end

author: johannesgerer


text:
% mri_exp_mult_mex_test
% test mri_exp_mult_mex

L = 10;
N = 30;
M = 20;
rng(0)
A = randn(N,L) + 1i * randn(N,L);
ur = randn(N,1);
ui = randn(N,1);
vr = randn(M,1);
vi = randn(M,1);
u = ur + 1i * ui;
v = vr + 1i * vi;

d1 = A' * exp(-ur * v.');
d2 = mri_exp_mult_mex(A, ur, v);
printf('double u real case error: %g%%', max_percent_diff(d1, d2))

d1 = A' * exp(-u * vr.');
d2 = mri_exp_mult_mex(A, u, vr);
printf('double v real case error: %g%%', max_percent_diff(d1, d2))

if 0 % test both complex, both real
	mri_exp_mult_mex(A, ur, vr);
	mri_exp_mult_mex(A, u, v);
end

if ~is_pre_v7
	u = single(u);
	v = single(v);
	ur = single(ur);
	vr = single(vr);
	A = single(A);

	d1 = A' * exp(-ur * v.');
	d2 = mri_exp_mult_mex(A, ur, v);
	printf('single u real case error: %g%%', max_percent_diff(d1, d2))

	d1 = A' * exp(-u * vr.');
	d2 = mri_exp_mult_mex(A, u, vr);
	printf('single v real case error: %g%%', max_percent_diff(d1, d2))
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
%> @brief struct_overlay Overlay default fields with input fields
%> @author Eric Cousineau <eacousineau@gmail.com>
%> @note From optOverlay() on Mathworks, modified
function [opts] = struct_overlay(opts_default, opts_in, options)

    % Simple cases
    if iscell(opts_default)
        opts_default = struct(opts_default{:});
    end
    if isempty(opts_in)
        opts = opts_default;
        return;
    end
    if iscell(opts_in)
        opts_in = struct(opts_in{:});
    end

	is_valid = @(o) isstruct(o) && length(o) == 1;
	assert(is_valid(opts_default), 'opts_default must be scalar struct or structable cell array');
	assert(is_valid(opts_in), 'opts_in must be scalar struct or structable cell array');

    persistent func_opts_default;
    if isempty(func_opts_default)
        func_opts_default = struct('Recursive', true, 'AllowNew', false);
    end
    if nargin < 3
        options = func_opts_default;
    else
        if iscell(options)
            options = struct(options{:});
        end
        options = struct_overlay(func_opts_default, options);
    end

    %TODO Tackle struct arrays?
    new_fields = fieldnames(opts_in);
    options_index = find(strcmp(new_fields, 'overlay_options'));
    indices = 1:length(new_fields);
    if ~isempty(options_index)
        options = struct_overlay(options, opts_in.overlay_options);
        indices(options_index) = [];
    end

    if ~options.AllowNew
        % Check to see if there are bad fields
        fields = fieldnames(opts_default);
        unoriginal = setdiff(new_fields(indices), fields);
        if ~isempty(unoriginal)
            error(['Fields are not in original: ', implode(unoriginal, ', ')]);
        end
    end

    opts = opts_default;

    for i = 1:length(indices)
        index = indices(i);
        field = new_fields{index};
        new_value = opts_in.(field);
        % Apply recursion
        if options.Recursive && isfield(opts, field)
            cur_value = opts.(field);
            % Both values must be proper option structs
            if is_valid(cur_value) && is_valid(new_value)
                new_value = struct_overlay(cur_value, new_value, options);
            end
        end
        opts.(field) = new_value;
    end
end

author: ayonga

repo: frost-dev

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


text:
function pcut = addSinCosCuts(p)

pcut = p;
if ~isempty(p.evalMap)
    sin_ = [];
    cos_ = [];
    for i = 1:length(p.evalMap)
        if isequal(p.evalMap{i}.fcn,'sin')
            sin_ = [sin_;p.evalMap{i}.variableIndex p.evalMap{i}.computes];
        elseif isequal(p.evalMap{i}.fcn,'cos')
            cos_ = [cos_;p.evalMap{i}.variableIndex p.evalMap{i}.computes];
        end
    end
    if ~isempty(sin_) && ~isempty(cos_) && ~isempty(intersect(sin_(:,1),cos_(:,1)))
        for i = 1:size(sin_,1)
            j = find(sin_(i,1)==cos_(:,1));
            if ~isempty(j)
                k1 = sin_(i,2);
                k2 = cos_(j,2);
                k3 = sin_(i,1);
                % x_k1 and x_k2 rpresent sin and cos of x_k3

                % Construct the SOCP model for 1 >= x_k1^2 + x_k2^2
                % i.e. norm([x_k1;x_k2]) <= 1
                % i.e. [1;x_k1;x_k2] in socp cone
                if 0
                    F_structemp = spalloc(3,length(p.c)+1,5);
                    F_structemp(:,1) = [1;0;0];
                    F_structemp(2,1+k1) = 1;
                    F_structemp(3,1+k2) = 1;
                    K.f = 0;
                    K.l = 0;
                    K.s = 0;
                    K.e = 0;
                    K.q = 3;
                    localModel1 = createNumericalModel(F_structemp,K);
                    pcut = mergeNumericalModels(pcut,localModel1);
                end

                % We also add simple linear cuts, in case the lower bound
                % solver doesn't support socps. FIX ME: Decide globally,
                % but make sure it is consistent with stand-alone envelope
                % y <= sqrt(2)-x, y>=-sqrt(2)+x,
                % y <= sqrt(2)+x, y>=-sqrt(2)-x
                % sqrt(2)-x-y, sqrt(2)+x-y, sqrt(2)-x+y, sqrt(2)+x+y
                F_structemp = spalloc(4,length(p.c)+1,5);
                F_structemp(:,1) = sqrt(2)*[1;1;1;1];
                F_structemp(1,1+k1) = -1;
                F_structemp(1,1+k2) = -1;
                F_structemp(2,1+k1) = 1;
                F_structemp(2,1+k2) = -1;
                F_structemp(3,1+k1) = -1;
                F_structemp(3,1+k2) = 1;
                F_structemp(4,1+k1) = 1;
                F_structemp(4,1+k2) = 1;
                K.f = 0;
                K.l = 4;
                K.s = 0;
                K.e = 0;
                K.q = 0;
                localModel2 = createNumericalModel(F_structemp,K);
                pcut = mergeNumericalModels(pcut,localModel2);
            end
        end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%dcm (Cba) to rotation vector (r^a)
%Cba: B to A transformation matrix
%r^a=Rotation vector defined in A. (Rotate A around r with mag(r) to obtain B)
%the magnitude of the returning rotation vector will be between [0 and PI]

%This is the implementaion of savage(3-31).
function rvec=dcm2rot(dcm)

sinPHI=0.5*sqrt((dcm(3,2)-dcm(2,3))^2+(dcm(1,3)-dcm(3,1))^2+(dcm(2,1)-dcm(1,2))^2);
cosPHI=0.5*(dcm(1,1)+dcm(2,2)+dcm(3,3)-1);

PHI=atan2(sinPHI,cosPHI);

if (cosPHI>=0)
    if (sinPHI==0)
        F=1;
    else
        F=PHI/sinPHI;
    end
    rvec=0.5*F*[dcm(3,2)-dcm(2,3);dcm(1,3)-dcm(3,1);dcm(2,1)-dcm(1,2)];
else
    sr_a=1-cosPHI;
    mu1=sqrt((dcm(1,1)-1)/sr_a+1);
    mu2=sqrt((dcm(2,2)-1)/sr_a+1);
    mu3=sqrt((dcm(3,3)-1)/sr_a+1);
    if (mu1>=mu2 && mu1>=mu3)
       u1=mu1*sign(dcm(3,2)-dcm(2,3));
       u2=1/2/u1*(dcm(1,2)+dcm(2,1))/sr_a;
       u3=1/2/u1*(dcm(1,3)+dcm(3,1))/sr_a;
    elseif (mu2>=mu3)
       u2=mu2*sign(dcm(1,3)-dcm(3,1));
       u3=1/2/u2*(dcm(2,3)+dcm(3,2))/sr_a;
       u1=1/2/u2*(dcm(1,2)+dcm(2,1))/sr_a;
    else
       u3=mu3*sign(dcm(2,1)-dcm(1,2));
       u1=1/2/u3*(dcm(1,3)+dcm(3,1))/sr_a;
       u2=1/2/u3*(dcm(2,3)+dcm(3,2))/sr_a;
    end
    rvec=PHI*[u1;u2;u3];
end

repo: ekfmonoslam


text:
function [rf rootIndex] = buildForest( rf_model )
%(treemap, nodestatus, nodeclass, bestvar, xbestsplit, ndbigtree)
    % extract all the fields
    treemap = double(rf_model.treemap);
    nodestatus = double(rf_model.nodestatus);
    nodeclass = double(rf_model.nodeclass);
    bestvar = double(rf_model.bestvar);
    xbestsplit = double(rf_model.xbestsplit);
    ndbigtree = double(rf_model.ndbigtree);

    [nrnodes ntree] = size(nodestatus);
    nnodes = sum(ndbigtree(1:ntree));
    rf = cell(1, nnodes);
    rootIndex = zeros(1, ntree);
    maskIndex = 0 * nodestatus;

    % first create nodes
    index = 1;
    for n = 1 : ntree
        rootIndex(n) = index;
        for jx = 1 : nrnodes
            if nodestatus(jx, n) ~= 0
                rf{index} = RFTreeNode( nodestatus(jx,n) == -1, nodeclass(jx, n),...
                    bestvar(jx, n), xbestsplit(jx, n) );
                maskIndex(jx, n) = index;
                index = index + 1;
            end
        end
    end

    assert( nnodes == index - 1 );

    % then create classification trees by linking the nodes
    for n = 1 : ntree
        for jx = 1 : nrnodes
            if nodestatus(jx, n) ~= 0
                nodeIndex = maskIndex(jx, n);
                treemat = treemap(:, 2*n-1 : 2*n);
                treemat = treemat(:);

                if treemat(2*jx-1) > 0
                    leftChildIndex = maskIndex( treemat(2*jx-1), n );
                    % rf{nodeIndex}.insertLeftChild( rf{leftChildIndex} );
                    rf{nodeIndex}.leftChild = leftChildIndex;
                end

                if treemat(2*jx) > 0
                    rightChildIndex = maskIndex( treemat(2*jx), n );
                    % rf{nodeIndex}.insertRightChild( rf{rightChildIndex} );
                    rf{nodeIndex}.rightChild = rightChildIndex;
                end
            end
        end
    end
end

function buildTree(nodes, r, treemat)

end

%     r = RFTreeNode;
%     r.bestVar = 'r';
%     n1 = RFTreeNode;
%     n1.bestVar = 'n1';
%     n2 = RFTreeNode;
%     n2.bestVar = 'n2';
%     n3 = RFTreeNode;
%     n3.bestVar = 'n3';
%     n3.isTerminal = true;
%     n4 = RFTreeNode;
%     n4.bestVar = 'n4';
%     n5 = RFTreeNode;
%     n5.bestVar = 'n5';
%     n5.isTerminal = true;
%     n6 = RFTreeNode;
%     n6.bestVar = 'n6';
%     n6.isTerminal = true;
%     n7 = RFTreeNode;
%     n7.bestVar = 'n7';
%     n7.isTerminal = true;
%     n8 = RFTreeNode;
%     n8.bestVar = 'n8';
%     n8.isTerminal = true;
%
%     r.insertLeftChild( n1 );
%     r.insertRightChild( n2 );
%     n1.insertLeftChild( n3 );
%     n1.insertRightChild( n4 );
%     n4.insertLeftChild( n7 );
%     n4.insertRightChild( n8 );
%     n2.insertLeftChild( n5 );
%     n2.insertRightChild( n6 );
%
% %     rf = cell(1, 8);
% %     rf{1} = n1;
% %     rf{2} = n2;
% %     rf{3} = n3;
% %     rf{4} = n4;
% %     rf{5} = n5;
% %     rf{6} = n6;
% %     rf{7} = n7;
% %     rf{8} = n8;
% %     rf{9} = r;
%
%     % rf = zeros(1, 8);
%     rf(1) = r;
%     rf(2) = n1;
%     rf(3) = n2;
%     rf(4) = n3;
%     rf(5) = n4;
%     rf(6) = n5;
%     rf(7) = n6;
%     rf(8) = n7;
%     rf(9) = n8;

author: zhangqianqianQQ


text:
function D = prtDistanceEuclidean(dataSet1,dataSet2)
%prtDistanceEuclidean   Euclidean distance
%
%   DIST = prtDistanceCityBlock(DS1,DS2) calculates the Euclidean distance
%   from all the observations in datasets DS1 to DS2, and ouputs a distance
%   matrix of size DS1.nObservations x DS2.nObservations. DS1 and DS2
%   should have the same number of features. DS1 and DS2 should be
%   prtDataSet objects.
%
%   For more information, see:
%
%   http://en.wikipedia.org/wiki/Euclidean_distance
%
% Example:
%
%   % Create 2 data sets
%   dsx = prtDataSetStandard('Observations', [0 0; 1 1]);
%   dsy = prtDataSetStandard('Observations', [1 0;2 2; 3 3]);
%   % Compute distance
%   distance = prtDistanceEuclidean(dsx,dsy)
%
% See also: prtDistanceCityBlock, prtDistanceChebychev
% prtDistanceMahalanobis, prtDistanceSquare, prtDistanceLnorm

[data1,data2] = prtUtilDistanceParseInputs(dataSet1,dataSet2);
D = prtDistanceLNorm(data1,data2,2);


text:
function amd_install
%AMD_INSTALL compile and install amd2 for use in MATLAB
%   Your current directory must be AMD/MATLAB for this function to work.
%
% Example:
%   amd_install
%
% See also amd, amd2.

% Copyright 1994-2007, Tim Davis, Patrick R. Amestoy, and Iain S. Duff.

% This orders the same matrix as the ANSI C demo, amd_demo.c.  It includes an

amd_make
addpath (pwd)
fprintf ('\nThe following path has been added.  You may wish to add it\n') ;
fprintf ('permanently, using the MATLAB pathtool command.\n') ;
fprintf ('%s\n\n', pwd) ;
amd_demo

repo: gspbox