text:
% script to specify network arch and train drdae on single machine
% make copies of this script for each architecture you want to try

%% setup paths for code. assumed this script runs in its own directory
codeDir = '/afs/cs.stanford.edu/u/awni/scr/noise_proj/audio_ufl/matlab_wd/drdae';
minFuncDir = '/afs/cs/u/amaas/scratch/matlab_trunk/parallel_proto/minFunc/';
baseDir = '/afs/cs.stanford.edu/u/awni/scr/noise_proj/aurora2/features/';
%% AMAAS setup
%codeDir = '/afs/cs.stanford.edu/u/amaas/scratch/audio/audio_repo/matlab_wd/drdae';
%% add paths
addpath(codeDir);
addpath(minFuncDir);

%% setup network architecture
eI = [];
% dimension of each input frame
eI.featDim = 14;
eI.dropout = 0;
% context window size of the input.
eI.winSize = 3;
% weight tying in hidden layers
% if you want tied weights, must have odd number of *hidden* layers
eI.tieWeights = 0;
% 2 hidden layers and output layer
eI.layerSizes = [512 eI.featDim];
% highest hidden layer is temporal
eI.temporalLayer = 0;
% dim of network input at each timestep (final size after window & whiten)
eI.inputDim = eI.featDim * eI.winSize;
% length of input sequence chunks.
eI.seqLen = [1 10 25 50 100];
% activation function
eI.activationFn = 'tanh';
% temporal initialization type
eI.temporalInit = 'rand';
% weight norm penaly
eI.lambda = 0;
% file containing whitening matrices for outputs
eI.targetWhiten = [codeDir '/aurora_whiten.mat'];
%% setup weight caching
saveDir = '/scail/group/deeplearning/speech/awni/aurora_results/';
eI.saveDir = [saveDir 'model_1hl512_ws3_small_test/'];
mkdir(eI.saveDir);
%% initialize weights
[stack_i, W_t_i] = initialize_weights(eI);
[theta] = rnn_stack2params(stack_i, eI, W_t_i);

[stack_new, W_t_new] = rnn_params2stack(theta,eI);
[theta_new] = rnn_stack2params(stack_new, eI, W_t_new);
%% Directory of features
eI.featInBase =baseDir; % '/afs/cs.stanford.edu/u/amaas/scratch/aurora2/features/';

%% load data
eI.useCache = 0;

%Get SNR List
snrList = {'clean1', 'N1_SNR5', 'N1_SNR10', 'N1_SNR15', 'N1_SNR20', ...
'clean2', 'N2_SNR5', 'N2_SNR10', 'N2_SNR15', 'N2_SNR20', ...
'clean3', 'N3_SNR5', 'N3_SNR10', 'N3_SNR15', 'N3_SNR20', };
eI.subdirs = snrList;
[data_cell, targets_cell] = load_aurora( baseDir, 'Mfc08_multiTR', snrList, -1, eI );

% dieplay mean as a whitening debug check
%disp(mean(data_cell{1},2));
%% setup minFunc
options.Diagnostics = 'on';
options.Display = 'iter';
options.MaxIter = 2000;
options.MaxFunEvals = 2500;
options.Corr = 50;
options.DerivativeCheck = 'off';
options.outputFcn = @save_callback;
%% run optimizer
minFunc(@drdae_obj, theta, options, eI, data_cell, targets_cell, false, false);

author: amaas

repo: rnn-speech-denoising


text:
function d = get_dtw_distance(d1, d2)
    SM = simmx(abs(d1), abs(d2));
%     [~, ~, D] = dpfast(1-SM);
    [~, ~, D] = dp(1-SM);
    d = D(size(D,1), size(D,2));
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function test_suite = test_TestSuiteInDir
%test_TestSuiteInDir Unit test for TestSuiteInDir class.

%   Steven L. Eddins
%   Copyright 2009 The MathWorks, Inc.

initTestSuite;

function test_constructor
this_test_path = fileparts(which(mfilename));
cwd_test_dir = fullfile(this_test_path, 'cwd_test');
suite = TestSuiteInDir(cwd_test_dir);

assertEqual(suite.Name, 'cwd_test');
assertEqual(suite.Location, cwd_test_dir);

function test_gatherTestCases
this_test_path = fileparts(which(mfilename));
cwd_test_dir = fullfile(this_test_path, 'cwd_test');
suite = TestSuiteInDir(cwd_test_dir);
suite.gatherTestCases();

assertEqual(numel(suite.TestComponents), 3);

author: ricedsp


text:
function W_down = get_W_down(W)
N = size(W, 2);
W_down = zeros(2, 2 * N^2);
for u2 = 0 : 1
    for y1 = 1 : N
        for y2 = 1 : N
            for u1 = 0 : 1
%                 P = 0.5 * W(mod(u1 + u2, 2) + 1, y1) * W(u2 + 1, y2);
%                 if P == 0
%                     W_down(u2 + 1, 2 * N * (y1 - 1) + 2 * (y2 - 1) + u1 + 1) = realmin;
%                 else
%                     W_down(u2 + 1, 2 * N * (y1 - 1) + 2 * (y2 - 1) + u1 + 1) = P;
%                 end
                W_down(u2 + 1, 2 * N * (y1 - 1) + 2 * (y2 - 1) + u1 + 1) = 0.5 * W(mod(u1 + u2, 2) + 1, y1) * W(u2 + 1, y2);
            end
        end
    end
end
end

author: YuYongRun

repo: PolarCodeDecodersInMatlab


text:
close all;
clear all;
clc;
rng('default');
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');
% Signal space
N = 256;
% Number of measurements
M = 32;
% Number of signals
Ss = [1, 2, 4, 8, 16, 32];
% Sparsity levels
Ks = 1:25;
num_trials = 100;
num_ks = length(Ks);
num_ss = length(Ss);
bp_success_with_k = zeros(num_ss, num_ks);

snr_threshold = 100;

for ns=1:num_ss
    % Current sparsity level
    S = Ss(ns);
    for nk=1:num_ks
        K = Ks(nk);
        num_bp_successes = 0;
        for nt=1:num_trials
            % Sensing matrix
            Phi = spx.dict.simple.gaussian_dict(M, N);
            X = model_1_data(N, K, S);
            % Measurement vectors
            Y = Phi * X;

            % OMP MMV solver instance (using P = 2 for P-SOMP)
            omp_solver = spx.pursuit.joint.OrthogonalMatchingPursuit(Phi, K, 2);
            % Solve the sparse recovery problem
            result = omp_solver.solve(Y);
            % Solution vectors
            X_OMP = result.Z;
            % Comparison
            cs = spx.commons.SparseSignalsComparison(X, X_OMP, K);
            snr = cs.cum_signal_to_noise_ratio;
            bp_success = snr > snr_threshold;
            num_bp_successes = num_bp_successes + bp_success;
            fprintf('S: %d, K=%d, trial=%d, residual omp: %e, SNR: %f dB\n'...
                , S, K, nt, cs.cum_difference_norm, snr);
        end
        bp_success_with_k(ns, nk) = num_bp_successes / num_trials;
    end
end

save ('bin/figure_1_spherical_dict_model_1_somp_success_with_k.mat');

author: indigits

repo: sparse-plex


text:
function Cn = correlation_image(Y,sz,d1,d2,flag_norm, K)

% construct correlation image based on neighboing pixels
% Y: raw data
% sz: define the relative location of neighbours. it can be scalar, 2
%       element vector or a binary matrix
%       scalar: number of nearest neighbours, either 4 or 8, default 4
%       2 element vector: [rmin, rmax], the range of neighbours. the
%           distance between the neighbour and the pixel is d, dmin <= r <
%           dmax.
%       matrix: a squared matrix (2k+1)*(2k+1)indicating the location of neighbours
% d1,d2: spatial dimensions
% flag_norm: indicate whether Y has been normalized and centered ( 1 is
%   yes, 0 is no)
% K:  scalar, the rank of the random matrix for projection

% Author: Eftychios A. Pnevmatikakis, Simons Foundation, 2015
% with modifications from Pengcheng Zhou, Carnegie Mellon University, 2015.
% It uses convolution and random projection for speeding up the
% computation.

%% preprocess the raw data
if ~exist('flag_norm', 'var') || isempty(flag_norm)
    flag_norm = false;
end

if ~exist('sz', 'var') || isempty(sz)
    sz = [0,1,0; 1,0,1; 0,1,0];
end

% center data
Y = bsxfun(@minus, double(Y), mean(Y, ndims(Y)));
if ~ismatrix(Y)
    [d1, d2, T] = size(Y);
else
    T = size(Y, 2);
end

if exist('K', 'var') && (~isempty(K))
    Y = double(reshape(Y, [], T))*randn(T, K);
    % centering
    mY = mean(Y,2);
    Y = bsxfun(@minus, Y, mY);        % normalizing
    flag_norm = false;
end
if ~flag_norm
    sY = sqrt(mean(Y.*Y, ndims(Y)));
    sY(sY==0) = 1; % avoid nan values
    Y = bsxfun(@times, Y, 1./sY);
end
if ismatrix(Y)
    Y = reshape(Y, d1, d2, []);
end

%% construct a matrix indicating location of the matrix
if  isscalar(sz)
    if sz == 8      % 8 nearest neighbours
        sz = [1,1,1; 1,0,1; 1,1,1];
    elseif sz==4
        sz = [0,1,0; 1,0,1; 0,1,0];
    end
elseif length(sz(:)) == 2
    % the specified neighbours has a distance within the domain [dmin,
    % dmax)
    sz = ceil(sz);
    dmin = min(sz); dmax = max(sz);
    rsub = (-dmax+1):(dmax-1);      % row subscript
    csub = rsub;      % column subscript
    [cind, rind] = meshgrid(csub, rsub);
    R = sqrt(cind.^2+rind.^2);
    sz = (R>=dmin) .* (R<dmax);
end

%% compute the correlation
Yconv = imfilter(Y, sz);        % sum over the neighbouring pixels
MASK = imfilter(ones(d1,d2), sz);   % count the number of neighbouring pixels
Cn = mean(Yconv.*Y, 3)./MASK;   % compute correlation and normalize

author: zhoupc


text:
%  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 2 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, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function klaudert
%KLAUDERT Unit test for function KLAUDER.

%	O. Lemoine - March 1996.

N=256;

% Admissibility condition
psi=klauder(N);
if abs(sum(psi))>sqrt(eps),
 error('klauder test 1 failed');
end

% Admissibility condition
psi=klauder(N,100,0.04);
if abs(sum(psi))>sqrt(eps),
 error('klauder test 2 failed');
end

% Admissibility condition
psi=klauder(N,1,0.49);
if abs(sum(psi))>sqrt(eps),
 error('klauder test 3 failed');
end

N=227;

% Admissibility condition
psi=klauder(N);
if abs(sum(psi))>sqrt(eps),
 error('klauder test 4 failed');
end

% Admissibility condition
psi=klauder(N,100,0.04);
if abs(sum(psi))>sqrt(eps),
 error('klauder test 5 failed');
end

% Admissibility condition
psi=klauder(N,1,0.49);
if abs(sum(psi))>sqrt(eps),
 error('klauder test 6 failed');
end

author: HeLiangHIT


text:
function t=getDCTTransform(im,N)
s=dctmtx(N);
t=s*im*s';

author: Sable

repo: mcbench-benchmarks


text:
function spm_dcm_local_minima(DCM)
% evaluates the free energy landscape around the posterior
% FORMAT: spm_dcm_local_minima(DCM)
% DCM - (invert) model structure
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dcm_local_minima.m 5892 2014-02-23 11:00:16Z karl $

% find dimension of greatest curvature
%==========================================================================
U  = spm_svd(DCM.Cp);
U  = U(:,end);
Cu = U'*DCM.Cp*U;
Su = spm_sqrtm(Cu);

% Free energy landscape
%==========================================================================
DCM.options.DATA = 0;
DCM.M.Nmax = 1;
DCM.name   = 'test';

N = 128;
s = linspace(-8,8,N);
for i = 1:N

    Ep           = spm_vec(DCM.Ep);
    Ep           = Ep + s(i)*U*Su*(U'*Ep);
    DCM.M.P      = spm_unvec(Ep,DCM.Ep);
    try
        [Qp,Cp,Eh,F] = spm_nlsi_GN(DCM.M,DCM.xU,DCM.xY);
        FF(i)        = F;
    catch
        FF(i)    = NaN;
    end
end

plot(s,FF)
title('Free energy','FontSize',16);
xlabel('parameter 1')
ylabel('parameter 2')
axis square

author: spm


text:
function plot_probabilities(sv,specs)

if nargin<2

    specs=[3,3];

end

myLimits=[-sqrt(eps),1+sqrt(eps)];

r0=specs(1);

c0=specs(2);

[f,the_regimes]=load_filters(sv);

cellform = iscell(f);

if cellform

    all_regimes=fieldnames(f{1}.smoothed_regime_probabilities);

    all_states=fieldnames(f{1}.smoothed_state_probabilities);

    regs=the_regimes{1};

    for ii=2:numel(f)

        [all_regimes,IA,IB]=intersect(all_regimes,...
            fieldnames(f{ii}.smoothed_regime_probabilities)); %#ok<ASGLU>

        regs=regs(IA);

        [all_states]=intersect(all_states,...
            fieldnames(f{ii}.smoothed_state_probabilities));

    end

    the_regimes=regs;

else

    all_regimes=fieldnames(f.smoothed_regime_probabilities);

    all_states=fieldnames(f.smoothed_state_probabilities);

end

utils.plot.multiple(@plotfuncr,all_regimes,...
    'smoothed regime probabilities',r0,c0);

utils.plot.multiple(@plotfuncs,all_states,...
    'smoothed state probabilities',r0,c0);

    function [tex,leg]=plotfuncr(vname)

        d=load_item('smoothed_regime_probabilities',vname);

        plot(d,'linewidth',2)

        ylim(myLimits)

        tex=vname;

        vname(1:numel('regime'))=[];

        vname=strrep(vname,'_','');

        tex=sprintf('%s(%s)',tex,the_regimes{str2double(vname)});

        leg='';

    end

    function [tex,leg]=plotfuncs(vname)

        d=load_item('smoothed_state_probabilities',vname);

        plot(d,'linewidth',2)

        ylim(myLimits)

        tex=vname;

        leg='';

    end

    function d=load_item(a,b)

        if cellform

            d=f{1}.(a).(b);

            for jj=2:numel(f)

                d=[d,f{jj}.(a).(b)]; %#ok<AGROW>

            end

        else

            d=f.(a).(b);

        end

    end

end

author: jmaih


text:
function V = base64decode(RES)
  % V = base64decode(RES)
  %
  % Inputs:
  %   RES  #RES long list of characters
  % Outputs:
  %   V  #V long list of bytes
  %
  % Faster implementation of:
  %  Documentation for matlab.net.base64decode
  %
  % See also: base64encode
  vals = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/=';
  equ = find(RES=='=')-1;
  if isempty(equ)
    equ = numel(RES);
  end
  [~,inds] = ismember(RES(1:equ),vals);
  bin = dec2bin_logical(inds-1,6);
  %assert(isequal(bin,dec2bin(inds-1,6)=='1'));

  bin = reshape(bin',1,[]);
  bin = bin(1:floor(numel(bin)/8)*8);
  bin = reshape(bin,8,[])';
  % This is very slow...
  %V = bin2dec(bin);
  V = (bin)*(2.^(7:-1:0))';
end

author: alecjacobson

repo: gptoolbox


text:
addpath(genpath('.'))
load('./Metadata/SUNRGBDMeta.mat')
% load('./Metadata/SUNRGBD2Dseg.mat')
%% change path
imageId = 1;
for i = 1:size(SUNRGBDMeta, 2)
    SUNRGBDMeta(i).depthpath = strrep(SUNRGBDMeta(i).depthpath, '/n/fs/sun3d/data/', '/home/siyuan/Documents/Dataset/SUNRGBD_ALL/');
    SUNRGBDMeta(i).rgbpath = strrep(SUNRGBDMeta(i).rgbpath, '/n/fs/sun3d/data/', '/home/siyuan/Documents/Dataset/SUNRGBD_ALL/');
end
data = SUNRGBDMeta(imageId);
[rgb,points3d,depthInpaint,imsize]=read3dPoints(data);

%% draw
figure,
imshow(data.rgbpath);
hold on;
for kk =1:length(data.groundtruth3DBB)
    rectangle('Position', [data.groundtruth3DBB(kk).gtBb2D(1) data.groundtruth3DBB(kk).gtBb2D(2) data.groundtruth3DBB(kk).gtBb2D(3) data.groundtruth3DBB(kk).gtBb2D(4)],'edgecolor','y');
    text(data.groundtruth3DBB(kk).gtBb2D(1),data.groundtruth3DBB(kk).gtBb2D(2),data.groundtruth3DBB(kk).classname,'BackgroundColor','y')
end
%% draw 3D
figure,
vis_point_cloud(points3d,rgb)
hold on;
for kk =1:length(data.groundtruth3DBB)
   vis_cube(data.groundtruth3DBB(kk),'r')
end

%%
% anno2d = SUNRGBD2Dseg(imageId);
% figure,
% imagesc(anno2d.seglabel);
% % category name in 37 categories list
% load('./Metadata/seg37list.mat');
% objectname37 = unique(anno2d.seglabel(:));
% objectname37 = seg37list(objectname37(objectname37~=0));
%
% figure,
% imagesc(anno2d.seglabelall);
% % category name of all categories
% objectnameall = anno2d.names

%% Example to read single data
% data = readframeSUNRGBD('/n/fs/sun3d/data/SUNRGBD/kv2/kinect2data/000002_2014-05-26_14-23-37_260595134347_rgbf000103-resize/','/n/fs/sun3d/data/');
% [rgb,points3d,depthInpaint,imsize]=read3dPoints(data);
% %% draw
% figure,
% imshow(data.rgbpath);
% hold on;
% for kk =1:length(data.groundtruth3DBB)
%     rectangle('Position', [data.groundtruth3DBB(kk).gtBb2D(1) data.groundtruth3DBB(kk).gtBb2D(2) data.groundtruth3DBB(kk).gtBb2D(3) data.groundtruth3DBB(kk).gtBb2D(4)],'edgecolor','y');
%     text(data.groundtruth3DBB(kk).gtBb2D(1),data.groundtruth3DBB(kk).gtBb2D(2),data.groundtruth3DBB(kk).classname,'BackgroundColor','y')
% end
% %% draw 3D
% figure,
% vis_point_cloud(points3d,rgb)
% hold on;
% for kk =1:length(data.groundtruth3DBB)
%    vis_cube(data.groundtruth3DBB(kk),'r')
% end

author: thusiyuan


text:
function normVal=dualQuatNorm(dq)
%%DUALQUATNORM Compute the norm of a dual quaternion. This is defined as
%           dq*dq^*, where dq^* is the conjugate of the dual quaternion.
%           The result is a dual quaternion with no hypercomplex components
%           (but it might have both a real and a dual component). Dual
%           quaternions are often used for simultaneously representing
%           orientation and position. A dual quaternion consists of two
%           parts dq=q1+eps*q2, where q1 and q2 are quaternions and eps is
%           a dual number. Dual numbers are such that eps^2=0 and have
%           various rules for multiplication with complex numbers. Note
%           that the norm does not depend on the handedness of the
%           quaternion algebra, even though the dualQuatMult and quatMult
%           functions do.
%
%INPUTS: dq A dual quaternion represented as a 4X2 matrix dq(:,1) is the
%           non-dual quaternion and dq(:,2) is the dual quaternion (the
%           ones times eps, the dual number). The elements of each
%           quaternion are ordered q(1,1)+i*q(2,1)+j*q(3,1)+k*q(4,1), where
%           i,j and k are the typical hypercomplex numbers.
%
%OUTPUTS: normVal The value of the norm of the dual quaternion. This is a
%                 dual quaternion with no complex parts, but it might have
%                 a dual part (normVal(1,:) can be nonzero). A unit dual
%                 quaternion has only a real part and no dual part
%                 (normVal(1,1)=1,normVal(2,1)=0) nor complex parts.
%
%Dual quaternions and common operations including computation of the norm
%of the dual quaternion are discussed in [1]. The norm of a standard
%(non-dual) quaternion is generally just taken to be the same as the norm
%of a typical 4X1 vector.
%
%REFERENCES:
%[1] B. Kenwright, "A beginners guide to dual-quaternions: What they are,
%    how they work, and how to use them for 3D character hierarchies," in
%    Proceedings of the 20th International Conference on Computer Graphics,
%    Visualization and Computer Vision, Prague, Czech Republic, 24-27 Jun.
%    2012.
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

normVal=dualQuatMult(dq,dualQuatConj(dq));

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [ y, x ] = cvx_check_dimlist( x, emptyok )

% CVX_CHECK_DIMLIST Verifies the input is a valid dimension list.

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

if nargin < 2 || emptyok,
    xmin = 0;
else
    xmin = 1;
end
if isa( x, 'cell' ),
    nel = numel( x );
    xnew = zeros( 1, nel );
    y = false;
    for k = 1 : nel,
        if ~isnumeric( x{k} ) || length( x{k} ) ~= 1, return; end
        xnew( k ) = x{k};
    end
    x = xnew;
end
if isnumeric( x ) && ndims( x ) <= 2 && ~all( size( x ) > 1 ) && isreal( x ) && ~any( isnan( x ) | isinf( x ) ) && all( x >= xmin ) && all( floor( x ) == x ),
    y = true;
else
    y = false;
end
if y && nargout > 1,
    x = [ x( : )', 1, 1 ];
    x = x( 1 : max( [ 2, find( x > 1 ) ] ) );
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function sp = spdimcc(n, d)
% SPDIMCC   Compute number of sparse grid points, Clenshaw-Curtis grid
%    SP = SPDIMCC(N,D) Computes the number of sparse grid points
%	   of the Clenshaw-Gurtis grid of dimension D and level N.
%    (Internal function)

% The function uses the formulas for N <= 7 given in A. Schreiber:
% Die Methode von Smoliak bei der multivariaten Interpolation,
% (2000), p.35/36.

% Author : Andreas Klimke, Universitaet Stuttgart
% Version: 1.0
% Date   : August 5, 2003

% ------------------------------------------------------------
% Sparse Grid Interpolation Toolbox
% Copyright (c) 2006 W. Andreas Klimke, Universitaet Stuttgart
% Copyright (c) 2007-2008 W. A. Klimke. All Rights Reserved.
% See LICENSE.txt for license.
% email: klimkeas@ians.uni-stuttgart.de
% web  : http://www.ians.uni-stuttgart.de/spinterp
% ------------------------------------------------------------

switch n
 case 0
	sp = 1;
 case 1
	sp = 2*d + 1;
 case 2
	sp = 2*d^2 + 2*d + 1;
 case 3
	sp = round((4*d^3 + 6*d^2 + 14*d)/3) + 1;
 case 4
	sp = round((2*d^4 + 4*d^3 + 22*d^2 + 20*d)/3) + 1;
 case 5
	sp = round((4*d^5 + 10*d^4 + 100*d^3 + 170*d^2 + 196*d)/15) + 1;
 case 6
	sp = round((4*d^6 + 12*d^5 + 190*d^4 + 480*d^3 + 1246*d^2 + 948* ...
							d)/45) + 1;
 case 7
	sp = round((8*d^7 + 28*d^6 + 644*d^5 + 2170*d^4 + 9632*d^3 + ...
							15442*d^2 + 12396*d)/315) + 1;
 otherwise
	sp = round((8*d^7 + 28*d^6 + 644*d^5 + 2170*d^4 + 9632*d^3 + ...
							15442*d^2 + 12396*d)/315) + 1;
	for m = 8:n
		seq = spgetseq(m,d);
		for k = 1:size(seq,1)
			temp = 1;
			for l = 1:d
				lval = seq(k,l);
				if lval == 0
				elseif lval < 3
					temp = temp * 2;
				else
					temp = temp * 2^double(lval-1);
				end
			end
			sp = sp + temp;
		end
	end
end

author: johannesgerer

repo: jburkardt-m


text:
function [m,varargout] = max(varargin)

funname = 'max';

narginchk(1, 2)
if nargin == 1;
    nargoutchk(0, 2)
    [m,varargout{1}] = max(varargin{1}.Z(:));

elseif nargin == 2;
    nargoutchk(0, 1)

    % which one of the two input arguments is a GRIDobj
    isGRIDobj = cellfun(@(x) isa(x,'GRIDobj'),varargin);

    % if both
    if all(isGRIDobj)
        validatealignment(varargin{1},varargin{2});
        M = builtin(funname,varargin{1}.Z,varargin{2}.Z);
    else
        varargin = varargin(~isGRIDobj + 1);
        if ~isscalar(varargin{2});
            validatealignment(varargin{1},varargin{2});
        end
        M = builtin(funname,varargin{1}.Z,varargin{2});
    end

    m = varargin{1};
    m.name = funname;
    m.Z = M;
    m.zunit = '';
end

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
%{
 * Copyright (C) 2020-2030, The Regents of The University of Michigan.
 * All rights reserved.
 * This software was developed in the Biped Lab (https://www.biped.solutions/)
 * under the direction of Jessy Grizzle, grizzle@umich.edu. This software may
 * be available under alternative licensing terms; contact the address above.
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 1. Redistributions of source code must retain the above copyright notice, this
 *    list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 *    this list of conditions and the following disclaimer in the documentation
 *    and/or other materials provided with the distribution.
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * The views and conclusions contained in the software and documentation are those
 * of the authors and should not be interpreted as representing official policies,
 * either expressed or implied, of the Regents of The University of Michigan.
 *
 * AUTHOR: Bruce JK Huang (bjhuang[at]umich.edu)
 * WEBSITE: https://www.brucerobot.com/
%}

function cost = checkCost(point, cons1, cons2)
    if point >= cons1 && point <= cons2
        cost = 0;
    else
        dis = min(abs(point-cons2), abs(point-cons1));
        cost = getCost(dis);
    end
end

author: UMich-BipedLab


text:
function M = CGDA(X, Y, d, lamda)

% X:  training samples
% Y:  training labels
% d:  reduced dimension d * N

[~, N] = size(X);

%% compute the block-wise graph (W) based on the labels

W = zeros(N, N);
num_labels = length(unique(Y));
L = [];

for i = 1 : num_labels
    ind = (find(Y == i));
    l = length(ind);
    L = [L, l];
    x = X(:, ind);
    w = (x' * x + lamda * eye(size(x' * x))) \ (x' * x);
    if i == 1
        W(1 : l, 1 : l) = w;
    else
        W(sum(L(1 : i-1)) +1 : sum(L(1 : i)), sum(L(1 : i-1)) +1 : sum(L(1 : i))) = w;
    end
end

W = W - diag(diag(W));
W = max(W, W');

%% embedding
M = CGDA_LPP(X, d, W); % M:  learned mapping

end

author: BehnoodRasti


text:
function [f] = betad(z)
%BETAD Dirichlet Beta function
%
%usage: f = betad(z)
%
%tested on version 5.3.1
%
%      This program calculates the Dirichlet Beta function
%      for the elements of Z using the Dirichlet deta function.
%      Z may be complex and any size.
%
%      Note: this is NOT the beta function defined by
%            Gamma(x)*Gamma(y)/Gamma(x+y)
%
%      Has zeros for z=(-odd integers),
%      and infinite number of zeros for z=1/2+i*y
%
%
%see also: Zeta, Deta, Eta, Lambda, Bern, Euler

%Paul Godfrey
%pgodfrey@conexant.com
%3-24-01

f=deta(z,2);

return


text:
function [F] = TriScatteredInterpVector(V,D)
  % TRISCATTEREDINTERPVECTOR simple extension of TriScatteredInterp to handle
  % vector values at each data point
  %
  % [F] = TriScatteredInterpVector(V,D)
  %
  % Inputs:
  %   V  #V by 2 list of data positions
  %   D  #V by m list of data vectors
  % Outputs:
  %   A  m-size cell of TriScatteredInterp classes corrsepnding to each
  %     coordinate of data vectors
  %   ev  function handle to evaluate data at given points (X,Y)
  %

  m = size(D,2);

  % make room in cell array
  A = cell(1,m);
  % loop over data coordinates
  for ii = 1:m
    % build interpolant for this coordinate
    interp_func = TriScatteredInterp(V,D(:,ii));
    A{ii} = interp_func;
  end

  % function handle to eval interps
  F = @EvalTriScatteredInterpVector;

  function DXY = EvalTriScatteredInterpVector(X,Y)
    % for convenience build the interpolantion evaluation function handle
    f = @(I) I(X,Y);
    DXY = cellfun(f,A, 'UniformOutput', false);
    DXY = reshape(cell2mat(DXY),[size(X) m]);
  end
end


text:
  pjob  = createParallelJob ( 'configuration', 'lyrics_ithaca_2009b' );
%
% Set the number of workers
%
  set(pjob, 'MinimumNumberOfWorkers', 3);
  set(pjob, 'MaximumNumberOfWorkers', 3);
  createTask ( pjob, @systolicLyricsFilter, 1, {'mjam.wav'} )

  submit(pjob);
  waitForState(pjob, 'finished');

  disp('Parallel Job is finished')
  results = getAllOutputArguments(pjob);
%
%  Clean up the job
%
  destroy ( pjob );

  [ tune, fs ] = wavread ( 'mjam.wav' );
  disp('Play the original version')
  sound(tune, fs)

  disp('Play the filtered version')
  sound(results{3}, fs)


text:
% Aligns the point cloud given a rotation matrix.
%
% Args:
%   points3d - Nx3 point cloud.
%   R - 3x3 rotation matrix.
%
% Returns:
%   points3d - Nx3 aligned point cloud.
%
% Author: Nathan Silberman (silberman@cs.nyu.edu)
function points3d = get_aligned_point_cloud(points3d, R)
  points3d = (R * points3d')';
end

author: zhirongw


text:
classdef Hopkins155 < handle

% The variables inside the ground-truth file are organized as follow:
% . width and height: dimensions (in pixels) of all the frames in the
%   video sequence.
% . points: number of tracked points P.
% . frames: number of frames F.
% . y: a matrix 3xPxF containing the homogeneous coordinates of the P
%   points in the F frames.
% . x: a matrix 3xPxF derived from y by normalizing the first two
%   components of each vector such that they belong to the interval
%   [-1;1].
% . K: the 3x3 normalization matrix used to pass from y to x
%   (x=K^(-1)*x).
% . s: a Px1 vector containing the ground-truth segmentation; for each
%   point it gives the index of the corresponding motion group.

properties(SetAccess=private)
    root_dir
    % names of examples
    example_names
    % paths of example data
    example_paths
    % data of each example
    examples
    % number of motions
    num_motions
end

methods

    function self = Hopkins155()
        self.root_dir = spx.data.local.hopkins155_dir();
        self.example_names = cell(0, 1);
        self.example_paths = cell(0, 1);
        self.load();
    end

    function load(self)
        listing = dir(self.root_dir);
        for i=1:length(listing)
            % we iterate through each directory
            entry = listing(i);
            if ~entry.isdir
                continue;
            end
            % ignore . and .. directories
            if strcmp(entry.name, '.')
                continue;
            end
            if strcmp(entry.name, '..')
                continue;
            end
            subdir = fullfile(self.root_dir, entry.name);
            sub_listing = dir(subdir);
            % the directory must have a file with _truth.mat extension
            valid_data = false;
            for j=1:length(sub_listing)
                subentry = sub_listing(j);
                fname = subentry.name;
                tmp = strfind(fname, '_truth.mat');
                if ~isempty(tmp)
                    valid_data = true;
                    break;
                end
            end
            if ~valid_data
                % we can ignore this directory
                continue;
            end
            self.example_names{end+1} = entry.name;
            self.example_paths{end+1} = fullfile(subdir, fname);
        end
        % data will be loaded in this array
        self.examples = cell(1, length(self.example_names));
        self.num_motions = zeros(1, length(self.example_names));
    end

    function result = num_examples(self)
        result = length(self.example_names);
    end

    function result = get_example(self, n)
        if n < 1 || n > self.num_examples()
            error('Invalid example number.');
        end
        result = self.examples{n};
        if isempty(result)
            % this example is not yet loaded.
            path = self.example_paths{n};
            data = load(path);
            % let's fill up more info
            s = data.s;
            % ensure that s is a row vector
            if iscolumn(s)
                s = s';
            end
            result.num_motions = max(s);
            % let us work through the labels to reorder them
            % such that points from each subspace are adjacent.
            mapping = [];
            counts = [];
            labels = [];
            for k=1:result.num_motions
                % entries for k-th cluster
                entries = find(s == k);
                mapping = [mapping entries];
                count = length(entries);
                counts = [counts count];
                labels = [labels k*ones(1, count)];
            end
            % labels of points
            result.labels = labels;
            % counts of each cluster
            result.counts = counts;
            x = data.x;
            % reorder x
            x = x(:, mapping, :);
            % number of points
            S = size(x, 2);
            % number of frames
            F = size(x, 3);
            result.num_frames = F;
            result.num_points = S;
            % Ambient dimension
            M = 2 * F;
            result.M  = M;
            self.num_motions(n) = result.num_motions;
            % reshape the data in x for clustering purpose
            % drop the third component of each point in x [homogeneous part]
            x = x(1:2,:,:);
            % now x is 2xPxF (P : points, F is frames)
            % bring frames dimension ahead and points dimension behind
            x = permute(x,[1 3 2]);
            % Now map the 2 coordinates of each point in each column
            X = reshape(x,M,S);
            result.X = X;
            % associate example name
            result.name = self.example_names{n};
            % clean up unnecessary fields
            %fields = {'x', 'y', 's', 'width', 'height', 'frames', 'K'};
            %rmfield(result, fields);
            % finally store it.
            self.examples{n} = result;
        end
    end

    function load_all_examples(self)
        n = self.num_examples;
        for i=1:n
            self.get_example(i);
        end
    end

    function result = get_all_examples(self)
        self.load_all_examples();
        result = self.examples;
    end

    function result = get_2_3_motions(self)
        self.load_all_examples();
        index2 = find(self.num_motions == 2);
        index3 = find(self.num_motions == 3);
        index = [index2 index3];
        result = self.examples(index);
    end

    function result = num_2_motions(self)
        result = sum(self.num_motions == 2);
    end

    function result = num_3_motions(self)
        result = sum(self.num_motions == 3);
    end

    function result = num_5_motions(self)
        result = sum(self.num_motions == 5);
    end

    function result = get_2_motions(self)
        self.load_all_examples();
        % identify 2 motions
        indices = find(self.num_motions == 2);
        result = self.examples(indices);
    end

    function result = get_3_motions(self)
        self.load_all_examples();
        % identify 3 motions
        indices = find(self.num_motions == 3);
        result = self.examples(indices);
    end

    function result = get_5_motions(self)
        self.load_all_examples();
        % identify 5 motions
        indices = find(self.num_motions == 5);
        result = self.examples(indices);
    end

    function describe(self)
        % ensure that all examples are loaded
        self.load_all_examples();
        fprintf('2 motions: %d\n', self.num_2_motions());
        fprintf('3 motions: %d\n', self.num_3_motions());
        fprintf('5 motions: %d\n', self.num_5_motions());
    end

end

end


text:
function dfupdatexlim(newminmax,updateplots)
%DFUPDATEXLIM Update the stored x axis min/max values

%   $Revision: 1.1.6.5 $  $Date: 2004/01/24 09:36:03 $
%   Copyright 2003-2004 The MathWorks, Inc.

minmax = [];                     % to become new x limits
oldminmax = dfgetset('xminmax'); % previous limits
ftype = dfgetset('ftype');
if nargin==0
   newminmax = [];
end

if isempty(newminmax) && isequal(ftype, 'icdf')
   % Default limits span most of the probability range
   minmax = [.01 .99];
elseif isempty(newminmax)
   % Update limits from datasets with a plotting flag on
   dsdb = getdsdb;
   ds = down(dsdb);
   while(~isempty(ds))
      if ds.plot == 1
         minmax = combineminmax(minmax,ds.xlim);
      end
      ds = right(ds);
   end

   % Update from fits with a plotting flag on
   fitdb = getfitdb;
   ft = down(fitdb);
   while(~isempty(ft))
      if ft.plot == 1
         minmax = combineminmax(minmax,xlim(ft));
      end
      ft = right(ft);
   end
else
   minmax = newminmax;
end

% Now update plot
dffig = dfgetset('dffig');
if ~isempty(minmax) && isequal(zoom(dffig,'getmode'),'off')
   ax = get(dffig,'CurrentAxes');
   islinscale = isequal(get(ax,'XScale'),'linear');
   if ~islinscale && any(minmax<=0)
      warning('stats:dfupdatexlim:NegativeDataIgnored',...
              'Negative data ignored.');
      minmax = [1e-6 1] * max(abs(minmax));
   end
   if isempty(newminmax) && ~isequal(ftype, 'icdf')
      % Adjust axis limits to include a margin around plotted points
      if islinscale
         dx = diff(minmax) * 0.01 * [-1 1];
         if all(dx==0), dx = [-1 1]; end
      else
         dlogx = .01 * diff(log(minmax));
         if dlogx==0, dlogx = 1; end
         dx = [minmax(1) * exp(-dlogx),  minmax(2) * exp(dlogx)] - minmax;
      end
   elseif minmax(1)==minmax(2)
      if islinscale
         dx = [-1 1];
      else
         dx = [minmax(1)/2, 2*minmax(1)];
      end
   else
      % Don't adjust the limits that were passed in or computed
      dx = 0;
   end
   oldxlim = get(ax,'XLim');
   newxlim = minmax + dx;
   if ~isequal(oldxlim,newxlim)
      set(ax,'XLim',newxlim);
      if nargin<2 || updateplots
         dfupdateallplots(false,true);
      end
   end
end
dfgetset('xminmax',minmax);

% ------------ Helper to combine old and new minmax values
function bothmm = combineminmax(oldmm,newmm)

if isempty(oldmm)
   bothmm = newmm;
elseif isempty(newmm)
   bothmm = oldmm;
else
   bothmm = [min(oldmm(1),newmm(1)) max(oldmm(2),newmm(2))];
end

author: zouchuhang

repo: LayoutNet


text:
% Script demonstrating usage of the cbpdndlms function.
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2017-04-29
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'Copyright' and 'License' files
% distributed with the library.

% Training images
S0 = zeros(512, 512, 2, 'single');
S0(:,:,1) = single(stdimage('lena.grey')) / 255;
S0(:,:,2) = single(stdimage('barbara.grey')) / 255;

%Reduce images size to speed up demo script
tmp = zeros(128, 128, 2, 'single');
for k = 1:size(S0,3),
  tmp(:,:,k) = imresize(S0(:,:,k), 0.25);
end
S0 = tmp;

% Filter input images and compute highpass images
npd = 16;
fltlmbd = 5;
[Sl, Sh] = lowpass(S0, fltlmbd, npd);

% Construct weight matrix and padded test image set
Shp = padarray(Sh, [7 7], 'post');
t = 0.5;
W = randn(size(Sh));
W(abs(W) > t) = 1;
W(abs(W) < t) = 0;
W = padarray(W, [7 7], 'post');
ShW = W .* Shp;

% Construct initial dictionary
D0 = zeros(8,8,32, 'single');
D0(3:6,3:6,:) = single(randn(4,4,32));

% Set up cbpdndl parameters
lambda = 0.05;
opt = [];
opt.Verbose = 1;
opt.MaxMainIter = 500;
opt.rho = 50*lambda + 0.5;
opt.sigma = size(Sh,3);
opt.AutoRho = 1;
opt.AutoRhoPeriod = 10;
opt.AutoSigma = 1;
opt.AutoSigmaPeriod = 10;
opt.XRelaxParam = 1.8;
opt.DRelaxParam = 1.8;

% Do standard dictionary learning and reconstruct
[D1, X1, optinf1] = cbpdndl(D0, ShW, lambda, opt);
DX1 = ifft2(bsxfun(@times, fft2(D1, size(X1,1), size(X1,2)), fft2(X1)), ...
           'symmetric');
Sr1 = squeeze(sum(DX1,3)) + padarray(Sl, [7 7], 'post');

% Do dictionary learning with additive mask simulation and reconstruct
opt.W = W;
[D2, X2, optinf2] = cbpdndlms(D0, ShW, lambda, opt);
DX2 = ifft2(bsxfun(@times, fft2(D2, size(X2,1), size(X2,2)), fft2(X2)), ...
           'symmetric');
Sr2 = squeeze(sum(DX2,3)) + padarray(Sl, [7 7], 'post');

% Display dictionaries
figure;
subplot(1,2,1);
imdisp(tiledict(D1));
title('Standard DL');
subplot(1,2,2);
imdisp(tiledict(D2));
title('DL with additive mask simulation');

% Display reconstructions
figure;
subplot(2,2,1);
imdisp(Sr1(:,:,1));
title('Standard DL');
subplot(2,2,2);
imdisp(Sr2(:,:,1));
title('DL with additive mask simulation');
subplot(2,2,3);
imdisp(Sr1(:,:,2));
title('Standard DL');
subplot(2,2,4);
imdisp(Sr2(:,:,2));
title('DL with additive mask simulation');

% Plot functional value evolution
figure;
subplot(1,2,1);
semilogx(optinf1.itstat(:,2), 'LineWidth', 2);
ylim([15, 45]);
xlabel('Iterations');
ylabel('Functional value');
title('Standard DL');
subplot(1,2,2);
semilogx(optinf2.itstat(:,2), 'LineWidth', 2);
ylim([15, 45]);
xlabel('Iterations');
ylabel('Functional value');
title('DL with additive mask simulation');

author: thfylsty


text:
function [] = run_me_first(varargin)

    if nargin < 1
        message_flag = true;
    else
        message_flag = varargin{1};
    end

    % Add folders to path.
    addpath(pwd);

    cd solver/;
    addpath(genpath(pwd));
    cd ..;

    cd auxiliary/;
    addpath(genpath(pwd));
    cd ..;

    cd plotter/;
    addpath(genpath(pwd));
    cd ..;

    cd data/;
    addpath(genpath(pwd));
    cd ..;

    cd applications/;
    addpath(genpath(pwd));
    cd ..;

    if message_flag
        nmflibrary_message();
    end

end

author: hiroyuki-kasai

repo: NMFLibrary


text: %SUM Dataset overload

author: marianux

repo: ecg-kit

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


text:
classdef TestPointPolygonTest
    %TestPointPolygonTest

    methods (Static)
        function test_contour_matrix
            cntr = [0 0; 1 0; 2 2; 3 3; 3 4];
            b = cv.pointPolygonTest(cntr, [2.3 2.4]);
            validateattributes(b, {'double'}, {'scalar'});

            % purely-integer points (both contour and query)
            b = cv.pointPolygonTest(int32(cntr), [2 2]);
            validateattributes(b, {'double'}, {'scalar'});
        end

        function test_contour_cellarray
            cntr = {[0 0], [1 0], [2 2], [3 3], [3 4]};
            b = cv.pointPolygonTest(cntr, [2.3 2.4]);
            validateattributes(b, {'double'}, {'scalar'});
        end

        function test_measure_dist
            cntr = [0 0; 1 0; 2 2; 3 3; 3 4];

            b = cv.pointPolygonTest(cntr, [2.3,2.4], 'MeasureDist',true);
            validateattributes(b, {'double'}, {'scalar', 'real'});

            b = cv.pointPolygonTest(cntr, [2.3,2.4], 'MeasureDist',false);
            validateattributes(b, {'double'}, {'scalar', 'integer'});
            assert(b == 0 || b == 1 || b == -1);
        end

        function test_error_argnum
            try
                cv.pointPolygonTest();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function as_m_i(A)
%
%  Generates of the data used in 'as_m'. Data are stored in global
%  variables.
%
%  Calling sequence:
%
%    as_m_i(A)
%
%  Input:
%
%    A         real, symmetric matrix.
%
%
%  LYAPACK 1.0 (Thilo Penzl, May 1999)

if nargin~=1
  error('Wrong number of input arguments.');
end

if norm(A-A','fro')~=0
  error('A is not symmetric!');
end

global LP_A

LP_A = A;


text:
function [inside] = ft_inside_vol(pos, vol)

% FT_INSIDE_VOL locates dipole locations inside/outside the source
% compartment of a volume conductor model.
%
% [inside] = ft_inside_vol(pos, vol, ...)
%
% where the input should be
%   pos      Nx3 matrix with dipole positions
%   vol      structure with volume conductor model
% and the output is
%   inside   list of dipoles inside the brain compartment
%            (1=inside, 0=outisde)
%
% Additional optional input arguments should be given in key value pairs
% and can include
%   <none>

% Copyright (C) 2003-2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% 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: ft_inside_vol.m 2885 2011-02-16 09:41:58Z roboos $

% determine the type of volume conduction model
switch ft_voltype(vol)

  % single-sphere or multiple concentric spheres
  case {'singlesphere' 'concentric'}
    if ~isfield(vol, 'source')
      % locate the innermost compartment and remember it
      [dum, vol.source] = min(vol.r);
    end
    if isfield(vol, 'o')
      % shift dipole positions toward origin of sphere
      tmp = pos - repmat(vol.o, size(pos,1), 1);
    else
      tmp = pos;
    end
    distance = sqrt(sum(tmp.^2, 2))-vol.r(vol.source);
    % positive if outside, negative if inside
    inside   = distance<0;

    % multi-sphere volume conductor model
  case 'multisphere'

    % nspheres = size(vol.r,1);
    % ndipoles = size(pos,1);
    % inside = zeros(ndipoles,1);
    % for sph=1:nspheres
    % for dip=1:ndipoles
    %   if inside(dip)
    %     % the dipole has already been detected in one of the other spheres
    %     continue
    %   end
    %   inside(dip) = (norm(pos(dip,:) - vol.o(sph,:)) <= vol.r(sph));
    % end
    % end
    % outside = find(inside==0);
    % inside  = find(inside==1);

    % this is a much faster implementation
    nspheres = size(vol.r,1);
    ndipoles = size(pos,1);
    inside = zeros(ndipoles,1);
    for sph=1:nspheres
      % temporary shift dipole positions toward origin
      if isfield(vol, 'o')
        tmp = pos - repmat(vol.o(sph,:), [ndipoles 1]);
      else
        tmp = pos;
      end
      flag = (sqrt(sum(tmp.^2,2)) <= vol.r(sph));
      inside = inside + flag;
    end
    inside  = inside>0;

    % realistic BEM volume conductor model
  case {'bem', 'dipoli', 'bemcp', 'asa', 'avo', 'nolte', 'neuromag'}
    if ~isfield(vol, 'source')
      % locate the innermost compartment and remember it
      vol.source = find_innermost_boundary(vol.bnd);
    end
    % use the specified source compartment
    pnt = vol.bnd(vol.source).pnt;
    tri = vol.bnd(vol.source).tri;
    % determine the dipole positions that are inside the brain compartment
    inside  = bounding_mesh(pos, pnt, tri);

    % unrecognized volume conductor model
  otherwise
    error('unrecognized volume conductor model');
end

% ensure that these are column vectors
inside(find(isnan(inside(:)))) = 0;
inside  = logical(inside(:));

author: PatternRecognition

repo: OpenBMI


text:
function [ n_data, tc, sigma ] = surten_values ( n_data )

%*****************************************************************************80
%
%% SURTEN_VALUES returns some values of the surface tension.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Lester Haar, John Gallagher and George Kell,
%    NBS/NRC Steam Tables:
%    Thermodynamic and Transport Properties and Computer Programs
%    for Vapor and Liquid States of Water in SI Units,
%    Hemisphere Publishing Corporation, Washington, 1984,
%    TJ270.H3, pages 267.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real TC, the temperature, in degrees Celsius.
%
%    Output, real SIGMA, the surface tension,
%    in Pascal * m = Newton / m.
%
  n_max = 14;

  sigma_vec = [ ...
     74.22E+00, ...
     72.74E+00, ...
     71.20E+00, ...
     69.60E+00, ...
     67.95E+00, ...
     58.92E+00, ...
     48.75E+00, ...
     37.68E+00, ...
     26.05E+00, ...
     14.37E+00, ...
      8.78E+00, ...
      3.67E+00, ...
      0.40E+00, ...
      0.00E+00 ];

  tc_vec = [ ...
      10.000E+00, ...
      20.000E+00, ...
      30.000E+00, ...
      40.000E+00, ...
      50.000E+00, ...
     100.000E+00, ...
     150.000E+00, ...
     200.000E+00, ...
     250.000E+00, ...
     300.000E+00, ...
     325.000E+00, ...
     350.000E+00, ...
     370.000E+00, ...
     373.976E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    tc = 0.0;
    sigma = 0.0;
  else
    tc = tc_vec(n_data);
    sigma = sigma_vec(n_data);
  end

  return
end


text:
function b = r8mat_gesl ( a, n, pivot, b, job )

%*****************************************************************************80
%
%% R8MAT_GESL solves a system factored by R8MAT_GEFA.
%
%  Discussion:
%
%    This is a simplified version of the LINPACK routine DGESL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real A(N,N), the LU factors from R8MAT_GEFA.
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer PIVOT(N), the pivot vector from R8MAT_GEFA.
%
%    Input, real B(N), the right hand side vector.
%
%    Input, integer JOB, specifies the operation.
%    0, solve A * x = b.
%    nonzero, solve A' * x = b.
%
%    Output, real B(N), the solution vector.
%

%
%  Solve A * x = b.
%
  if ( job == 0 )
%
%  Solve PL * Y = B.
%
    for k = 1 : n - 1

      l = pivot(k)

      if ( l ~= k )
        temp = b(l);
        b(l) = b(k);
        b(k) = temp;
      end

      b(k+1:n) = b(k+1:n) + a(k+1:n,k) * b(k);

    end
%
%  Solve U * X = Y.
%
    for k = n : -1 : 1
      b(k) = b(k) / a(k,k);
      b(1:k-1) = b(1:k-1) - a(1:k-1,k) * b(k);
    end
%
%  Solve A' * X = B.
%
  else
%
%  Solve U' * Y = B.
%
    for k = 1 : n
      b(k) = ( b(k) - b(1:k-1) * a(1:k-1,k) ) / a(k,k);
    end
%
%  Solve ( PL )' * X = Y.
%
    for k = n - 1 : -1 : 1

      b(k) = b(k) + b(k+1:n) * a(k+1:n,k);

      l = pivot(k);

      if ( l ~= k )
        temp = b(l);
        b(l) = b(k);
        b(k) = temp;
      end

    end

  end

  return
end


text:
%% (Internal) Wave definitions for ECG delineators used in the toolbox
%
%   Fields used in the payload variables generated by ECGtask_ECG_delineation
%
%
% See also ECGtask_ECG_delineation
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 09/05/2016
% Birthdate  : 09/05/2016
% Copyright 2008-2016
%

% ecg-kit struct field format definition
cAnnotationOutputFields = { 'Pon' 'P' 'Poff' 'Ptipo' 'QRSon' 'qrs' 'Q' 'R' 'S' 'QRSoff' 'Ton' 'T' 'Tprima' 'Toff' 'Ttipo' };
cAnnotationSLRfields =    { 'Pon' 'P' 'Poff'         'QRSon' 'qrs'             'QRSoff' 'Ton' 'T' 'Tprima' 'Toff' };

% conversion from HL7-aECG to ecg-kit format
cHL7aECG_translation = ...
    { 'MDC_ECG_WAVC_PWAVE',   'P', 'Pon', 'Poff'; ...
      'MDC_ECG_WAVC_QWAVE',   'Q', '', ''; ...
      'MDC_ECG_WAVC_RWAVE',   'R', '', ''; ...
      'MDC_ECG_WAVC_SWAVE',   'S', '', ''; ...
      'MDC_ECG_WAVC_QRSWAVE', 'qrs', 'QRSon', 'QRSoff'; ...
      'MDC_ECG_WAVC_QSWAVE',  'qrs', 'QRSon', 'QRSoff'; ...
      'MDC_ECG_WAVC_TWAVE',   'T', 'Ton', 'Toff' };

% conversion from ecg-kit format to Physionet's ecgcodes.h (wfdb_ECGCODES_H)
cMIT_wave_translation = { ...
                    'Pon'    39; ...
                    'P'      24; ...
                    'Poff'   40; ...
                    'Ptipo'  nan; ...
                    'QRSon'  39; ...
                    'qrs'    1; ...
                    'Q'      nan; ...
                    'R'      nan; ...
                    'S'      nan; ...
                    'QRSoff' 40; ...
                    'Ton'    39; ...
                    'T'      27; ...
                    'Tprima' 27; ...
                    'Toff'   40; ...
                    'Ttipo'  nan; ...
                    };


text:
function b = r8po_inverse ( n, r )

%*****************************************************************************80
%
%% R8PO_INVERSE computes the inverse of a matrix factored by R8PO_FA.
%
%  Discussion:
%
%    The R8PO storage format is appropriate for a symmetric positive definite
%    matrix and its inverse.  (The Cholesky factor of a R8PO matrix is an
%    upper triangular matrix, so it will be in R8GE storage format.)
%
%    Only the diagonal and upper triangle of the square array are used.
%    This same storage scheme is used when the matrix is factored by
%    R8PO_FA, or inverted by R8PO_INVERSE.  For clarity, the lower triangle
%    is set to zero.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 February 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real R(N,N), the Cholesky factor, in R8GE storage, returned by R8PO_FA.
%
%    Output, real B(N,N), the inverse matrix, in R8PO storage.
%
  b(1:n,1:n) = r(1:n,1:n);
%
%  Compute Inverse ( R ).
%
  for k = 1 : n

    b(k,k) = 1.0E+00 / b(k,k);
    b(1:k-1,k) = -b(1:k-1,k) * b(k,k);

    for j = k + 1 : n
      t = b(k,j);
      b(k,j) = 0.0E+00;
      b(1:k,j) = b(1:k,j) + t * b(1:k,k);
    end

  end
%
%  Compute Inverse ( R ) * ( Inverse ( R ) )'.
%
  for j = 1 : n

    for k = 1 : j - 1
      t = b(k,j);
      b(1:k,k) = b(1:k,k) + t * b(1:k,j);
    end

    b(1:j,j) = b(1:j,j) * b(j,j);

  end

  return
end


text:
function C = rdivide(A, B)
%./    CHEBOP right division.
%   C = A./B, where ones of a A is a CHEBOP and is B a scalar is equivalent
%   to A*(1/B), repectively.
%
%   All other instances returns an error.
%
% See also CHEBOP/MTIMES, CHEBOP/MRDIVIDE.

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

if ( isnumeric(B) )
    C = mtimes(A, 1./B);
else
    error('CHEBOP:RDIVIDE:NotSupported', '%s./%s is not supported.', ...
        class(A), class(B));
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function errors = test_gsp_estimate_vertex_time_psd

errors = 0;

%Time parameter
T = 200;fs=1;
% Graph parameters
N = 100;
%Graph
G = gsp_sensor(N);
% G1 = gsp_2dgrid(10);
G = gsp_jtv_graph(G,T,fs);
G = gsp_compute_fourier_basis(G);

%wave filter
alpha=1;
beta = 0.7;
[g, ft] = gsp_jtv_design_damped_wave(G, alpha,beta);

Nx = 5;
x2 = (rand(N,T,1,Nx)>0.99).*randn(N,T,1,Nx);
X2 = gsp_jtv_filter_synthesis(G,g,ft,x2);

param.L = 100;
param.use_fast = 0;
t1 = tic;
param.estimator = 'TVA';
psd_TVA1 = gsp_estimate_vertex_time_psd(G,X2,param);
time1 = toc(t1)

t2 = tic;
param.estimator = 'TVA';
param.use_fast = 1;
psd_TVA2 = gsp_estimate_vertex_time_psd(G,X2,param);
time2 = toc(t2)

A1 = gsp_filter_evaluate(psd_TVA1,G.e);
A2 = gsp_filter_evaluate(psd_TVA2,G.e);

gsp_assert_test(A1,A2,1e-10,'ESTIMATE_PSD: TVA')

end

repo: gspbox


text:
function err = save_nifti(hdr,niftifile)
% err = save_nifti(hdr,niftifile)
%
% Pixel data should be in hdr.vol
%
% Handles data structures with more than 32k cols by setting
% hdr.dim(2) = -1 and hdr.glmin = ncols. This is FreeSurfer specific,
% for handling surfaces. The exception to this is when the total
% number of spatial voxels equals 163842, then the volume is
% reshaped to 27307x1x6xnframes. This is for handling the 7th
% order icosahedron used by FS group analysis.
%
% $Id: save_nifti.m,v 1.14 2011/10/17 19:58:37 nicks Exp $

%
% save_nifti.m
%
% Original Author: Doug Greve
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/10/17 19:58:37 $
%    $Revision: 1.14 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

err = 1;
if(nargin ~= 2)
  fprintf('err = save_nifti(hdr,niftifile)\n');
  return;
end

ext = niftifile((strlen(niftifile)-2):strlen(niftifile));
if(strcmpi(ext,'.gz'))
  gzip_needed = 1;
  niftifile = niftifile(1:strlen(niftifile)-3);
  %fprintf('First, saving to %s before compressing\n',niftifile);
else
  gzip_needed = 0;
end

% Check for ico7
sz = size(hdr.vol);
if(sz(1) == 163842)
  %fprintf('save_nifti: ico7 reshaping\n');
  dim = [27307 1 6 size(hdr.vol,4)];
  hdr.vol = reshape(hdr.vol, dim);
end

fp = fopen(niftifile,'w');
if(fp == -1)
  fprintf('ERROR: could not open %s\n',niftifile);
  return;
end

hdr.data_type = [hdr.data_type(:)' repmat(' ',[1 10])];
hdr.data_type = hdr.data_type(1:10);

hdr.db_name = [hdr.db_name(:)' repmat(' ',[1 18])];
hdr.db_name = hdr.db_name(1:18);

hdr.dim = ones(1,8);
hdr.dim(1) = 4;
hdr.dim(2) = size(hdr.vol,1);
hdr.dim(3) = size(hdr.vol,2);
hdr.dim(4) = size(hdr.vol,3);
hdr.dim(5) = size(hdr.vol,4);

% This is to accomodate structures with more than 32k cols
% FreeSurfer specific. See also mriio.c.
if(hdr.dim(2) > 2^15)
  hdr.glmin = hdr.dim(2);
  hdr.dim(2) = -1;
end

hdr.pixdim = [hdr.pixdim(:)' repmat(0,[1 8])];
hdr.pixdim = hdr.pixdim(1:8);

hdr.descrip = [hdr.descrip(:)' repmat(' ',[1 80])];
hdr.descrip = hdr.descrip(1:80);

hdr.aux_file = [hdr.aux_file(:)' repmat(' ',[1 24])];
hdr.aux_file = hdr.aux_file(1:24);

hdr.intent_name = [hdr.intent_name(:)' repmat(' ',[1 16])];
hdr.intent_name = hdr.intent_name(1:16);

hdr.magic = [hdr.magic(:)' repmat(0,[1 4])];
hdr.magic = hdr.magic(1:4);

if(isempty(hdr.regular))  hdr.regular  = ' '; end
if(isempty(hdr.dim_info)) hdr.dim_info = ' '; end
if(isempty(hdr.slice_code)) hdr.slice_code = ' '; end
if(isempty(hdr.xyzt_units)) hdr.xyzt_units = ' '; end % should be err

hdr.vox_offset = 352; % not 348
fwrite(fp,348,'int');
fwrite(fp,hdr.data_type,    'char');
fwrite(fp,hdr.db_name,      'char');
fwrite(fp,hdr.extents,      'int');
fwrite(fp,hdr.session_error,'short');
fwrite(fp,hdr.regular,      'char');
fwrite(fp,hdr.dim_info,     'char');
fwrite(fp,hdr.dim,          'short');
fwrite(fp,hdr.intent_p1,    'float');
fwrite(fp,hdr.intent_p2,    'float');
fwrite(fp,hdr.intent_p3,    'float');
fwrite(fp,hdr.intent_code,  'short');
fwrite(fp,hdr.datatype,     'short');
fwrite(fp,hdr.bitpix,       'short');
fwrite(fp,hdr.slice_start,  'short');
fwrite(fp,hdr.pixdim,       'float');
fwrite(fp,hdr.vox_offset,   'float');
fwrite(fp,hdr.scl_slope,    'float');
fwrite(fp,hdr.scl_inter,    'float');
fwrite(fp,hdr.slice_end,    'short');
fwrite(fp,hdr.slice_code,   'char');
fwrite(fp,hdr.xyzt_units,   'char');
fwrite(fp,hdr.cal_max,      'float');
fwrite(fp,hdr.cal_min,      'float');
fwrite(fp,hdr.slice_duration,'float');
fwrite(fp,hdr.toffset,       'float');
fwrite(fp,hdr.glmax,         'int');
fwrite(fp,hdr.glmin,         'int');
fwrite(fp,hdr.descrip,       'char');
fwrite(fp,hdr.aux_file,      'char');
fwrite(fp,hdr.qform_code,    'short');
fwrite(fp,hdr.sform_code,    'short');
fwrite(fp,hdr.quatern_b,     'float');
fwrite(fp,hdr.quatern_c,     'float');
fwrite(fp,hdr.quatern_d,     'float');
fwrite(fp,hdr.quatern_x,     'float');
fwrite(fp,hdr.quatern_y,     'float');
fwrite(fp,hdr.quatern_z,     'float');
fwrite(fp,hdr.srow_x,        'float');
fwrite(fp,hdr.srow_y,        'float');
fwrite(fp,hdr.srow_z,        'float');
fwrite(fp,hdr.intent_name,   'char');
fwrite(fp,hdr.magic,         'char');

% Pad to get to 352 bytes (header size is 348)
fwrite(fp,0,'char');
fwrite(fp,0,'char');
fwrite(fp,0,'char');
fwrite(fp,0,'char');

npix = prod(size(hdr.vol));
switch(hdr.datatype)
 case   2, nitemswritten = fwrite(fp,hdr.vol,'uchar'); % dont use char
 case   4, nitemswritten = fwrite(fp,hdr.vol,'short');
 case   8, nitemswritten = fwrite(fp,hdr.vol,'int');
 case  16, nitemswritten = fwrite(fp,hdr.vol,'float');
 case  64, nitemswritten = fwrite(fp,hdr.vol,'double');
 case 512, nitemswritten = fwrite(fp,hdr.vol,'ushort');
 case 768, nitemswritten = fwrite(fp,hdr.vol,'uint');
 otherwise,
  fprintf('WARNING: data type %d not supported, but writing as float',...
	  hdr.datatype);
  nitemswritten = fwrite(fp,hdr.vol,'float');
  return;
end
fclose(fp);

if(npix ~= nitemswritten)
  fprintf('ERROR: tried to write %d, but only wrote %d',npix,nitemswritten);
  return;
end

if(gzip_needed)
  cmd = sprintf('gzip -f %s', niftifile);
  %fprintf('Compressing with\n');
  %fprintf('   %s\n',cmd);
  unix(cmd);
end

err = 0;
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
close all
clear

pos_rob = [-0.5,-0.5;0.5,0.5];
pos_ref =  [1.1,1;-1.1,-1];
pos_circ = [-0.5,-0.5;0.5,0.5]/2;

% circular constraint sets
zbf_circ = @(i,pos_rob) 1-(pos_rob(i,:)-pos_circ(i,:))*(pos_rob(i,:)-pos_circ(i,:))';
% nominal control
knom = 1;
unom_rob = @(i,pos_rob) -knom*(pos_rob(i,:)-pos_ref(i,:));
% bounded distance costraints
ld = 0.5; ud = 1;
zbf_lb = @(i,j,pos_rob) (pos_rob(i,:)-pos_rob(j,:))*(pos_rob(i,:)-pos_rob(j,:))'-ld^2;
zbf_ub = @(i,j,pos_rob) ud^2-(pos_rob(i,:)-pos_rob(j,:))*(pos_rob(i,:)-pos_rob(j,:))';
% extended K-class function
akcf = 1;
ekcf = @(h) akcf*h;
% control barrier function (less or equal)
cbfa_circ = @(i,pos_rob) pos_rob(i,:)-pos_circ(i,:);
cbfb_circ = @(i,pos_rob) ekcf(zbf_circ(i,pos_rob));
cbfa_lb = @(i,j,pos_rob) -(pos_rob(i,:)-pos_rob(j,:));
cbfb_lb = @(i,j,pos_rob) ekcf(zbf_lb(i,j,pos_rob));
cbfa_ub = @(i,j,pos_rob) -cbfa_lb(i,j,pos_rob);
cbfb_ub = @(i,j,pos_rob) ekcf(zbf_ub(i,j,pos_rob));

% simulation
dt = 0.01;
T = 10;
loop = 0;
color_list = ['b','g','m','r','c'];
for t=0:dt:T
    loop = loop+1;
    % control barrier function - quadratic program
    for i=1:2
        fqp = -unom_rob(i,pos_rob)';
        aqp = cbfa_circ(i,pos_rob);
        bqp = cbfb_circ(i,pos_rob);
        for j=1:2
            if i~=j
                aqp = [aqp;cbfa_lb(i,j,pos_rob)];
                bqp = [bqp;cbfb_lb(i,j,pos_rob)];
                aqp = [aqp;cbfa_ub(i,j,pos_rob)];
                bqp = [bqp;cbfb_ub(i,j,pos_rob)];
            end
        end
        [uqp,~,is_solved] = quadprog(eye(2),fqp,aqp,bqp);
        uqp_rob(i,:) = uqp';
    end
    % update
    pos_rob = pos_rob+dt*uqp_rob;
    % plot
    for i=1:2
        plot(pos_rob(i,1),pos_rob(i,2),[color_list(i) 'o']); hold on
        circle_(pos_circ(i,:),1,'color',color_list(i));
        plot(pos_ref(i,1),pos_ref(i,2),[color_list(i) 'd']);
        quiver(pos_rob(i,1),pos_rob(i,2),uqp_rob(i,1),uqp_rob(i,2))
    end
    hold off
    drawnow
end

repo: paper-simulation


text:
function K = covLINone(hyp, x, z, i)

% Linear covariance function with a single hyperparameter. The covariance
% function is parameterized as:
%
% k(x^p,x^q) = (x^p'*x^q + 1)/t^2;
%
% where the P matrix is t2 times the unit matrix. The second term plays the
% role of the bias. The hyperparameter is:
%
% hyp = [ log(t) ]
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2010-09-10.
%
% See also COVFUNCTIONS.M.

if nargin<2, K = '1'; return; end                  % report number of parameters
if nargin<3, z = []; end                                   % make sure, z exists
xeqz = isempty(z); dg = strcmp(z,'diag');                       % determine mode

it2 = exp(-2*hyp);                                                  % t2 inverse

% precompute inner products
if dg                                                               % vector kxx
  K = sum(x.*x,2);
else
  if xeqz                                                 % symmetric matrix Kxx
    K = x*x';
  else                                                   % cross covariances Kxz
    K = x*z';
  end
end

if nargin<4                                                        % covariances
  K = it2*(1+K);
else                                                               % derivatives
  if i==1
    K = -2*it2*(1+K);
  else
    error('Unknown hyperparameter')
  end
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% spherror() - chancenter() sub function to compute minimum distance
%               of Cartesian coordinates to a sphere
%
% Author: Scott Makeig, CNL / Salk Institute, 2000

% Copyright (C) Scott Makeig, CNL / Salk Institute, 2000
%
% 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 2 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function wobble = spherror(center,x,y,z)

% 03/14/02 corrected wobble calculation -lf

x = x - center(1);  % center the points at (0,0,0)
y = y - center(2);
z = z - center(3);
radius = (sqrt(x.^2+y.^2+z.^2)); % distances from the center
wobble = std(radius-mean(radius)); % test if xyz values are on a sphere


text:
%% Copyright (C) 2016-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 sinhint (@var{x})
%% Symbolic sinhint function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = sinhint (x)
%%   @result{} y = (sym) Shi(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 = sinhint(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('Shi', x);
end

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

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

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

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

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

repo: octsympy