text:
function [H_filt, h_filt] = arraySHTfiltersMeas_regLSHD(H_array, order_sht, grid_dirs_rad, w_grid, nFFT, amp_threshold)
%ARRAYSHTFILTERSMEAS_REGLSHD Generate SHT filters based on measured responses
%(regularized least-squares in the SHD)
%
%   Generate the filters to convert microphone signals from a spherical
%   microphone array to SH signals, based on a least-squares solution with
%   a constraint on noise amplification, using Tikhonov regularization. The
%   method formulates the LS problem in the spherical harmonic domain, by
%   first expanding the array response to SH coefficients, and then solving
%   the LS problem, similar to
%
%       Jin, C.T., Epain, N. and Parthy, A., 2014.
%       Design, optimization and evaluation of a dual-radius spherical microphone array.
%       IEEE/ACM Transactions on Audio, Speech, and Language Processing, 22(1), pp.193-204.
%
%   Inputs:
%       H_array:    nFFT x nMics x nMeasurement_dirs matrix of measured
%           responses, in the frequency domain with length nFFT
%       order_sht:  order of SH signals to generate
%       grid_dirs:  nMeasurement_dirs x 2 matrix of [azi elev] of the
%           measurement directions, in rads
%       w_grid:     nMeasurement_dirs x 1 vector of weights for
%           weighted-least square solution, based on the importance or area
%           around each measurement point (leave empty if not known, or
%           not important)
%       nFFT:       number of FFT points in the responses
%       amp_threshold:  max allowed amplification for filters, in dB
%
%   Outputs:
%       H_filt: nSH x nMics x (nFFT/2+1) returned filters in the frequency
%           domain (half spectrum up to Nyquist)
%       h_filt: nFFT x nMics x nSH impulse responses of the above filters
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% ARRAYSHTFILTERSMEAS_REGLSHD.M - 5/10/2016
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nGrid = size(H_array,3);
nMics = size(H_array,2);
if order_sht>sqrt(nMics)-1
    warning('Set order too high for the number of microphones, should be N<=sqrt(Q)-1')
    order_sht = floor( sqrt(nMics)-1 );
end
nBins = nFFT/2+1;
if isempty(w_grid)
    w_grid = ones(nGrid,1);
end

% SH matrix at grid directions
order_array = floor(sqrt(nGrid)/2-1);
aziElev2aziPolar = @(dirs) [dirs(:,1) pi/2-dirs(:,2)]; % function to convert from azimuth-inclination to azimuth-elevation
Y_grid = sqrt(4*pi) * getSH(order_array, aziElev2aziPolar(grid_dirs_rad), 'real')'; % SH matrix for grid directions

% compute inverse matrix
a_dB = amp_threshold;
alpha = 10^(a_dB/20);
beta = 1/(2*alpha);
W_grid = diag(w_grid);
H_filt = zeros((order_sht+1)^2, nMics, nBins);
% compute first the SHT of the array response
for kk=1:nBins
    tempH = squeeze(H_array(kk,:,:));
    H_nm(kk,:,:) = tempH * W_grid * Y_grid'* inv(Y_grid*W_grid*Y_grid');
end
% compute the inverse matrix in the SHD with regularization
for kk=1:nBins
    tempH_N = squeeze(H_nm(kk,:,:));
    tempH_N_trunc = tempH_N(:,1:(order_sht+1)^2);
    H_filt(:,:,kk) = tempH_N_trunc' * inv(tempH_N*tempH_N' + beta^2*eye(nMics));
end

if nargout>1
    % time domain filters
    h_filt = H_filt;
    h_filt(:,:,end) = abs(h_filt(:,:,end));
    h_filt = cat(3, h_filt, conj(h_filt(:,:,end-1:-1:2)));
    h_filt = real(ifft(h_filt, [], 3));
    h_filt = fftshift(h_filt, 3);
end

author: polarch

repo: Spherical-Array-Processing


text:
classdef TestNiBlackThreshold
    %TestNiBlackThreshold

    properties (Constant)
        im = fullfile(mexopencv.root(),'test','sudoku.jpg');
    end

    methods (Static)
        function test_1
            img = cv.imread(TestNiBlackThreshold.im, 'Grayscale',true);
            out = cv.niBlackThreshold(img, -0.2);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_2
            img = cv.imread(TestNiBlackThreshold.im, 'Grayscale',true);
            out = cv.niBlackThreshold(img, 0.2, ...
                'Type','Binary', 'BlockSize',7, 'MaxValue',255);
            validateattributes(out, {class(img)}, {'size',size(img)});
        end

        function test_error_argnum
            try
                cv.niBlackThreshold();
                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 pass = test_periodicbvp
% TAD, 10 Jan 2014

tol = 1e-8;

%% Building blocks
dom = [-pi pi];
[Z, I, D, C] = linop.primitiveOperators(dom);
[z, E, s] = linop.primitiveFunctionals(dom);
x = chebfun('x', dom);
c = sin(x.^2);
C = operatorBlock.mult(c);
El = E(dom(1));
Er = E(dom(end));

%% Solve a linear system
L = [ D^2, -I, sin(x); C, D, chebfun(0,dom); functionalBlock.zero(dom), El, 4 ];
f = [sin(x); chebfun(0,dom); 1 ];
L = addbc(L,'periodic');

%%

type = {@chebcolloc2, @chebcolloc1, @ultraS, @chebcolloc2, @chebcolloc1, @ultraS};
prefs = cheboppref;
prefs.bvpTol = 1e-14;

w = [];
for k = 1:6
    prefs.discretization = type{k};
    w = linsolve(L,f,prefs);

    %%
    % check the ODEs
    w1 = w{1};  w2 = w{2};  w3 = w{3};
    f1 = f{1};  f2 = f{2}; f3 = f{3};
    residual1 = diff(w1,2) - w2 + sin(x).*w3 - f1;
    err(k,1) = norm( residual1 );
    residual2 = c.*w1 + diff(w2) + 0 - f2;
    err(k,2) = norm( residual2 );
    err(k,3) = abs( 0 + feval(w2, dom(1)) + 4*w3 - f3 );

    %%
    % check the BCs
    v = w{2};  u = w{1};
    Du = D*u;  Dv = D*v;
    err(k,4) = abs( u(pi) - u(-pi) );
    err(k,5) = abs( v(pi) - v(-pi) );
    err(k,6) = abs( Du(pi) - Du(-pi) );

    %%
    % check continuity
    err(k,7) = feval(u, pi/2, 'left') - feval(u, pi/2, 'right');
    err(k,8) = feval(v, pi/2, 'left') - feval(v, pi/2, 'right');
    err(k,9) = feval(Du, pi/2, 'left') - feval(Du, pi/2, 'right');
    err(k,10) = feval(u, -pi/2, 'left') - feval(u, -pi/2, 'right');
    err(k,11) = feval(v, -pi/2, 'left') - feval(v, -pi/2, 'right');
    err(k,12) = feval(Du, -pi/2, 'left') - feval(Du, -pi/2, 'right');

    if ( k == 2 )
        % introduce breakpoints
        f = [abs(cos(x)); 0*x; 1 ];
    end

end

%% Test for TRIGCOLLOC.

% Domain.
dom = [0 2*pi];

% Differentiation operators.
D = operatorBlock.diff(dom);
D2 = operatorBlock.diff(dom, 2);

% Multiplication operators.
a1 = chebfun(@(x) 1 + sin(2*x), dom);
A1 = operatorBlock.mult(a1);
a0 = chebfun(@(x) 1 + cos(x), dom);
A0 = operatorBlock.mult(a0);

% Linop.
L = linop(D2 + A1*D + A0);

% Rhs.
f = chebfun(@(x) cos(x), dom);

% Solve it with TRIGCOLLOC.
prefs.discretization = @trigcolloc;
u = linsolve(L, f, prefs);
u = u{1};

err(1, 13) = norm(diff(u,2) + a1.*diff(u) + a0.*u - f);
err(2, 13) = abs(u(2*pi) - u(0));
err(3, 13) = abs(feval(diff(u), 2*pi) - feval(diff(u), 0));

%%
pass = err < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function M = euclideanfactory(m, n)
% Returns a manifold struct to optimize over m-by-n matrices.
%
% function M = euclideanfactory(m, n)
%
% Returns M, a structure describing the Euclidean space of m-by-n matrices,
% equipped with the standard Frobenius distance and associated trace inner
% product, as a manifold for Manopt.
%
% m and n in general can be vectors to handle multidimensional arrays.
% If either of m or n is a vector, they are concatenated as [m, n].
%
% Using this simple linear manifold, Manopt can be used to solve standard
% unconstrained optimization problems, for example in replacement of
% Matlab's fminunc.
%
% See also: euclideancomplexfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors: Bamdev Mishra, May 4, 2015.
% Change log:
%
%   July 5, 2013 (NB):
%       Added egred2rgrad, ehess2rhess, mat, vec, tangent.
%   May 4, 2015 (BM):
%       Added functionality to handle multidimensional arrays.

    % The size can be defined using both m and n, or simply with m.
    % If m is a scalar, then n is implicitly 1.
    % This mimicks the use of built-in Matlab functions such as zeros(...).
    if ~exist('n', 'var') || isempty(n)
        if numel(m) == 1
            n = 1;
        else
            n = [];
        end
    end

    dimensions_vec = [m(:)', n(:)']; % We have a row vector.

    M.name = @() sprintf('Euclidean space R^(%s)', num2str(dimensions_vec)); % BM: okay.

    M.dim = @() prod(dimensions_vec);% BM: replacing m*n;

    M.inner = @(x, d1, d2) d1(:).'*d2(:); % BM: okay.

    M.norm = @(x, d) norm(d(:), 'fro');% BM: replacing norm(d, 'fro');

    M.dist = @(x, y) norm(x(:) - y(:), 'fro');% BM: replacing norm(x-y, 'fro');

    M.typicaldist = @() sqrt(prod(dimensions_vec));% BM: replacing sqrt(m*n);

    M.proj = @(x, d) d; % BM: okay.

    M.egrad2rgrad = @(x, g) g; % BM: okay.

    M.ehess2rhess = @(x, eg, eh, d) eh; % BM: okay.

    M.tangent = M.proj;

    M.exp = @exp;
    function y = exp(x, d, t)
        if nargin == 3
            y = x + t*d; % BM: okay.
        else
            y = x + d; % BM: okay.
        end
    end

    M.retr = M.exp;

	M.log = @(x, y) y-x; % BM: okay.

    M.hash = @(x) ['z' hashmd5(x(:))]; % BM: okay.

    M.rand = @() randn(dimensions_vec);% BM: replacing randn(m, n);

    M.randvec = @randvec;
    function u = randvec(x) %#ok<INUSD>
        u = randn(dimensions_vec);% BM: replacing randn(m, n);
        u = u / norm(u(:), 'fro');% BM: replacing u / norm(u, 'fro');
    end

    M.lincomb = @matrixlincomb;

    M.zerovec = @(x) zeros(dimensions_vec);% BM: replacing zeros(m, n);

    M.transp = @(x1, x2, d) d;

    M.pairmean = @(x1, x2) .5*(x1+x2); % BM: okay.

    M.vec = @(x, u_mat) u_mat(:); % BM: okay.
    M.mat = @(x, u_vec) reshape(u_vec, dimensions_vec);% BM: replacing reshape(u_vec, [m, n]);
    M.vecmatareisometries = @() true;

end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function [result] = smallMatrix(size)
    result =10^(-50)*ones(size, 'double', 'gpuArray');
end

author: jiweil

repo: Hierarchical-Neural-Autoencoder


text:
function x = p24_start ( n )

%*****************************************************************************80
%
%% P24_START returns a starting point for optimization for problem 24.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables X.
%
%    Output, real X(N), a starting point for the optimization.
%
  x = ( linspace ( -5.12, +5.12, n ) )';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function gf=gradloglikGaP(x, varargin)
% f=loglikGaP(x, varargin)
% complete log likellihood of the GaP model
% Vxt = varargin{1};      %data
% sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
% alpha = varargin{3}; %parameters of the Gamma prior on the blinking
% beta = varargin{4}; %parameters of the Gamma prior on the blinking
% peval = varargin{5}; %parameters
% x(1:end-2*peval.ncomp) is Hkt

Vxt = varargin{1};      %data
sigpsf = varargin{2};  %std deviation of the PSF gaussian approx
alpha = varargin{3}; %parameters of the Gamma prior on the blinking
beta = varargin{4}; %parameters of the Gamma prior on the blinking
peval = varargin{5}; %parameters

% % % Hkt_linear=x(1:end-peval.ncomp*2); % intensities
cx=x(end-peval.ncomp*2+1:end-peval.ncomp); %x-coordinates of the centers
cy=x(end-peval.ncomp+1:end); % y-coordinates of the centers
% % %
% % % Hkt_linear=x; % intensities
Hkt_linear=varargin{6};
% % % cy=varargin{7};
% % %

sigpsf_vec=repmat(sigpsf,peval.ncomp,1); %all psfs same sigma
cxy_vec=[cx'+1, cy'+1];
% cxy_vec=[cx', cy'];
a_vec=1./(sigpsf_vec.^2*2*pi); % all normalised to 1

Hkt=reshape(Hkt_linear, peval.ncomp, peval.nt);
% generate PSFs from given parameters:
Wxkpix=gauss2dmultislice([peval.nx, peval.ny, peval.ncomp], cxy_vec, sigpsf_vec, a_vec);
Wxkpix=normalizePSF(Wxkpix); %normalize PSFs to 1
Wxk=reshape(Wxkpix,peval.nx*peval.ny, peval.ncomp);
[Wxkbg,Hktbg]=addbg(Wxk, Hkt, peval.bg);
P=Wxkbg*Hktbg; %current approximation

%linear grasdient shifted by cx
xxvc = lineargrad([peval.nx, peval.ny, peval.ncomp], cx, 'xx');
yyvc = lineargrad([peval.nx, peval.ny, peval.ncomp], cy, 'yy');

% dW/dcx:
Wxtcx=1/sigpsf^2*xxvc.*Wxk;
% dW/dcy:
Wxtcy=1/sigpsf^2*yyvc.*Wxk;

% d(log(L))/dHkt:
% % % gfHkt=(alpha-1)*1./Hkt - 1/beta +
% Wxk'*(Vxt./P-eye(peval.nx*peval.ny,peval.nt));
gfHkt= Wxk'*(Vxt./P)-1; %without background (->not Wxkgb) and d(log(L)/dcx)
% d(log(L))/dcx:
gfcx=diag(Wxtcx'*(Vxt./P-ones(peval.nx*peval.ny, peval.nt))*Hkt');

% d(log(L))/dcy:
gfcy=diag(Wxtcy'*(Vxt./P-ones(peval.nx*peval.ny, peval.nt))*Hkt');

% % % gf = [reshape(gfHkt',1,peval.nt*peval.ncomp), gfcx', gfcy'];
%gf = [reshape(gfHkt,1,peval.nt*peval.ncomp)];
gf = [gfcx', gfcy'];
end
function Wnorm=normalizePSF(W)
sw=size(W);
Wr=reshape(W, sw(1)*sw(2),sw(3));
q=squeeze(sum(Wr,1));
Wrnorm=Wr./repmat(q,sw(1)*sw(2),1);
Wnorm=reshape(Wrnorm,sw(1), sw(2), sw(3));
end
function xxvc = lineargrad(sizevec, cx, dir)
switch dir
    case 'xx'
%         xxp=double(xx(sizevec, 'true')); %linear function - pixels
        xxp=double(xx(sizevec, 'corner')); %linear function - pixels
    case 'yy'
%         xxp=double(yy(sizevec, 'true')); %linear function - pixels
        xxp=double(yy(sizevec, 'corner')); %linear function - pixels
    otherwise
        error('Wrong dir')
end

xxv=reshape(xxp,sizevec(1)*sizevec(2),sizevec(3)); %linear function - vector
xxvc=xxv-repmat(cx,sizevec(1)*sizevec(2),1);
% xxp=double(xx([peval.nx, peval.ny, peval.ncomp], 'true')); %linear function - pixels
%yyp=double(yy([peval.nx, peval.ny, peval.ncomp], 'true'));
% xxv=reshape(xxp,peval.nx*peval.ny,peval.ncomp); %linear function - vector
%yyv=reshape(yyp,peval.nx*peval.ny,peval.ncomp);
% xxvc=xxv-repmat(cx,peval.ncomp,1);
%yyvc=yyv-repmat(cy,peval.ncomp,1);
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function DEM_demo_song_priors
% Demo for a bird songs: In this example, we simulate local field potential
% using the prediction error from the song-bird example below. We look at
% these responses under natural stimuli and after removing the second
% level of the hierarchy to show it is necessary for veridical perception.
% We then repeat but omitting dynamical priors by forsaking generalised
% coordinates
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEM_demo_song_priors.m 6030 2014-05-31 13:09:24Z karl $

% hierarchical non-linear generative model (dynamic & chaotic)
%==========================================================================

% timing
%--------------------------------------------------------------------------
N        = 128;                      % length of stimulus (bins)
dt       = 1/64;                     % time bin (seconds)

% correlations
%--------------------------------------------------------------------------
M(1).E.s = 1;
M(1).E.K = exp(-2);

% level 1
%--------------------------------------------------------------------------
% P(1): Prandtl number
% P(2): 8/3
% P(3): Rayleigh number

P        = [10; 8/3];
x        = [0.9; 0.8; 30];
f        = '[-P(1) P(1) 0; (v(1) - 4 - x(3)) -1 0; x(2) 0 -P(2)]*x/16;';
M(1).f   = f;
M(1).g   = inline('x([2 3])','x','v','P');
M(1).x   = x;
M(1).pE  = P;
M(1).V   = exp(1);
M(1).W   = exp(8);

% level 2
%--------------------------------------------------------------------------
P        = [10; 8/3];
x        = [0.9; 0.8; 30];
f        = '[-P(1) P(1) 0; (32 - x(3)) -1 0; x(2) 0 -P(2)]*x/128;';
M(2).f   = f;
M(2).g   = inline('x(3)','x','v','P');
M(2).x   = x;
M(2).pE  = P;
M(2).V   = exp(8);
M(2).W   = exp(8);

% create data
%==========================================================================

% create innovations & add causes
%--------------------------------------------------------------------------
DEM      = spm_DEM_generate(M,N);

% DEM estimation and display
%==========================================================================
DEM.M(1).x = [1; 1; 8];
DEM.M(2).x = [1; 1; 8];

% canoncial DEM
%--------------------------------------------------------------------------
DEMc   = DEM;

% without second level
%==========================================================================
DEMa   = DEM;
DEMa.M = DEMa.M(1);

% without generlised coordinates
%==========================================================================
DEMb   = DEM;
DEMb.M(1).E.n = 1;

% deconvolve
%--------------------------------------------------------------------------
DEMa   = spm_DEM(DEMa);
DEMb   = spm_DEM(DEMb);
DEMc   = spm_DEM(DEMc);

% show songs and prediction error (ERP)
%==========================================================================
spm_DEM_qU(DEMc.qU,DEMc.pU)
colormap('pink')

% Sonograms
%--------------------------------------------------------------------------
subplot(3,2,5)
spm_DEM_play_song(DEMc.pU ,N*dt);
title('simulus','Fontsize',18)
axis square

subplot(3,2,6)
spm_DEM_play_song(DEMc.qU ,N*dt);
title('percept','Fontsize',18)
axis square
drawnow

% Sonograms
%--------------------------------------------------------------------------
spm_figure('GetWin','Figure 1');
clf, colormap('pink')

subplot(3,2,1)
spm_DEM_play_song(DEMc.qU ,N*dt);
title('percept (right click on image)','Fontsize',16)

subplot(3,2,3)
spm_DEM_play_song(DEMa.qU,N*dt);
title('no structural priors','Fontsize',16)

subplot(3,2,5)
spm_DEM_play_song(DEMb.qU,N*dt);
title('no dynamical priors','Fontsize',16)

% LFPs
%--------------------------------------------------------------------------
subplot(3,2,2)
spm_DEM_EEG(DEMc,dt,[1 2],1);
title('LFP','Fontsize',16)

subplot(3,2,4)
spm_DEM_EEG(DEMa,dt,[1 2],1);
title('LFP','Fontsize',16)

subplot(3,2,6)
spm_DEM_EEG(DEMb,dt,[1 2],1);
title('LFP','Fontsize',16)
drawnow, disp(' '),disp('Click sonograms to play songs'),disp(' ')

author: spm


text:
function d =  testing(a,d)

  %K = get_kernel(a.child,d,a.X); %% OLD CODE

  if a.use_kernels
    K = test(a.child, d);  K=K.X;
  else
    K=d.X;
  end

  if(a.use_b & ~a.use_kernels) K=[K ; ones(1,size(K,2))]; end;
  Yest = K'*a.alpha;

  if a.algorithm.use_signed_output==1
    Yest=sign(Yest);
  end
  d=set_x(d,Yest);
  d=set_name(d,[get_name(d) ' -> ' get_name(a)]);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function kcoeff = kwiener_train(X,Y,s)
% FUNCTION kcoeff = kwiener_train(X,Y,param)
%   AUTHOR:    Makoto Yamada
%              (myamada@ism.ac.jp)
%   DATE:       12/25/07
%
%  DESCRIPTION:
%
%   This function compute the kernel Wiener filter (kernel Dependency
%   Estimation) using Canonical Correlation Analysis (CCA)
%   and outputs it in the abstract object "kcoeff".
%   The "kcoeff" object can then be used for finding pre-images.
%
%  INPUTS:
%
%            X:     "X" is a (n times N) dimensional input matrix,
%                   where "n" is the vector dimension and N is the number
%                    of samples.
%
%            Y:     "Y" is a (m times N) dimensional output matrix.
%                    where "m" is the vector dimension and N is the number
%                    of samples.
%
%            s:      "s" is a parameter for Gaussian kernel,
%                    exp(-norm(x - y)^2/s).
%
%  OUTPUTS:
%
%    kcoeff   :     The Kernel Wiener Filter object.
%
%
%  Example:
%              load usps; %loading USPS data
%              s = 256*0.7; %Gaussian Kernel Parameter
%              %Training Kernel Wiener Filter
%              kcoeff = kwiener_train(X,Y,s);
%
%  Reference:  1. J. Weston et.al. ``Kernel Dependency Estimation,'' NIPS
%                 2003
%
%              2. M. Yamada and M. R. Azimi-Sadjadi, ``Kernel Wiener Filter
%                 using Canonical Correlation Analysis Framework,''
%                 IEEE SSP'05, Bordeaux, France, July 17-20, 2005.
%
%              3. C. Cortes et.al. ``A General Regression Technique for
%                 Learning Transductions,'' ICML, Bonn, Germany, Aug 7-11.
%
%              4. M. Yamada and M. R. Azimi-Sadjadi, ``Kernel Wiener Filter
%                 with Distance Constraint,'' ICASSP,  Toulouse, France,
%                 May 14-19, 2006
%
%              5. Zhe Chen et.al, ``Correlative Learning: A Basis for Brain and Adaptive
%                 Systems,'' Wiley, Oct. 2007 (Section 4.6)

n = size(X,2); %Number of Training samples

%Kernel Computation (Gaussian Kernel)
disp('Computing the Gram Matrix of X');
Kx = kernelg(X,X,s);
disp('Computing the Gram Matrix of Y');
Ky = kernelg(Y,Y,s);

%Centering Matrix
J = eye(n) - 1/n*ones(n);
%Centering Kernel Gram Matrix
A = J*Kx*J;

%Computing the Kernel Wiener Filter (Kernel Dependency Estimation)
%Computing the first kernel Canonical Coordinate.
[F D] = eig(A);
Lambda = sqrt(abs(real(diag(D))));
[YY I]  = sort(-Lambda);

%Normalizing Eigenvectors.
for k = 1:n
    if Lambda(k) ~= 0
        F(:,k)=F(:,k)/Lambda(k);
    end
end

Lambda = -YY;
F      = F(:,I);
kcoeff.Df   = F;                 %First kernel Canonical mapping function.
kcoeff.lambda = Lambda;          %Eigenvalue of kernel CCA.

%Computing the second kernel Canonical Coordinate.
ry = rank(Ky);
[Uy Wy Vy] = svd(Ky);
invKy = Vy(:,1:ry)*inv(Wy(1:ry, 1:ry))*Uy(:,1:ry)';
kcoeff.Dg = invKy*Kx*J*kcoeff.Df; %Second kernel Canonical mapping function.

kcoeff.Kxt = kcoeff.Df'*Kx;

%Computing rank of kernel Wiener filter.
[temp1,temp2] = find(kcoeff.lambda > 0);

if sum(temp2) == size(kcoeff.Kxt,2)
    kcoeff.rnk = sum(temp2) -1;
else
    kcoeff.rnk = sum(temp2);
end

disp('Finished Kernel Wiener Estimation');

author: Sable

repo: mcbench-benchmarks


text:
function [idx, new_mpc, new_xgd, new_sd] = addstorage(storage, mpc, xgd, sd)
%ADDWIND   Adds storage units and corresponding xGen/StorageData to existing data.
%
%   [IDX, NEW_MPC] = ADDSTORAGE(STORAGE, MPC)
%   [IDX, NEW_MPC, NEW_XGD, NEW_SD] = ADDSTORAGE(STORAGE, MPC)
%   [IDX, NEW_MPC, NEW_XGD, NEW_SD] = ADDSTORAGE(STORAGE, MPC, XGD)
%   [IDX, NEW_MPC, NEW_XGD, NEW_SD] = ADDSTORAGE(STORAGE, MPC, XGD, SD)
%
%   Given a StorageUnitData structure, or the name of a file containing such
%   a structure, this function adds the specified storage generators to an
%   existing MATPOWER case, xGenData struct and StorageData struct.
%
%   Inputs:
%       STORAGE : a StorageUnitData struct or the name of an M-file
%                   or MAT-file that returns one, with the following fields
%           .gen     : rows to be appended to the GEN matrix from MPC
%           .gencost : (optional) rows to be added to the GENCOST matrix
%                      from MPC, default is zero cost
%           .xgd_table : xGenData table struct or filename providing data for
%                        the storage units being added. See LOADSTORAGEDATA
%                        for more information on the xGenData table format.
%           .sd_table : StorageData table struct or filename providing data
%                       for the storage units being added. See LOADSTORAGEDATA
%                       for more information on the StorageData table format.
%       MPC : MATPOWER case struct to which storage generators will be added
%       XGD : (optional) xGenData struct corresponding to the generators
%             already in MPC, to which the new xGenData for the storage
%             units will be added.
%       SD : (optional) StorageData struct corresponding to the generators
%            already in MPC, to which the new StorageData for the storage
%            units will be added.
%
%   Output:
%       IDX : Generator indices of newly added storage units.
%       NEW_MPC : MPC with storage units appended to MPC.GEN and MPC.GENCOST
%                 MPC.GENFUEL (= 'ess').
%       NEW_XGD : XGD with xGenData for new storage units appended.
%       NEW_SD  : SD with StorageData for new storage units appended.
%
%   See also LOADSTORAGEDATA, LOADXGENDATA.

%   MOST
%   Copyright (c) 2013-2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MOST.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/most for more info.

%% define named indices into data matrices
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;

%% define fuel type for storage
STORAGE_FUEL = 'ess';

%% input arg handling
if nargin < 4
    sd = [];
    if nargin < 3
        xgd = [];
    end
end
if ischar(storage)
    infile = sprintf(' in file: ''%s''', storage);
    storage = loadgenericdata(storage, 'struct', 'gen', 'storage');
else
    infile = '';
end

%% add to MPC
ns = size(storage.gen, 1);          %% number of storage units being added
if isfield(storage, 'gencost')
    storage_gencost = storage.gencost;
else        %% use zero cost by default
    storage_gencost = repmat([POLYNOMIAL 0 0 2 0 0], ns, 1);
end
[new_mpc, idx] = addgen2mpc(mpc, storage.gen, storage_gencost, STORAGE_FUEL);

%% handle xGenData and StorageData
if nargout > 2      %% output NEW_XGD, NEW_SD requested
    %% xGenData
    if isfield(storage, 'xgd_table')
        storage_xgd = loadxgendata(storage.xgd_table, storage.gen);
    else
        error('addstorage: missing XGD_TABLE field in STORAGE');
    end

    if isempty(xgd)     %% no input XGD provided
        new_xgd = storage_xgd;
    else                %% input XGD provided
        %% append rows of every field in xgd
        new_xgd = xgd;
        fields = fieldnames(xgd);
        for f = 1:length(fields)
            ff = fields{f};
            %% dims of wind_xgd fields already checked by loadxgendata
            if size(xgd.(ff), 1) ~= size(mpc.gen, 1)
                error('addstorage: # of rows in XGD.%s (%d) does not match MPC.GEN (%d)', ...
                    ff, size(xgd.(ff), 1), size(mpc.gen, 1));
            end
            new_xgd.(ff) = [xgd.(ff); storage_xgd.(ff)];
        end
    end

    %% StorageData
    if isfield(storage, 'sd_table')
        storage_sd = loadstoragedata(storage.sd_table, storage.gen);
    else
        error('addstorage: missing SD_TABLE field in STORAGE');
    end

    storage_sd.UnitIdx = idx;       %% add UnitIdx field

    if isempty(sd)     %% no input SD provided
        new_sd = storage_sd;
    else                %% input SD provided
        %% find number of storage units in original mpc
        ns0 = length(find(strcmp(mpc.genfuel, STORAGE_FUEL)));
        new_sd = sd;
        fields = fieldnames(sd);
        for f = 1:length(fields)    %% append rows of every field in sd
            ff = fields{f};
            add_ns  = size(storage_sd.(ff), 1);
            orig_ns = size(        sd.(ff), 1);
            if orig_ns ~= 1 || add_ns ~= 1 || sd.(ff) ~= storage_sd.(ff)
                %% at least one of them is not 1-d or else they don't match
                if orig_ns == 1     %% original is 1-d
                    %% expand it
                    sd.(ff) = ones(ns0, 1) * sd.(ff);
                end
                if add_ns == 1      %% new is 1-d
                    %% expand it
                    storage_sd.(ff) = ones(ns, 1) * storage_sd.(ff);
                end
                %% concatenate them
                new_sd.(ff) = [sd.(ff); storage_sd.(ff)];
            % else (both are 1-d and they match)
            %   do nothing, just use old sd.(ff)
            end
        end
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
% Choose a png file from the internet that represents something from your
% home country, your life or just something you generally like, and plot
% its polyline silhouette.

author: odedstein

repo: sgi-introduction-course


text:
% Originally Coded By: Niall Mangan
% Copyright 2017, All Rights Reserved
% Code by Niall Mangan for paper "Inferring biological networks by sparse
% identification of nonlinear dynamics"
% by N. M. Mangan S. L. Brunton, J. L. Proctor, and J. N. Kutz
% Original Code: https://github.com/niallmm/iSINDy
%%
% This file is modified to compare the performance of i-SINDy on
% insufficient data.
% Modified by:K, 2019/07/16
% Find equations for yeast glycolysis state variable 7
%%
clc;clear all; close all;

addpath('./utils');
addpath('./bioutils');

% define libarary parameters
laurentorder = 0;
polyorder = 3;
usesine = 0;
dyorder = 1;

% clear variables from other solutions.

clear Theta Thetastring Xi indTheta lambdavec numterms errorv indopt
clear indTheta1 Xi1 numterms1 nT

% Here we load the previously simulated data
load('TrainingData.mat')

% Define the new data length
percent=0.005;new_length=round(percent*length(xt));

% Shuffel the original data
Sequence=randperm(size(xt,1));
xt=xt(Sequence,:);
dxt=dxt(Sequence,:);

% Assign the value to the new variables
Data=xt(1:new_length,:);dData=dxt(1:new_length,:);

% Define the number of states
n=size(Data,2);

% pool Data  (i.e., build library of nonlinear time series)
[Theta, Thetastring] = poolDatady(Data,n,polyorder,usesine, laurentorder, dData(:,7), dyorder);
% %initial lambda value, which is the value used for soft thresholding in ADM

tol = 2e-3;
lambda = 8e-3;

jj = 1; % counter
num= 1; % initialize the number of nonzero terms found for the lambda
errorvec= 0;

MaxIter = 1e3;

% for now calculate null space using null function
nT = null(Theta);

% Define the parameters for plooting
plottag=2;

[indTheta1, Xi1, numterms1] = ADMinitvary(nT,lambda,MaxIter,tol, plottag);

Thetastring(indTheta1)'
n0Xi = Xi1(Xi1~=0);
n0Xi/n0Xi(end)
save('Results/7th_state_variable.mat')

author: dynamicslab

repo: SINDy-PI


text:
function p = sigmoid(log_odds)
% SIGMOID: Inverse of the logit function.
%          This is a one-to-one mapping from log odds to probability.
%          i.e. it maps the real line to the interval (0,1).
%
%   p = sigmoid(log_odds)

assert(nargin==1)

p = 1 ./ (1 + exp(-log_odds));

author: nesl


text:
function bmu_colors=som_bmucolor(bmus, m, colors)

% SOM_BMUCOLOR Returns the colors of the bmus according to a map colorcode
%
% bmu_colors=som_bmucolor(bmus, msize, colors);
%
% INPUT ARGUMENTS ([]'s are optional)
%
% bmus   (matrix) Nx1 vector of BMU indexes
% msize  (map struct, topol struct or 1x2 vector)
%          gives the map grid size
% colors (matrix) colormap(s): munits x 3 x d matrix of RGB vectors
%
% OUTPUT ARGUMENTS
%
% bmu_colors (Nx3xd matrix) color of the data point according to its BMU's
%              color(s).
%
% Idea is to get a color for each data point that links it to its BMU.
%
% EXAMPLE
%
% We want to show how an time series is projected  to a map. Instead of
% a trajectory, we use 'color linking'
%
% map=som_make(multi_dim_signal);
% bmus=som_bmu(map,multi_dim_signal);
% Colors=som_bmucolor(bmus, map, som_colorcode(map,'rgb1'));
% colorsignal(Colors, multi_dim_signal);
%
% See also SOM_COLORCODE.

% Copyright (c) 1999-2000 by the SOM toolbox programming team.
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0alpha Johan 170699

%% Check arguments %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

error(nargchk(3, 3, nargin))   % check no. of input args is correct

% Check map grid size

if vis_valuetype(m,{'1x2'}),
  msize=m;
else
  [tmp,ok,tmp]=som_set(m);
  if isstruct(m) && all(ok)        % check m type
    switch m.type
    case 'som_topol'
      msize=m.msize;
      lattice=m.lattice;
    case 'som_map'
      msize=m.topol.msize;
      lattice=m.topol.lattice;
    otherwise
      error('Invalid map or topol struct.');
    end
  end
end

if length(msize)>2
  error('Only 2D maps allowed!');
end

n=prod(msize)

% Check colorcode size

if ~vis_valuetype(colors,{'nx3xdimrgb','nx3rgb'})
  error('Colorcode matrix not valid!');
end

% Check bmu vector

if ~vis_valuetype(bmus,{'nx1'}),
  error('Need a column vector of BMU indexes!');
else
  bmus=round(bmus);
  if max(bmus) > n || min(bmus) < 1
    error('BMU indexes exeed the map size!')
  end
end

%% Action %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

bmu_c=colors(bmus,:,:);

%% Build output %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

bmu_colors=squeeze(bmu_c);

author: ilarinieminen

repo: SOM-Toolbox


text:
function pts = cuboid_prm2pts(prm, vp)
% converts five parameters of cuboid to eight 2D points
%
% prm: [x, y, w, h, d]
% x,y: 2D coordinates of corner1 (top,left,front corner)
% w: length (in pixels) from corner1 to corner2
% h: length (in pixels) from corner1 to corner3
% d: length (in pixels) from corner1 to corner5
%
% pts: [8 x 2] 2D coordinates of eight points on cuboid
%   5---6   1------2
%  /   /|   |      |
% 1---2 8   3------4    ....
% |   |/     \    /
% 3---4       7--8
%
% 1: front,top,left
% 2: front,top,right
% ...
% 8: back,bottom,right

pts = zeros(8,2);

pts(1,:) = [prm(1) prm(2)];
pts(2,:) = towardsright(pts(1,:), vp, prm(3));
pts(3,:) = towardsdown(pts(1,:), vp, prm(4));
pts(5,:) = towardsback(pts(1,:), vp, prm(5));
pts(4,:) = line_intersect(pts(2,:), vp{1}, pts(3,:), vp{2});
pts(6,:) = line_intersect(pts(2,:), vp{3}, pts(5,:), vp{2});
pts(7,:) = line_intersect(pts(3,:), vp{3}, pts(5,:), vp{1});
pts(8,:) = mean(...
    [line_intersect(pts(4,:), vp{3}, pts(6,:), vp{1}); ...
     line_intersect(pts(7,:), vp{2}, pts(6,:), vp{1}); ...
     line_intersect(pts(4,:), vp{3}, pts(7,:), vp{2})], 1);

author: zouchuhang

repo: LayoutNet


text:
function pred = ml_predicthdca(trials, model)
% Prediction function for Hierarchical Discriminant Component Analysis
% Prediction = ml_predicthdca(Trials, Model)
%
% In:
%   Trials  : the data a matrix, as in ml_predict
%
%   Model   : predictive model as produced by ml_trainhdca
%
% Out:
%   Prediction  : discrete probability distribution, formatted as
%                 {'disc' [NxC] [Cx1]}, with element #2 being the per-class probability and
%                 element #3 the original target values per class
%                 thus, the expected target values are Prediction{2}*Prediction{3}
%
% Examples:
%   targets might look like this: [-1 -1 1 -1 1 -1 -1 1 -1 -1 1 -1 -1 1 ...]'
%
%   model = ml_trainhdca(data,targets)
%   p = ml_predicthdca(data, model); expectation = p{2}*p{3};
%   now expectation might look like this: [-0.6 -0.9 +0.4 -0.7 +0.8 -0.1 +0.5 +1.0 -0.9 +1.0 -1.0 -1.0 +1.0 ...]'
%
% See also:
%   ml_trainhdca
%
%                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                           2010-04-03

if isfield(model,'voted')
    % dispatch to the voter for multi-class classification
    pred = ml_predictvote(trials,model);
else
    % for each modality range...
    for m=1:length(model.modality_ranges)
        range = model.modality_ranges{m};

        % for each sub-block...
        for b=size(trials,2):-1:1
            % extract per-block features
            blocktrials = reshape(trials(range,b,:),[],size(trials,3))';
            blockpredictions{b} = blocktrials*model.blockmodels{m}{b}.w' - model.blockmodels{m}{b}.b';
        end
        % extract per-modality features
        layertrials = cat(2,blockpredictions{:});
        rangepredictions{m} = layertrials*model.rangemodels{m}.w' - model.rangemodels{m}.b'; %#ok<AGROW>
    end
    % perform top-level prediction
    toptrials = cat(2,rangepredictions{:});
    pred = ml_predictlda(toptrials,model.topmodel);
end

repo: REST


text:
function varargout = grErode(varargin)
%GRERODE Morphological erosion on graph.
%
%   LBL2 = grErode(EDGES, LBL1)
%   Each label of the graph is assigned the smallest label of its
%   neighbours, or it keeps the same label this one is smaller.
%
%   Example
%   grErode
%
%   See also
%   grDilate, grOpen, grClose

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2006-01-20
% Copyright 2006-2022 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas)

if length(varargin) == 2
    edges   = varargin{1};
    lbl     = varargin{2};
elseif length(varargin) == 3
    edges   = varargin{2};
    lbl     = varargin{3};
else
    error('Wrong number of arguments in "grErode"');
end

lbl2 = zeros(size(lbl));

uni = unique(edges(:));
for n = 1:length(uni)
    neigh = grAdjacentNodes(edges, uni(n));
    lbl2(uni(n)) = min(lbl([uni(n); neigh]));
end

varargout{1} = lbl2;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function dat=func_ar(dat, order, varargin)
% func_ar (Feature extraction) :
%
% This function calculates the autoregression(AR) parameter.
%
% Example:
% [out] = func_ar(dat, 7, {'method','arburg'})
%
% Returs:
%     dat    - Data structure, segmented
%     order  - Order of AR setting
% Option: models for obtatining AR parameter
%     method - 'aryule'(default), 'arburg', 'arcov', 'armcov'
%

opt=opt_cellToStruct(varargin{:});
opt=struct('method',opt.method);

if isempty(dat)
    warning('[OpenBMI] Warning! data is empty.');
end

if isempty(order)
    warning('[OpenBMI] Order is not exist.');
end

if isempty(opt.method) %method selection
   opt.method='aryule';
end

[T, nEvents , nChans]= size(dat.x);

temp_ar= [];
for i= 1:nChans*nEvents,
  ar= feval(opt.method, dat.x(:,i), order);
  temp_ar(:,i)= ar(2:end)';
end

dat.x= temp_ar;

author: PatternRecognition

repo: OpenBMI


text:
function Y = TT_weingarten(V, Z, ind)
% Weingarten map computation for the fixed TT-rank manifold.
%
% NOTE: this manifold requires the use of a modified version of TTeMPS_1.1,
% which is packaced with Manopt and can be found in
% manopt/manopt/manifolds/ttfixedrank/TTeMPS_1.1
%
% Y = TT_weingarten(V, Z, ind)
%
% This function implements the efficient way of computing the Weingarten
% map for the Riemannian submanifold of fixed TT-rank tensors embedded in
% its natural embedding Euclidean space, as described in the paper:
%
%   Psenka and Boumal,
%   Second-order optimization for tensors with fixed tensor-train rank,
%   Optimization workshop at NeurIPS 2020.
%
% Y is the output tangent vector, V is a tangent vector (which contains
% needed information of base point X), Z is a tensor (embedding space).
%
% TT_weingarten supports the following formats for Z:
% 1) Full
% 2) TT-format (any rank)
% 3) sparse, with non-zero indices indexed by ind (see fixedTTrankfactory).
%
% See also: fixedTTrankfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Michael Psenka, Nov. 24, 2020.
% Contributors: Nicolas Boumal
% Change log:

    % Preliminary tangent vector set-up for Y
    d = V.order;
    r = V.rank;
    n = V.size;

    Y = V; % all properties except for dU cores inherited from V
    normV = norm(V);
    V = (1 / normV) * V;
    Y.dU = cell(1, d);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Begin calculating variational components
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Pre-calculate all R-components as needed for tangent space conversion
    R = cell(1, d-1);
    % Intermediary used to calculate R-terms for re-orthogonalization
    Xr = V.U;

    for k = d:-1:2

        sz = size(Xr{k});

        if length(sz) == 2
            sz = [sz, 1]; %#ok<AGROW>
        end

        [Q, R{k - 1}] = qr_unique(unfold(Xr{k}, 'right')');
        Xr{k} = reshape(Q', [size(Q, 2), sz(2), sz(3)]);

        Xr{k - 1} = tensorprod_ttemps(Xr{k - 1}, R{k - 1}, 3);
    end

    % Calculate V.dU under the original parametrization of tangent space

    dUR = cell(1, d);

    for k = 1:(d - 1)
        dUR{k} = unfold(V.dU{k}, 'left') / R{k}';
        dUR{k} = reshape(dUR{k}, [r(k), n(k), r(k + 1)]);
    end

    dUR{d} = V.dU{d};

    % Calculate ~X^tV_{\ge k} efficiently
    XtVg = cell(1, d);

    XtVg{d} = conj(unfold(V.V{d}, 'right')) * unfold(V.dU{d}, 'right').';
    XtVgTmp = conj(unfold(V.V{d}, 'right')) * unfold(V.U{d}, 'right').';

    for k = (d - 1):-1:2
        tmp = tensorprod_ttemps(V.V{k}, XtVg{k + 1}', 3);
        XtVg{k} = conj(unfold(tmp, 'right')) * unfold(V.U{k}, 'right').';

        tmp2 = tensorprod_ttemps(V.V{k}, XtVgTmp', 3);
        XtVg{k} = XtVg{k} + conj(unfold(tmp2, 'right')) * unfold(dUR{k}, 'right').';
        XtVgTmp = conj(unfold(tmp2, 'right')) * unfold(V.U{k}, 'right').';
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%% CASES FOR TYPE OF EUCLIDEAN GRADIENT Z %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if isa(Z, 'TTeMPS') || isa(Z, 'TTeMPS_tangent_orth')% Z either TT-tensor or tangent vec

        if isa(Z, 'TTeMPS_tangent_orth')
            Z = tangent_to_TTeMPS(Z); % efficient to treat both as TTeMPS case
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%               DIAGONAL TERMS                 %%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        ZVr = cell(1, d); % stores values of V_<^T(Z<k>)
        Zr = cell(1, d);  % stores values of X_<^T(Z<k>)
        ZVl = cell(1, d); % stores values of (Z<k>)^T(V_>)
        Zl = cell(1, d);  % stores values of (Z<k>)^T(X_>)

        ZVr{d}  = conj(unfold(Z.U{d}, 'right')) * unfold(V.dU{d}, 'right').';
        XtVgTmp = conj(unfold(Z.U{d}, 'right')) * unfold(V.U{d},  'right').';

        for k = (d - 1):-1:2
            tmp = tensorprod_ttemps(Z.U{k}, ZVr{k + 1}', 3);
            ZVr{k} = conj(unfold(tmp, 'right')) * unfold(V.U{k}, 'right').';

            tmp2 = tensorprod_ttemps(Z.U{k}, XtVgTmp', 3);
            ZVr{k} = ZVr{k} + conj(unfold(tmp2, 'right')) * unfold(dUR{k}, 'right').';
            XtVgTmp = conj(unfold(tmp2, 'right')) * unfold(V.U{k}, 'right').';
        end

        for k=1:d
            ZVr{k} = ZVr{k}';
        end

        Zr{d} = conj(unfold(V.V{d}, 'right')) * unfold(Z.U{d}, 'right').';

        for k = (d - 1):-1:2
            tmp = tensorprod_ttemps(V.V{k}, Zr{k + 1}', 3);
            Zr{k} = conj(unfold(tmp, 'right')) * unfold(Z.U{k}, 'right').';
        end

        % Computation of ZVl and Zl

        ZVl{1} = unfold(dUR{1}, 'left')' * unfold(Z.U{1}, 'left');
        Zl{1}  = unfold(V.U{1}, 'left')' * unfold(Z.U{1}, 'left');

        for k = 2:(d - 1)
            % (V_k-1)(U_k)
            tmp = tensorprod_ttemps(V.U{k}, ZVl{k - 1}', 1);
            ZVl{k} = unfold(tmp, 'left')' * unfold(Z.U{k}, 'left');
            % + (X_k-1)(dU_k)
            tmp = tensorprod_ttemps(dUR{k}, Zl{k - 1}', 1);
            ZVl{k} = ZVl{k} + unfold(tmp, 'left')' * unfold(Z.U{k}, 'left');
            % update Zr to keep up w/ recursive definition
            tmp = tensorprod_ttemps(V.U{k}, Zl{k - 1}', 1);
            Zl{k} = unfold(tmp, 'left')' * unfold(Z.U{k}, 'left');
        end

        ZZ = cell(1, d); % final cores for normal (I \otimes X)Z(X^T)
        Zv = cell(1, d); % final cores for (I \otimes X)Z(V^T)
        vZ = cell(1, d); % final cores for (I \otimes V)Z(X^T)

        % contract to first core
        ZZ{1} = tensorprod_ttemps(Z.U{1}, Zr{2}, 3);
        % contract to inner cores
        for k = 2:(d - 1)
            res = tensorprod_ttemps(Z.U{k}, Zl{k - 1}, 1);
            ZZ{k} = tensorprod_ttemps(res, Zr{k + 1}, 3);
        end

        % contract to last core
        ZZ{d} = tensorprod_ttemps(Z.U{d}, Zl{d - 1}, 1);

        Zv{1} = tensorprod_ttemps(Z.U{1}, ZVr{2}, 3);

        for k = 2:(d - 1)
            res = tensorprod_ttemps(Z.U{k}, Zl{k - 1}, 1);
            Zv{k} = tensorprod_ttemps(res, ZVr{k + 1}, 3);
        end

        Zv{d} = tensorprod_ttemps(Z.U{d}, Zl{d - 1}, 1);

        vZ{1} = tensorprod_ttemps(Z.U{1}, Zr{2}, 3);

        for k = 2:(d - 1)
            res = tensorprod_ttemps(Z.U{k}, ZVl{k - 1}, 1);
            vZ{k} = tensorprod_ttemps(res, Zr{k + 1}, 3);
        end

        vZ{d} = tensorprod_ttemps(Z.U{d}, ZVl{d - 1}, 1);

        % Left unfold everything to apply additional operations
        for k = 1:d
            ZZ{k} = unfold(ZZ{k}, 'left');
            Zv{k} = unfold(Zv{k}, 'left');
            vZ{k} = unfold(vZ{k}, 'left');
        end
        % %%%%%%%%%%%%%%%%%%%%%%%%% TESTING %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % 1st and last cases special, so they're computed outside the loop
        UL  = unfold(V.U{1}, 'left');
        dUL = unfold(dUR{1}, 'left');

        % right variational term without (I - UU^T) projection
        rightV = (Zv{1} - ZZ{1} * XtVg{2}) / R{1};

        Y.dU{1} = -dUL * UL' * ZZ{1} + rightV - UL * (UL' * rightV);

        for k = 2:(d - 1)
            UL  = unfold(V.U{k}, 'left');
            dUL = unfold(dUR{k}, 'left');

            % right variational term without (I - UU^T) projection
            rightV = (Zv{k} - ZZ{k} * XtVg{k + 1}) / R{k};

            Y.dU{k} = vZ{k} - UL * (UL' * vZ{k}) - dUL * UL' * ZZ{k} ...
                       + rightV - UL * (UL' * rightV);
        end

        Y.dU{d} = vZ{d};

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%                  CROSS TERMS                 %%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        for k = 1:(d - 1)
            U = unfold(V.U{k}, 'left');
            ZZ{k} = ZZ{k} - U * (U' * ZZ{k});
            ZZ{k} = reshape(ZZ{k}, [r(k), n(k), r(k + 1)]);
        end

        ZZ{d} = reshape(ZZ{d}, [r(d), n(d), 1]);

        % Double loop to represent P_k DP_m, 1 <= k < d
        for k = 1:(d - 1)

            UL = unfold(V.U{k}, 'left');
            dUL = unfold(dUR{k}, 'left');

            for m = 1:(k - 1)
                VtZ = crossTermVariational(k, m, ZZ, dUR{m}, V);
                projUVtZ = VtZ - UL * (UL' * VtZ);

                Y.dU{k} = Y.dU{k} - projUVtZ;
            end

            for m = (k + 1):d
                ZtX = crossTermMatrixRight(k + 1, m, ZZ, V);
                Y.dU{k} = Y.dU{k} + dUL * ZtX;
            end

        end

        % Final terms for k = d
        for m = 1:(d - 1)
            VtZ = crossTermVariational(d, m, ZZ, dUR{m}, V);
            Y.dU{d} = Y.dU{d} - VtZ;
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%FINAL RESHAPE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        for k = 1:d
            Y.dU{k} = reshape(Y.dU{k}, [r(k), n(k), r(k + 1)]);
        end

    elseif ~exist('ind', 'var')% Z is a full array

        %{
        Note that you can split up the full expression for the Weingarten map in
        the following way:

        1) Diagonal terms(P_i DP_i terms, or non - cross terms)
        a) This can be split into(left variational) + (right variational),
        which correspond to two terms in product rule expansion
        2) Cross terms(P_i DP_j, we show mathematically equivalent expression
        in companion paper.)
        a) This likewise can be split into(i < j) + (i > j)
        %}

        % Cell that represents the LEFT VARIATIONAL SIDE
        Zv = cell(1, d);
        % Cell that represents the RIGHT VARIATIONAL SIDE
        vZ = cell(1, d);

        % Initialization step
        Zv{d} = Z(:);
        vZ{d} = Z(:);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%               DIAGONAL TERMS                 %%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        %%%%%%%%%%%%%%%               LEFT VARIATIONAL               %%%%%%%%%%%%%%%%%%

        % First, calculate where the right side is the variational matrix
        % Note that the calculation here is split between VR^-1 and XX^tVR^-1
        % in (I - XX^t)VR^{-1} through xx and xx2 terms respectively

        % (d-1) term done separately to initialize the temp variable handoff
        zz = reshape(Z, [prod(n(1:(d - 1))), n(d)]);
        xx = transpose(unfold(V.dU{d}, 'right'));
        xxTmp = transpose(unfold(V.U{d}, 'right')); % temp var to help calc xx terms
        xx2 = transpose(unfold(V.V{d}, 'right'));

        % Z-matrix-k-variational, represents the xx calculation for Zv{k}
        Zmkv = zz * xx;
        vZ{d - 1} = zz * xx2;
        Zv{d - 1} = (Zmkv - vZ{d - 1} * XtVg{d}) / R{d - 1};

        % auxillery term for computing anything with variational matrix
        ZZtmp = zz * xxTmp;

        % Main calcualtion for right size of Zv
        for k = (d - 2):-1:1
            zz = reshape(Zmkv, [prod(n(1:k)), n(k + 1) * r(k + 2)]);
            xx = transpose(unfold(V.U{k + 1}, 'right'));
            Zmkv = zz * xx;

            zzTmp = reshape(ZZtmp, [prod(n(1:k)), n(k + 1) * r(k + 2)]);
            xxTmp = transpose(unfold(dUR{k + 1}, 'right'));
            Zmkv = Zmkv + zzTmp * xxTmp;
            ZZtmp = zzTmp * xx;

            % NOTE: vZ is only used here to store these values for later use
            % It stores the standard right side calc for the orthogonal projector
            zz2 = reshape(vZ{k + 1}, [prod(n(1:k)), n(k + 1) * r(k + 2)]);
            xx2 = transpose(unfold(V.V{k + 1}, 'right'));
            vZ{k} = zz2 * xx2;

            Zv{k} = (Zmkv - vZ{k} * XtVg{k + 1}) / R{k};
        end

        % Standard left side calculations for Zv
        for k = 2:d

            for i = 1:(k - 1)
                Z_i = reshape(Zv{k}, [r(i) * n(i), prod(n(i + 1:k)) * r(k + 1)]);
                X_i = unfold(V.U{i}, 'left');
                Zm = X_i' * Z_i;
                Zv{k} = Zm;
            end

            Zv{k} = reshape(Zv{k}, [r(k) * n(k), r(k + 1)]);
        end

        % Final U core projector mult for Zv
        for k = 1:(d - 1)
            U = unfold(V.U{k}, 'left');
            Zv{k} = Zv{k} - U * (U' * Zv{k});
        end

        %%%%%%%%%%%%%%%               RIGHT VARIATIONAL               %%%%%%%%%%%%%%%%%%

        % represents standard projection, also needed as intermediary for this section
        ZZ = vZ;
        % Start left side calculation for vZ
        for k = 2:d

            % Do first in loop separately to initialize intermediary variable
            Z_i = reshape(vZ{k}, [n(1), prod(n(2:k)) * r(k + 1)]);
            dU_i = unfold(dUR{1}, 'left');
            Zm = dU_i' * Z_i;
            vZ{k} = Zm;

            X_i = unfold(V.U{1}, 'left');
            ZZ{k} = X_i' * Z_i;

            for m = 2:(k - 1)
                Z_i = reshape(vZ{k}, [r(m) * n(m), prod(n(m + 1:k)) * r(k + 1)]);
                X_i = unfold(V.U{m}, 'left');
                Zm = X_i' * Z_i;

                Z_tmp = reshape(ZZ{k}, [r(m) * prod(n(m)), prod(n(m + 1:k)) * r(k + 1)]);
                dU_i = unfold(dUR{m}, 'left');

                vZ{k} = Zm + dU_i' * Z_tmp;

                ZZ{k} = X_i' * Z_tmp;
            end

            vZ{k} = reshape(vZ{k}, [r(k) * n(k), r(k + 1)]);
            ZZ{k} = reshape(ZZ{k}, [r(k) * n(k), r(k + 1)]);
        end

        % Final U core projector mult for vZ

        % 1st core separate since vZ{1} = 0 up until now. Purely for efficiency
        U  = unfold(V.U{1},  'left');
        dU = unfold(V.dU{1}, 'left');
        vZ{1} = -dU * U' * ZZ{1};

        for k = 2:(d - 1)
            U  = unfold(V.U{k},  'left');
            dU = unfold(V.dU{k}, 'left');
            vZ{k} = vZ{k} - U * (U' * vZ{k});

            vZ{k} = vZ{k} - dU * U' * ZZ{k};

        end

        for k = 1:d - 1
            Y.dU{k} = Zv{k} + vZ{k};
        end

        Y.dU{d} = vZ{d};

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%                  CROSS TERMS                 %%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        for k = 1:(d - 1)
            % now only need ZZ as dU cores of normal projection of V
            U = unfold(V.U{k}, 'left');
            ZZ{k} = ZZ{k} - U * (U' * ZZ{k});
            ZZ{k} = reshape(ZZ{k}, [r(k), n(k), r(k + 1)]);
        end

        ZZ{d} = reshape(ZZ{d}, [r(d), n(d), 1]);

        % Double loop to represent P_k DP_m, 1 <= k < d
        for k = 1:(d - 1)

            UL  = unfold(V.U{k},  'left');
            dUL = unfold(V.dU{k}, 'left');

            for m = 1:(k - 1)
                VtZ = crossTermVariational(k, m, ZZ, dUR{m}, V);

                projUVtZ = VtZ - UL * (UL' * VtZ);

                Y.dU{k} = Y.dU{k} - projUVtZ;
            end

            for m = (k + 1):d
                ZtX = crossTermMatrixRight(k + 1, m, ZZ, V);
                Y.dU{k} = Y.dU{k} + dUL * ZtX;
            end

        end

        % Final terms for k = d
        for m = 1:(d - 1)
            VtZ = crossTermVariational(d, m, ZZ, dUR{m}, V);

            Y.dU{d} = Y.dU{d} - VtZ;
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FINAL RESHAPE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        for k = 1:d
            Y.dU{k} = reshape(Y.dU{k}, [r(k), n(k), r(k + 1)]);
        end

    else
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SPARSE CASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % permuted U, V, and dU cells for the efficient computation in C
        CU = cell(1, d);
        CV = cell(1, d);
        CdUR = cell(1, d);

        for k = 1:d
            CU{k} = permute(V.U{k}, [1 3 2]);
            CV{k} = permute(V.V{k}, [1 3 2]);
            CdUR{k} = permute(dUR{k}, [1 3 2]);
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%               DIAGONAL TERMS                 %%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        [ZZ, vZ, Zv] = weingarten_omega(n, r, CU, CV, CdUR, ind.', Z);

        for k = 1:d

            ZZ{k} = reshape(ZZ{k}, [r(k), r(k + 1), n(k)]);
            ZZ{k} = unfold(permute(ZZ{k}, [1 3 2]), 'left');

            vZ{k} = reshape(vZ{k}, [r(k), r(k + 1), n(k)]);
            vZ{k} = unfold(permute(vZ{k}, [1 3 2]), 'left');

            Zv{k} = reshape(Zv{k}, [r(k), r(k + 1), n(k)]);
            Zv{k} = unfold(permute(Zv{k}, [1 3 2]), 'left');
        end

        % 1st and last cases special, so they're computed outside the loop
        UL  = unfold(V.U{1}, 'left');
        dUL = unfold(dUR{1}, 'left');

        % right variational term without (I - UU^T) projection
        rightV = (Zv{1} - ZZ{1} * XtVg{2}) / R{1};
        Y.dU{1} = -dUL * UL' * ZZ{1} + rightV - UL * (UL' * rightV);

        for k = 2:(d - 1)
            UL  = unfold(V.U{k}, 'left');
            dUL = unfold(dUR{k}, 'left');

            % right variational term without (I - UU^T) projection
            rightV = (Zv{k} - ZZ{k} * XtVg{k + 1}) / R{k};
            Y.dU{k} = vZ{k} - UL * (UL' * vZ{k}) - dUL * UL' * ZZ{k} + rightV - UL * (UL' * rightV);
            % norm(vZ{k} - UL * (UL' * vZ{k}) + (Zv{k} - UL * (UL' * Zv{k})) / R{k}) / norm(V)
            % norm(V)
        end

        Y.dU{d} = vZ{d};

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%                  CROSS TERMS                 %%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        for k = 1:(d - 1)
            U = unfold(V.U{k}, 'left');
            ZZ{k} = ZZ{k} - U * (U' * ZZ{k});
            ZZ{k} = reshape(ZZ{k}, [r(k), n(k), r(k + 1)]);
        end

        ZZ{d} = reshape(ZZ{d}, [r(d), n(d), 1]);

        % Double loop to represent P_k DP_m, 1 <= k < d
        for k = 1:(d - 1)

            UL  = unfold(V.U{k}, 'left');
            dUL = unfold(dUR{k}, 'left');

            for m = 1:(k - 1)
                % XtZ = crossTermMatrixLeft(k, m, ZZ, V);
                VtZ = crossTermVariational(k, m, ZZ, dUR{m}, V);
                projUVtZ = VtZ - UL * (UL' * VtZ);

                Y.dU{k} = Y.dU{k} - projUVtZ; %+ dUL * XtZ;

                % Y.dU{k} = Y.dU{k} - (projUk * kron(eye(n(k)), Vllk') - dUkL * Xlek') * Zkm;
            end

            for m = (k + 1):d
                ZtX = crossTermMatrixRight(k + 1, m, ZZ, V);
                Y.dU{k} = Y.dU{k} + dUL * ZtX;

                %         Y.dU{k} = Y.dU{k} + dUkL * Zkm' * Xg{k+1};
            end

        end

        % Final terms for k = d
        for m = 1:(d - 1)
            VtZ = crossTermVariational(d, m, ZZ, dUR{m}, V);

            Y.dU{d} = Y.dU{d} - VtZ;

            %     Y.dU{d} = Y.dU{d} - kron(eye(n(d)), Vlld') * Zkm;
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%FINAL RESHAPE%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        for k = 1:d
            Y.dU{k} = reshape(Y.dU{k}, [r(k), n(k), r(k + 1)]);
        end

        Y = (normV) * Y;
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Helper methods
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Efficient ZtX term for P_k DP_m, m > k
% Zp here is the projection of Z (just dU cores)
% V is only passed for access to U and V cores
function ZtX = crossTermMatrixRight(k, m, Zp, V)
    ZtX = conj(unfold(Zp{m}, 'right')) * unfold(V.V{m}, 'right').';

    for p = (m - 1):-1:k
        tmp = tensorprod_ttemps(V.U{p}, (ZtX)', 3);
        ZtX = conj(unfold(tmp, 'right')) * unfold(V.V{p}, 'right').';
    end

end

% Efficient XtZ term for P_k DP_m, m < k
% Zp here is the projection of Z
% V is only passed for access to U and V cores
function XtZ = crossTermMatrixLeft(k, m, Zp, V) %#ok<DEFNU>

    XtZ = unfold(V.U{m}, 'left')' * unfold(Zp{m}, 'left');

    for p = (m + 1):k
        XtZ = tensorprod_ttemps(V.U{p}, (XtZ)', 1);
        XtZ = unfold(XtZ, 'left')' * unfold(V.V{p}, 'left');
    end

end

% Efficient VtZ (more accurately kron(eye(n(k)), V') * Z) term for P_k DP_m, m < k
% Zp here is the projection of Z, dURm is normally parametrized dU{m}
% V is only passed for access to U and V cores
function VtZ = crossTermVariational(k, m, Zp, dURm, V)

    VtZ = unfold(dURm, 'left')' * unfold(Zp{m}, 'left');

    for p = (m + 1):(k - 1)
        VtZ = tensorprod_ttemps(V.U{p}, (VtZ)', 1);
        VtZ = unfold(VtZ, 'left')' * unfold(V.V{p}, 'left');
    end

    VtZ = tensorprod_ttemps(V.V{k}, VtZ, 1);
    VtZ = unfold(VtZ, 'left');

end

author: NicolasBoumal

repo: manopt


text:
function varargout=rotqr2ro(varargin)
%
% For calling details please see v_rotqr2ro.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_rotqr2ro.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: rotqr2ro.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_rotqr2ro(varargin{:});
else
    v_rotqr2ro(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [T,irow,icol,DS]=CDTs(TrNum);
DatabasePath='ORL\s'; % - Path of the testing database
T=[];
DS=0;

escImages=1:TrNum; % Escape Images
for td = 1:40
       Files = dir(strcat(DatabasePath,int2str(td)));
for i = 1:size(Files,1)
  x=0;
  y=1;
    for j=1:length(escImages)
        if  strcmp(Files(i).name,strcat(int2str(escImages(j)),'.pgm'))
            x=1;
        else
            x=x;
        end
    end
    if not(strcmp(Files(i).name,'.')|strcmp(Files(i).name,'..')|strcmp(Files(i).name,'Thumbs.db')|x)
        DS = DS + 1; % Number of all images in the training database
        img = imread(strcat(DatabasePath,int2str(td),'/',Files(i).name));
        [irow icol dim] = size(img);
        if dim>1
            img=rgb2gray(img);
        end
        temp = reshape(img,irow*icol,1);   % Reshaping 2D images into 1D image vectors
    T = [T temp];
    end
end
end


text:
%%**********************************************************************
%% NTpred: Compute (dX,dy,dZ) for NT direction.
%%
%% compute SVD of Xchol*Zchol via eigenvalue decompostion of
%%     Zchol * X * Zchol' = V * diag(sv2) * V'.
%% compute W satisfying W*Z*W = X.
%%     W = G'*G,  where G = diag(sqrt(sv)) * (invZchol*V)'
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%**********************************************************************

function [par,dX,dy,dZ,coeff,L,hRd] = ...
    NTpred(blk,At,par,rp,Rd,sigmu,X,Z,Zchol,invZchol)

global schurfun schurfun_par
%%
%% compute NT scaling matrix
%%
[par.W,par.G,par.sv,par.gamx,par.gamz,par.dd,par.ee,par.ff] = ...
    NTscaling(blk,X,Z,Zchol,invZchol);
%%
%% compute schur matrix
%%
m = length(rp);
schur = sparse(m,m);
UU = []; EE = []; Afree = [];
%%
for p = 1:size(blk,1)
    pblk = blk(p,:);
    if strcmp(pblk{1},'l')
        [schur,UU,EE] = schurmat_lblk(blk,At,par,schur,UU,EE,p,par.dd);
    elseif strcmp(pblk{1},'q');
        [schur,UU,EE] = schurmat_qblk(blk,At,par,schur,UU,EE,p,par.dd,par.ee);
    elseif strcmp(pblk{1},'s')
        if isempty(schurfun{p})
            schur = schurmat_sblk(blk,At,par,schur,p,par.W);
        elseif ischar(schurfun{p})
            if ~isempty(par.permZ{p})
                Wp = par.W{p}(par.permZ{p},par.permZ{p});
            else
                Wp = par.W{p};
            end
            schurtmp = feval(schurfun{p},Wp,Wp,schurfun_par(p,:));
            schur = schur + schurtmp;
        end
    elseif strcmp(pblk{1},'u')
        Afree = [Afree, At{p}']; %#ok
    end
end
%%
%% compute rhs
%%
[rhs,EinvRc,hRd] = NTrhsfun(blk,At,par,X,Z,rp,Rd,sigmu);
%%
%% solve linear system
%%
[xx,coeff,L] = linsysolve(par,schur,UU,Afree,EE,rhs);
%%
%% compute (dX,dZ)
%%
[dX,dy,dZ] = NTdirfun(blk,At,par,Rd,EinvRc,xx,m);
%%**********************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [coeff,resp,modelCI,medianHR] = computeModelCoefficientsTime_HN(X,Y,censoring,seed)
% -------------------------------------------------------------------------
% function [coeff,resp,modelCI] = computeModelCoefficients(X,Y,imbalance,batchNum)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes the final model logistic regression coefficients
% using bootstrap resampling, and it associated multivariable model
% response and bootstrap confidence intervals. See ref. [1] for more
% details. This function uses logistic regression utilities from DREES
% <http://www.cerr.info/drees>, as well as a rounding function written by
% Francois Beauducel available at:
% <http://www.mathworks.com/matlabcentral/fileexchange/26212-round-with-significant-digits>
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% [2] Vallieres, M. et al. (2015).
% -------------------------------------------------------------------------
% INPUTS:
% - X: Matrix of size [nInst X nFeat], specifying the numerical data of the
%      features of the input features, where 'nInst' refers to the number
%      of instances in X, and 'nFeat' to the number of features in X.
%      Each column is a different feature.
% - Y: Column vector of size [nInst X 1] specifying the outcome status
%      (1 or 0) for all instances.
% - imbalance: String specifying the type of imbalance-adjustement strategy
%              employed. Either 'IABR' for imbalance-adjusted bootstrap
%              resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%              logistic regression (see ref.[2]).
% - batchNum: (optional input). If present, integer that specifies the
%             batch number for parallelization purposes.
% -------------------------------------------------------------------------
% OUTPUTS:
% - coeff: Column vector of size [nCoeff+1 X 1] specifying the final
%          logistic regression coefficients. Last entry specify the offset
%          of the model.
% - resp: Column vector of size [nInst X 1] specifying the linear
%         multivariable model response when 'coeff' is applied to 'X'.
% - modelCI: Column vector of size [nInst X 2] specifying the 95%
%            confidence interval on the multivariable model response as
%            defined by the 2.5 (modelCI(i,1)) and the 97.5 (modelCI(i,2))
%            percentiles, for the ith instance. See ref. [1] for more
%            details.
% -------------------------------------------------------------------------
% AUTHOR(S):
% - Martin Vallieres <mart.vallieres@gmail.com>
% - DREES development team <http://www.cerr.info/drees> (logistic regression)
% - Francois Beauducel (roundsd.m)
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: May 2015
% - Revision I: July 2015 (including imbalance-adjusted logistic regression)
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
% --> Copyright 2010, Joseph O. Deasy, on behalf of the DREES development team.
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
%
%    _______________________________________________________________
%
% --> Copyright (c) 2015, François Beauducel
%     All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% 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
% -------------------------------------------------------------------------

% INITIALIZATION
warning off
nBoot = 1000;
alpha = 0.05;
bound = 1; % One standard error
order = size(X,2);
coeff = zeros(order,nBoot);
tol = 10^6; % To avoid extremely large coefficients
nInst = numel(Y);

% RANDOM NUMBER GENERATOR SEED
rng(seed);

% COMPUTING OVER ALL BOOTSTRAP SAMPLES
respBoot = zeros(size(X,1),nBoot);
modelCI = zeros(size(X,1),2);
for n = 1:nBoot
    average = inf;
    while average > tol
        bootSam = ceil(nInst .* rand(nInst,1));
        Xtrain = X(bootSam,:); Ytrain = Y(bootSam,1);

        coeff_temp = coxphfit(Xtrain,Ytrain,'censoring',censoring(bootSam),'baseline',0);

        coeff_temp(isnan(coeff_temp)) = 0;
        coeff(:,n) = coeff_temp;
        average = mean(abs(coeff(:,n)));
    end
    [respBoot(:,n)] = responseCox(X,coeff(:,n));
end
SE_coeff = bound.*(std(coeff')')./sqrt(nBoot);
coeff = mean(coeff')';

% CALCULATING THE BOOTSTRAP CONFIDENCE INTERVALS OF THE FINAL MODEL
for i = 1:size(X,1)
    modelCI(i,1)= prctile(respBoot(i,:),alpha/2*100);
    modelCI(i,2)= prctile(respBoot(i,:),(1-alpha/2)*100);
end

% ROUNDING COEFFICIENTS

% According to their standard errors, as in ref. [1]
% Note: This type of rounding seems too strong and may not be desirable, as it significantly affects the specificity of models
%       --> kept here as comments before further investigations
% for i = 1:order+1
%     if coeff(i) > 0
%         coeff(i) = roundsd(coeff(i),ceil(log10(abs(coeff(i)/roundsd(SE_coeff(i),1)))) + 1);
%     else
%         coeff(i) = roundsd(coeff(i),ceil(log10(abs(coeff(i)/roundsd(SE_coeff(i),1)))));
%     end
% end

% Rounding to 4 significant digits
% --> seems to preserve the predictive properties of models
coeff = roundsd(coeff,4);

% MULTIVARIABLE RESPONSE OF THE FINAL MODEL
[resp] = responseCox(X,coeff);
medianHR = median(resp);

end

author: mvallieres

repo: radiomics


text:
function box_display_test06 ( )

%*****************************************************************************80
%
%% BOX_DISPLAY_TEST06 compares several methods that reach degree 4 in X, 8 in Y.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 January 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BOX_DISPLAY_TEST06:\n' );
  fprintf ( 1, '  Plot index sets for various anisotropic methods\n' );
  fprintf ( 1, '  that reach degree 4 in X, 8 in Y.\n' );
  fprintf ( 1, '\n' );

  m = 10;
  n = 10;
  title_string = '';
%
%  Total degree <= 8.
%
  box_display ( m, n, @(x,y)x+y<=-1, @(x,y)2*x+y<=8, title_string );
  filename = 'degree48_total.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Created "%s".\n', filename );
%
%  Maximum degree <= 8.
%
  box_display ( m, n, @(x,y)x+y<=-1, @(x,y)max(2*x,y)<=8, title_string );
  filename = 'degree48_max.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Created "%s".\n', filename );
%
%  Hyperbolic cross.
%
  box_display ( m, n, @(x,y)x+y<=-1, @(x,y)(2*x+1)*(y+1)<=9, title_string );
  filename = 'degree48_hyper.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Created "%s".\n', filename );
%
%  Clenshaw-Curtis sparse grid.
%
% box_display ( m, n, @(x,y)x+y<=-1, @cc_aniso, title_string );
% filename = 'degree48_cc.png';
% print ( '-dpng', filename );
% fprintf ( 1, '  Created "%s".\n', filename );
%
%  Smolyak: Log2(2*x) + Log2(y) <= Log2(L)
%
  box_display ( m, n, @(x,y)x+y<=-1, @(x,y)log2(max(2*x,1))+log2(max(y,1))<=log2(8.1), title_string );
  filename = 'degree48_smolyak.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Created "%s".\n', filename );

  return
end

author: johannesgerer


text:
function OutputName = Recognition(TestImage, m, A, Eigenfaces)
%-------------Project the selected test image and all of the training
%images into Eigenfaces space. Compare the Euclidean distances between them and find the
%  index of image who gets minmum Euclidean distances.
ProjectedImages = [];
Train_Number = size(A,2);
for i = 1 : Train_Number
    temp = Eigenfaces' * A(:,i); % Projection of centered images into facespace
    ProjectedImages = [ProjectedImages temp];
end

%-------------Project the test image you selected into Eigenfaces space-------------
InputImage = imread(TestImage);
temp = InputImage(:,:,1);

[irow icol] = size(temp);
InImage = reshape(temp',irow*icol,1);
Difference = double(InImage)-m;
Projected_TestImage = Eigenfaces'*Difference; % Test image feature vector

%----------------------- Calculate Euclidean distances and find the
%  index of image of minmum Euclidean distances--------------------
Euc_dist = [];
for i = 1 : Train_Number
    q = ProjectedImages(:,i);
    temp = ( norm( Projected_TestImage - q ) )^2;
    Euc_dist = [Euc_dist temp];
end

[Euc_dist_min , Recognized_index] = min(Euc_dist);
OutputName = strcat(int2str(Recognized_index),'.jpg');

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%% ===== Secondary sources ===============================================
conf = SFS_config;
conf.secondary_sources.size = 3;

% === linear ===
conf.secondary_sources.geometry = 'line';
conf.secondary_sources.number = 21;
x0 = secondary_source_positions(conf);
figure;
figsize(540,404,'px');
draw_loudspeakers(x0,conf);
axis([-2 2 -2 1]);
pause(1)
print_png('secondary_sources_linear.png');

conf.secondary_sources.logspread = 3.5;
x0 = secondary_source_positions(conf);
figure;
figsize(540,404,'px');
draw_loudspeakers(x0,conf);
axis([-2 2 -2 1]);
pause(1)
print_png('secondary_sources_linear_log.png');

% === circular ===
conf.secondary_sources.geometry = 'circle';
conf.secondary_sources.number = 56;
x0 = secondary_source_positions(conf);
figure;
figsize(540,404,'px');
draw_loudspeakers(x0,conf);
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('secondary_sources_circle.png');

% === box shaped ===
conf.secondary_sources.geometry = 'box';
conf.secondary_sources.number = 84;
x0 = secondary_source_positions(conf);
figure;
figsize(540,404,'px');
draw_loudspeakers(x0,conf);
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('secondary_sources_box.png');

% === box shaped with smoothed edges ===
conf.secondary_sources.geometry = 'rounded-box';
conf.secondary_sources.number = 84;
conf.secondary_sources.corner_radius = 0.3;
x0 = secondary_source_positions(conf);
figure;
figsize(540,404,'px');
draw_loudspeakers(x0,conf);
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('secondary_sources_rounded-box.png');

% === spherical array ===
conf.secondary_sources.geometry = 'sphere'; % or 'spherical'
conf.secondary_sources.number = 225;
x0 = secondary_source_positions(conf);
figure;
figsize(540,404,'px');
draw_loudspeakers(x0,conf);
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('secondary_sources_sphere.png');

% === arbitrary shaped arrays ===
% create a stadium like shape by combining two half circles with two linear
% arrays
% first getting a full circle with 56 loudspeakers
conf.secondary_sources.geometry = 'circle';
conf.secondary_sources.number = 56;
conf.secondary_sources.x0 = [];
x0 = secondary_source_positions(conf);
% store the first half cricle and move it up
x01 = x0(2:28,:);
x01(:,2) = x01(:,2) + ones(size(x01,1),1)*0.5;
% store the second half circle and move it down
x03 = x0(30:56,:);
x03(:,2) = x03(:,2) - ones(size(x03,1),1)*0.5;
% create a linear array
conf.secondary_sources.geometry = 'linear';
conf.secondary_sources.number = 7;
conf.secondary_sources.size = 1;
x0 = secondary_source_positions(conf);
% rotate it and move it left
R = rotation_matrix(pi/2);
x02 = [(R*x0(:,1:3)')' (R*x0(:,4:6)')'];
x02(:,1) = x02(:,1) - ones(size(x0,1),1)*1.5;
x02(:,7) = x0(:,7);
% rotate it the other way around and move it right
R = rotation_matrix(-pi/2);
x04 = [(R*x0(:,1:3)')' (R*x0(:,4:6)')'];
x04(:,1) = x04(:,1) + ones(size(x0,1),1)*1.5;
x04(:,7) = x0(:,7);
% combine everything
conf.secondary_sources.geometry = 'custom';
conf.secondary_sources.x0 = [x01; x02; x03; x04];
% if we gave the conf.secondary_sources.x0 to the secondary_source_positions
% function it will simply return the defined x0 matrix
x0 = secondary_source_positions(conf);
figure;
figsize(540,404,'px');
draw_loudspeakers(x0,conf);
axis([-2 2 -2.5 2.5]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('secondary_sources_arbitrary.png');
conf.plot.realloudspeakers = true;
figure;
figsize(540,404,'px');
draw_loudspeakers(x0,conf);
axis([-2 2 -2.5 2.5]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('secondary_sources_arbitrary_realloudspeakers.png');

%% ===== Monochromatic sound fields ======================================
% === WFS 3D ===
conf = SFS_config;
conf.dimension = '3D';
conf.secondary_sources.size = 3;
conf.secondary_sources.number = 225;
conf.secondary_sources.geometry = 'sphere';
% [P,x,y,z,x0,win] = sound_field_mono_wfs_25d(X,Y,Z,xs,src,fconf);
sound_field_mono_wfs([-2 2],[-2 2],0,[0 -1 0],'pw',800,conf);
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('sound_field_wfs_3d_xy.png');
sound_field_mono_wfs([-2 2],0,[-2 2],[0 -1 0],'pw',800,conf);
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('sound_field_wfs_3d_xz.png');
sound_field_mono_wfs(0,[-2 2],[-2 2],[0 -1 0],'pw',800,conf);
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('sound_field_wfs_3d_yz.png');
conf.resolution = 100;
sound_field_mono_wfs([-2 2],[-2 2],[-2 2],[0 -1 0],'pw',800,conf);
print_png('sound_field_wfs_3d_xyz.png');

% === WFS 2.5D ===
% simulating 2.5D WFS with circular array and a point source
conf = SFS_config;
conf.dimension = '2.5D';
conf.plot.useplot = true;
conf.plot.normalisation = 'center';
% [P,x,y,z,x0] = sound_field_mono_wfs(X,Y,Z,xs,src,f,conf);
[P,~,~,~,x0] = sound_field_mono_wfs([-2 2],[-2 2],0,[0 2.5 0],'ps',800,conf);
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('sound_field_wfs_25d.png');
% plotting WFS with all secondary sources
x0_all = secondary_source_positions(conf);
[~,idx] = secondary_source_selection(x0_all,[0 2.5 0],'ps');
x0_all(:,7) = zeros(1,size(x0_all,1));
x0_all(idx,7) = x0(:,7);
conf.plot.realloudspeakers = true;
plot_sound_field(P,[-2 2],[-2 2],0,x0_all,conf);
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('sound_field_wfs_25d_with_all_sources.png');
% simulating 2.5D NFCHOA with circular array and a plane wave
conf = SFS_config;
conf.dimension = '2.5D';
% sound_field_mono_nfchoa(X,Y,Z,xs,src,f,conf);
sound_field_mono_nfchoa([-2 2],[-2 2],0,[0 -1 0],'pw',800,conf);
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('sound_field_nfchoa_25d.png');

% === 2D local WFS with box shaped array and circular virtual array ===
X = [-1 1];
Y = [-1 1];
Z = 0;
xs = [1 -1 0];
src = 'pw';
f = 7000;
conf = SFS_config;
conf.resolution = 1000;
conf.dimension = '2D';
conf.secondary_sources.geometry = 'box';
conf.secondary_sources.number = 4*56;
conf.secondary_sources.size = 2;
conf.localwfs_vss.size = 0.4;
conf.localwfs_vss.center = [0 0 0];
conf.localwfs_vss.geometry = 'circular';
conf.localwfs_vss.number = 56;
sound_field_mono_localwfs_vss(X,Y,Z,xs,src,f,conf);
axis([-1.1 1.1 -1.1 1.1]);
print_png('sound_field_localwfs_2d.png');

% === stereo setup ===
conf = SFS_config;
conf.plot.normalisation = 'center';
x0 = [-1 2 0 0 -1 0 1;1 2 0 0 -1 0 1];
% [P,x,y,z] = sound_field_mono(X,Y,Z,x0,src,D,f,conf)
sound_field_mono([-2 2],[-1 3],0,x0,'ps',[1 1],800,conf)
print_png('sound_field_stereo.png');

%% ===== spatio-temporal snapshots of the sound field ====================
conf = SFS_config;
conf.dimension = '2.5D';
conf.plot.useplot = true;
% sound_field_imp_nfchoa(X,Y,Z,xs,src,t,conf)
[p,x,y,z,x0] = sound_field_imp_nfchoa([-2 2],[-2 2],0,[0 2 0],'ps',0.005,conf);
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('sound_field_imp_nfchoa_25d.png');
conf.plot.usedb = true;
plot_sound_field(p,[-2 2],[-2 2],0,x0,conf);
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('sound_field_imp_nfchoa_25d_dB.png');
conf.plot.useplot = false;
conf.t0 = 'source';
t_40cm = 0.4/conf.c; % time to travel 40 cm in s
t0 = 0.0005; % start time of focused source in s
[p_ps,~,~,~,x0_ps] = ...
    sound_field_imp_wfs([-2 2],[-2 2],0,[1.9 0 0],'ps',t0+t_40cm,conf);
[p_pw,~,~,~,x0_pw] = ...
    sound_field_imp_wfs([-2 2],[-2 2],0,[1 -2 0],'pw',t0-t_40cm,conf);
[p_fs,~,~,~,x0_fs] = ...
    sound_field_imp_wfs([-2 2],[-2 2],0,[0 -1 0 0 1 0],'fs',t0,conf);
plot_sound_field(p_ps+p_pw+p_fs,[-2 2],[-2 2],0,[x0_ps; x0_pw; x0_fs],conf)
hold;
scatter(0,0,'k','x');   % origin of plane wave
scatter(1.9,0,'k','o'); % point source
scatter(0,-1,'k','o');  % focused source
hold off;
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('sound_field_imp_multiple_sources_dB.png');

%% ===== custom grids ====================================================
conf = SFS_config;
conf.dimension = '3D';
conf.secondary_sources.number = 225;
conf.secondary_sources.geometry = 'sphere';
conf.resolution = 100;
conf.plot.normalisation = 'center';
X = randi([-2000 2000],125000,1)/1000;
Y = randi([-2000 2000],125000,1)/1000;
Z = randi([-2000 2000],125000,1)/1000;
sound_field_mono_wfs(X,Y,Z,[0 -1 0],'pw',800,conf);
print_png('sound_field_wfs_3d_xyz_custom_grid.png');
conf.plot.usedb = true;
conf.dimension = '2.5D';
conf.secondary_sources.number = 64;
conf.secondary_sources.geometry = 'circle';
sound_field_imp_nfchoa(X,Y,0,[0 2 0],'ps',0.005,conf);
axis([-2 2 -2 2]);
set(gca, 'XTick', -2:1:2); set(gca, 'YTick', -2:1:2);
print_png('sound_field_imp_nfchoa_25d_dB_custom_grid.png');

%% ===== modal windows ===================================================
conf = SFS_config;
conf.dimension = '2.5D';
conf.secondary_sources.number = 16;
conf.secondary_sources.geometry = 'circular';
conf.secondary_sources.size = 3;
conf.resolution = 300;
conf.plot.usedb = true;
conf.t0 = 'source';
X = [-2,2];
Y = [-2,2];
Z = 0;
conf.modal_window = 'rect';  % default
sound_field_imp_nfchoa(X,Y,Z,[0 -1 0],'pw',0,conf);
print_png('sound_field_imp_nfchoa_25d_dB_rect.png');
conf.modal_window = 'max-rE';
sound_field_imp_nfchoa(X,Y,Z,[0 -1 0],'pw',0,conf);
print_png('sound_field_imp_nfchoa_25d_dB_max-rE.png');
conf.modal_window = 'kaiser';
conf.modal_window_parameter = 1.0;
sound_field_imp_nfchoa(X,Y,Z,[0 -1 0],'pw',0,conf);
print_png('sound_field_imp_nfchoa_25d_dB_kaiser.png');
conf.modal_window = 'tukey';
conf.modal_window_parameter = 0.5;
sound_field_imp_nfchoa(X,Y,Z,[0 -1 0],'pw',0,conf);
print_png('sound_field_imp_nfchoa_25d_dB_tukey.png');

%% ===== impulse response of a spatial audio system ======================
conf = SFS_config;
conf.t0 = 'source';
X = [0 0 0];
phi = 0;
xs = [2.5 0 0];
src = 'ps';
t = (1:1000)/conf.fs*1000;
hrtf = dummy_irs(conf);
[ir,~,delay] = ir_wfs(X,phi,xs,src,hrtf,conf);
figure;
figsize(540,404,'px');
plot(t,ir(1:1000,1),'-g');
hold on;
offset = round(delay*conf.fs);
plot(t,ir(1+offset:1000+offset,1),'-b');
hold off;
xlabel('time / ms');
ylabel('amplitude');
print_png('impulse_response_wfs_25d.png');
X = [0 0 0];
head_orientation = [0 0];
xs = [2.5 0 0];
src = 'ps';
conf = SFS_config;
conf.N = 1000;
conf.t0 = 'source';
time_response_wfs(X,xs,src,conf)
axis([0 25 -0.005 0.025]);
print_png('impulse_response_wfs_25d_imp.png');

%% ===== frequency response of a spatial audio system ====================
X = [0 0 0];
head_orientation = [pi/2 0];
xs = [0 2.5 0];
src = 'ps';
hrtf = dummy_irs(conf);
conf = SFS_config;
conf.ir.usehcomp = false;
conf.wfs.usehpre = false;
[ir1,x0] = ir_wfs(X,head_orientation,xs,src,hrtf,conf);
conf.wfs.usehpre = true;
conf.wfs.hprefhigh = aliasing_frequency(x0,conf);
ir2 = ir_wfs(X,head_orientation,xs,src,hrtf,conf);
[a1,p,f] = spectrum_from_signal(norm_signal(ir1(:,1)),conf);
a2 = spectrum_from_signal(norm_signal(ir2(:,1)),conf);
figure;
figsize(540,404,'px');
semilogx(f,20*log10(a1),'-b',f,20*log10(a2),'-r');
axis([10 20000 -80 -40]);
set(gca,'XTick',[10 100 250 1000 5000 20000]);
legend('w/o pre-filter','w pre-filter');
xlabel('frequency / Hz');
ylabel('magnitude / dB');
print_png('frequency_response_wfs_25d.png');
% alternative variant
X = [0 0 0];
xs = [0 2.5 0];
src = 'ps';
conf = SFS_config;
freq_response_wfs(X,xs,src,conf);
axis([10 20000 -40 0]);
print_png('frequency_response_wfs_25d_mono.png');

author: sfstoolbox

repo: sfs-matlab


text:
function lp=dir_lpdf(x,a)
%DIR_LPDF   Log probability density function of uniform Dirichlet
%           distribution
%
%       Description:
%       LP = DIR_LPDF(X, A) returns the lpdf of Dirichlet distribution
%       with A at X
%
% Copyright (c) 2000 Aki Vehtari

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

lp=gammaln(sum(a))-sum(gammaln)+sum(log(x).^(a-1));

author: gpstuff-dev

repo: gpstuff

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


text:
%MAIN_simpleHarmonicOscillator.m
%
% This script runs a simulation of a simple harmonic oscillator

% Use: EoM_Single_Pendulum to write the equations of motion

m = 1.0;  % (kg)  mass
k = 1.0;  % (N/m) spring constant

tSpan = [0,10]; %Simulation time interval

x0 = 0.5;  % (m) initial position
v0 = 0;  % (m/s) initial velocit
z0 = [x0;v0];

userFunc = @(t,z)simpleHarmonicOscillatorDynamics(t,z,m,k);

options = odeset(...
    'AbsTol',1e-8,...
    'RelTol',1e-8,...
    'Vectorized','on');

% Run the simulation!
sol = ode45(userFunc,tSpan,z0,options);

% Break apart solution for plotting
nPlot = 1000;
time = linspace(tSpan(1),tSpan(2),nPlot);
z = deval(sol,time); %Evaluate solution from ode45 at points in time
x = z(1,:);
v = z(2,:);

% Plotting!

figure(111);clf;

subplot(2,1,1);
plot(time,x,'k-','LineWidth',2);
xlabel('time (s)')
ylabel('position (m)');

subplot(2,1,2);
plot(time,v,'k-','LineWidth',2);
xlabel('time (s)')
ylabel('velocity (m/s)');

author: MatthewPeterKelly

repo: dscTutorials


text:
%ZIPLOAD   Load compressed data file created with ZIPSAVE
%
%	[data] = zipload( filename )
%	filename: string variable that contains the name of the
%	compressed file (do not include '.zip' extension)
%	Use only with files created with 'zipsave'
%	pkzip25.exe has to be in the matlab path. This file is a compression utility
%	made by Pkware, Inc. It can be dowloaded from:	http://www.pkware.com
%	Or directly from ftp://ftp.pkware.com/pk250c32.exe, for the Windows 95/NT version.
%	This function was tested using 'PKZIP 2.50 Command Line for Windows 9x/NT'
%	It is important to use version 2.5 of the utility. Otherwise the command line below
%	has to be changed to include the proper options of the compression utility you
%	wish to use.
%	This function was tested in MATLAB Version 5.3 under Windows NT.
%	Fernando A. Brucher - May/25/1999
%
%	Example:
%		[loadedData] = zipload('testfile');
%--------------------------------------------------------------------

function [data] = zipload( filename )

%--- Decompress data file by calling pkzip (comand line command) ---
% 	Options used:
%	'extract' = decompress file
%	'silent' = no console output
%	'over=all' = overwrite files

%eval( ['!pkzip25 -extract -silent -over=all ', filename, '.zip'] )
eval( ['!pkzip25 -extract -silent -over=all ', filename, '.zip'] )

%--- Load data from decompressed file ---
%	try, catch takes care of cases when pkzip fails to decompress a
%	valid matlab format file

try
   tmpStruc = load( filename );
   data = tmpStruc.data;
catch, return, end

%--- Delete decompressed file ---

delete( [filename,'.mat'] )

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [S,D,J] = sample_discrete(D,n)
  % SAMPLE_DISCRETE  Sample a discrete distribution given by D n times. This
  % employs the "Alias Table Method". Ensuring that the running time is O(n +
  % m)
  %
  % Inputs:
  %   D  m list of probabilities summing to 1
  %   n  number of samples to output
  % Outputs:
  %   S  n list of indices into D
  %   D  m list of alias table heights
  %   J  m list of "upper half" reindices
  %

  % Assume column vector
  D = D(:);
  % https://hips.seas.harvard.edu/blog/2013/03/03/the-alias-method-efficient-sampling-with-many-discrete-outcomes/
  m = numel(D);
  assert(abs(sum(D) - 1)<1e-15);
  % We want to sort based on above or below 1/n
  D = D*m;
  % now we can sort based on above or below 1
  % Below queue: http://www.alecjacobson.com/weblog/?p=3933
  % and preallocate up to 2*m (extra m to be safe, I don't think it's needed)
  B = find(D<1);
  Bl = numel(B);
  B(end+1:2*m) = 0;
  % After queue, and preallocate up to m
  A = find(D>=1);
  Al = numel(A);
  A(end+1:m) = 0;
  % Upper half list (initialize with this index because above includes "equal")
  J = (1:m)';

  % While below and above still have elements
  while Bl>0 && Al>0
    % pop below
    b = B(Bl);
    Bl = Bl - 1;
    % Look above
    a = A(Al);
    % a going to be b's upper half
    J(b) = a;
    % lob off enough to fit D(b) to 1
    %D(a) = D(a) - (1.0 - D(b));
    % Reduce rounding error (see comments above)
    D(a) = (D(a) + D(b)) - 1.0;
    if D(a) < 1
      % pop off above
      Al = Al - 1;
      % push onto below
      Bl = Bl+1;
      B(Bl) = a;
    % else still on above
    end
  end

  R1 = floor(rand(n,1)*m)+1;
  R2 = rand(n,1);
  I = (1:m)';
  % select those that are in upper half
  upper = R2>=D(R1);
  S = I(R1);
  S(upper) = J(R1(upper));

end

author: alecjacobson

repo: gptoolbox


text:
function cdf_discrete_test ( )

%*****************************************************************************80
%
%% CDF_DISCRETE_TEST tests CDF_DISCRETE_VALUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 May 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CDF_DISCRETE_TEST\n' );
  fprintf ( 1, '  CDF_DISCRETE_VALUE evaluates the CDF associated with a\n' );
  fprintf ( 1, '  discrete histogram.\n' );
%
%  Set up the discrete histogram from sample data.
%
  s_num = 5;
  s_min = 0.0;
  s_max = 10.0;
  s = [ 2.0, 2.0, 4.0, 5.0, 9.0 ];
  fprintf ( 1, '\n' );
  fprintf ( 1, '  S_MIN = %g\n', s_min );
  fprintf ( 1, '  S_MAX = %g\n', s_max );
  r8vec_print ( s_num, s, '  Sample data:' );

  [ x_num, x, y ] = setup_discrete ( s_num, s, s_min, s_max );

  r8vec2_print ( x_num, x, y, '  Discrete histogram data:' );
%
%  Evaluate the discrete CDF.
%
  v_num = 21;
  v = linspace ( s_min, s_max, v_num );
  cdf = cdf_discrete_value ( x_num, x, y, v_num, v );
  r8vec2_print ( v_num, v, cdf, '  Discrete CDF table:' );
%
%  Plot the discrete CDF.
%
  plot ( v, cdf, 'bo-', 'Linewidth', 3 )
  grid on
  xlabel ( '<--- X --->' )
  ylabel ( '<--- CDF(X) --->' )
  title ( 'Discrete CDF function' )

  filename = 'cdf_discrete_test.png';
  print ( '-dpng', filename )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Discrete CDF plotted as "%s".\n', filename );

  return
end

repo: jburkardt-m


text:
function proj = DetectorPointScatterCorrection(proj, geo, ScCalib, gpuids)
%% Detector Point Scatter Correction
% Reference: Improved scatter correction using adaptive scatter kernel superposition
% Date: 2021-03-26
% Author: Yi Du (yi.du@hotmail.com)

%% Empirical values from reference paper
% unit: cm-2
% a0 = 3.43; (in refrence paper, but incorrect)
a0 = 1;

a1 = 0.000309703536035;
a1 = str2double(ScCalib.CalibrationResults.Globals.DetectorScatterModel.PScFit0.Text);

% unit: cm-1
a2 = 0.546566915157327;
a2 = str2double(ScCalib.CalibrationResults.Globals.DetectorScatterModel.PScFit1.Text);

% unit: cm-1
a3 = 0.311272841141691;
a3 = str2double(ScCalib.CalibrationResults.Globals.DetectorScatterModel.PScFit2.Text);

% unit: cm-1
a4 = 0.002472148007134;
a4 = str2double(ScCalib.CalibrationResults.Globals.DetectorScatterModel.PScFit3.Text);

a5 = -12.6606856375944;
a5 = str2double(ScCalib.CalibrationResults.Globals.DetectorScatterModel.PScFit4.Text);

% for amplitude normalization
CoverSPR = 0.04;

% unit: mm
offset=geo.offDetector;

% grid unit: mm
us = ((-geo.nDetector(1)/2+0.5):1:(geo.nDetector(1)/2-0.5))*geo.dDetector(1); % + offset(1);
vs = ((-geo.nDetector(2)/2+0.5):1:(geo.nDetector(2)/2-0.5))*geo.dDetector(2); % + offset(2);
% unit mm - > cm
% unit converter: 1/10 for mm-> cm
mm2cm = 1/10;
us = us * mm2cm;
vs = vs * mm2cm;

%% Downsampling
% about 10 mm in axial direction
dus = downsample(us, 26);
% about 4 mm in axial direction
dvs = downsample(vs, 10);

ds_rate = 8;
dus = decimate(us, ds_rate);
dvs = decimate(vs, ds_rate);

%% Grid mesh
[uu,vv] = meshgrid(us,vs); %detector
[duu, dvv] = meshgrid(dus,dvs); %detector

%% Scatter convolution kernel
grid = sqrt(duu.^2 + dvv.^2);
% a0 is the normalization factor
hd = a0*( a1* exp(-a2 * grid) + a3 * (exp( -a4 * ( grid - a5).^3 )));

% a0 = CoverSPR/sum(hd(:));

% normalized to 4% SPR coverage
hd = CoverSPR/sum(hd(:)) .* hd;

%% GPU based
% reset(gpuDevice(1));
for ii=1:length(gpuids)
    g = gpuDevice(gpuids.devices(ii));
    reset(g);
end
gproj = gpuArray(single(proj));

%% 2D Convolution with downsampling and upsampling
for ii = 1:size(proj,3)
    % CPU version
    %{
    page = interp2(uu, vv, proj(:,:,ii), duu, dvv);
    % gross scatter distribution
    sc = conv2(page, hd, 'same');
    % upsample the scatter distribution to the same grid level as the
    % measured intensity
    scpage = interp2(duu, dvv, sc, uu, vv, 'spline');
    % primary = measure - scatter
    proj(:,:,ii) = proj(:,:,ii) - scpage;
    %}

    % GPU version
    page = interp2(uu, vv, gproj(:,:,ii), duu, dvv);
    % gross scatter distribution
    sc = gather(conv2(page, hd, 'same'));
    % upsample the scatter distribution to the same grid level as the
    % measured intensity
    scpage = interp2(duu, dvv, sc, uu, vv, 'spline');
    % primary = measure - scatter
    gproj(:,:,ii) = gproj(:,:,ii) - scpage;
end

proj = single(gather(gproj));
% Reset GPU
for ii=1:length(gpuids)
    g = gpuDevice(gpuids.devices(ii));
    reset(g);
end

%% Cutoff for over-correction
proj(proj<0) = NaN;
for ii = 1:size(proj,3)
    % default fast extrapolation: robust for noise and holes at boundaries
    proj(:,:,ii) = single(inpaint_nans(double(proj(:,:,ii)), 2));
end
proj(proj<0) = eps;

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function tetrahedron_arbq_rule_test ( )

%*****************************************************************************80
%
%% TETRAHEDRON_ARBQ_RULE_TEST tests the TETRAHEDRON_ARBQ_RULE library.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    10 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TETRAHEDRON_ARBQ_RULE_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TETRAHEDRON_ARBQ_RULE library.\n' );

  degree = 8;
  n = tetrahedron_arbq_size ( degree );
  header = 'tetrahedron08';

  tetrahedron_arbq_rule_test01 ( degree, n );

  tetrahedron_arbq_rule_test02 ( degree, n, header );

  tetrahedron_arbq_rule_test03 ( degree, n, header );

  tetrahedron_arbq_rule_test04 ( degree, n );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TETRAHEDRON_ARBQ_RULE_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


function Omega = generateQ(m)

% generate an orthonormal matrix.

G     = randn(m);
[Q,R] = qr(G);
% normalize to positive entry in the diagonal
In    = diag(sign(diag(R)));
Omega = Q  * In;


text:
function aedat = ReorientatePolarity(aedat, transpose, flipX, flipY)

%{
For the frame data, if present, apply transpose, flipX and flipY operations
in that order according to the boolean inputs.
%}

dbstop if error

if ~isfield(aedat, 'data')
    disp('No data found')
    return
end
if ~isfield(aedat.data, 'polarity')
    disp('No polarity data to reorientate')
    return
end
if transpose
    temp                    = aedat.data.polarity.x;
    aedat.data.polarity.x   = aedat.data.polarity.y;
    aedat.data.polarity.y   = temp;
end

if flipX
    aedat.data.polarity.x   = aedat.info.deviceAddressSpace(1) - aedat.data.polarity.x - 1;
end

if flipY
    aedat.data.polarity.y   = aedat.info.deviceAddressSpace(2) - aedat.data.polarity.y - 1;
end

author: panpanfei


text:
function RND = mnbinrnd(W, K, M, N)
%MNBINRND       Maszle's negative binomial random numbers
%
%   RND = MNBINRND(W, K, M, N)
%
%   returns an M x N matrix of random numbers chosen from a negative
%   binomial distribution with mean W and aggregation K.
%
%   This is an alternate form of the negative binomial commonly
%   employed in biology to describe aggregated count data.
%   W and K are related to the traditional parametrization by the
%   relationships:
%
%       P = K/(K + W),  and K = R,
%
%   with the distinction that K is any positive real number, not just
%   positive integers.  K < 1 yields a highly skewed distribution.
%
%   The size of RND is the common size of W and K if both are matrices.
%   If either parameter is a scalar, the size of RND is the size of
%   the other parameter.
%
%   Alternatively, RND = NBINRND(W, K, M, N) returns an M x N matrix.
%
%   Requires stats toolbox.
%
%   See also MNBINPDF

%----------------------------------------------------------------------
% Copyright (c) 1999.  Don R. Maszle.  All rights reserved.
%
%   -- Revisions -----
%        Date:  4 March 1999
%      Author:  Don R. Maszle
%      E-mail:  maze@sparky.berkeley.edu
%   -- SCCS  ---------
%----------------------------------------------------------------------

%-- Perform standard argument checks for random generators

if (nargin < 2)
    error('Requires at least two parameters.');
end

if (nargin == 2)
    [iError nRows nCols] = rndcheck(2,2,K,W);
    if (max(size(K)) == 1)
      K = K(ones(nRows,1),ones(nCols,1));
    end
    if (max(size(W)) == 1)
      W = W(ones(nRows,1),ones(nCols,1));
    end
end

if (nargin == 3)
    [iError nRows nCols] = rndcheck(3,2,K,W,M);
    K = K(ones(M(1),1),ones(M(2),1));
    W = W(ones(M(1),1),ones(M(2),1));

end

if (nargin == 4)
    [iError nRows nCols] = rndcheck(4,2,K,W,M,N);
    K = K(ones(M,1), ones(N,1));
    W = W(ones(M,1), ones(N,1));
end

if (iError > 0)
    error('Size information is inconsistent.');
end

% From "Non-uniform random variate generation", Luc Devroye, New York,
% Springer-Verlag, 1986.  (Thanks to Charles N. Haas, Drexel U. for
% pointing me to this source.)
%
% As Devroye derives, the negative binomial can be expressed as a
% Poisson of a Gamma with parameters K and (1-P)/P, where P = K/(K + W).

RND = poissrnd(gamrnd(K, W./K, M, N));

author: Sable

repo: mcbench-benchmarks